Starlight-polarization-based tomography of the magnetized ISM: PASIPHAE’s line-of-sight inversion method

Astronomy

&

Astrophysics

A&A 670, A164 (2023)

https://doi.org/10.1051/0004-6361/202244625

Starlight-polarization-based tomography of the magnetized ISM:

ASIPHAE

’s line-of-sight inversion method

V. Pelgrims

, G. V. Panopoulou

, K. Tassis

, V. Pavlidou

, A. Basyrov

, D. Blinov

, E. Gjerløw

S. Kiehlmann

, N. Mandarakas

, A. Papadaki

, R. Skalidis

, A. Tsouros

, R. M. Anche

H. K. Eriksen

, T. Ghosh

, J. A. Kypriotakis

, S. Maharana

, E. Ntormousi

, T. J. Pearson

S. B. Potter

, A. N. Ramaprakash

, A. C. S. Readhead

, and I. K. Wehus

Institute of Astrophysics, Foundation for Research and Technology-Hellas, N. Plastira 100, Vassilika Vouton, 71110 Heraklion,

Greece

Department of Physics, and Institute for Theoretical and Computational Physics, University of Crete, Voutes University campus,

70013 Heraklion, Greece

e-mail:

pelgrims@physics.uoc.gr

Hubble Fellow, California Institute of Technology, MC350-17, 1200 East California Boulevard, Pasadena, CA 91125, USA

Institute of Theoretical Astrophysics, University of Oslo, PO Box 1029 Blindern, 0315 Oslo, Norway

Institute of Computer Science, Foundation for Research and Technology-Hellas, 71110 Heraklion, Greece

Department of Astronomy/Steward Observatory, Tucson, AZ 85721-0065, USA

Inter-University Centre for Astronomy and Astrophysics, Post bag 4, Ganeshkhind, Pune 411007, India

School of Physical Sciences, National Institute of Science Education and Research, HBNI, Jatni 752050, Odisha, India

Scuola Normale Superiore di Pisa, piazza dei Cavalieri 7, 56126 Pisa, Italy

Cahill Center for Astronomy and Astrophysics, California Institute of Technology, 1216 E California Blvd, Pasadena, CA 91125,

USA

Department of Physics, University of Johannesburg, PO Box 524, Auckland Park 2006, South Africa

South African Astronomical Observatory, PO Box 9, Observatory, 7935 Cape Town, South Africa

Received 29 July 2022 / Accepted 28 December 2022

ABSTRACT

We present the first Bayesian method for tomographic decomposition of the plane-of-sky orientation of the magnetic field with the use

of stellar polarimetry and distance. This standalone tomographic inversion method presents an important step forward in reconstructing

the magnetized interstellar medium (ISM) in three dimensions within dusty regions. We develop a model in which the polarization

signal from the magnetized and dusty ISM is described by thin layers at various distances, a working assumption which should be

satisfied in small-angular circular apertures. Our modeling makes it possible to infer the mean polarization (amplitude and orientation)

induced by individual dusty clouds and to account for the turbulence-induced scatter in a generic way. We present a likelihood function

that explicitly accounts for uncertainties in polarization and parallax. We develop a framework for reconstructing the magnetized ISM

through the maximization of the log-likelihood using a nested sampling method. We test our Bayesian inversion method on mock data,

representative of the high Galactic latitude sky, taking into account realistic uncertainties from

Gaia

and as expected for the optical

polarization survey P

ASIPHAE

according to the currently planned observing strategy. We demonstrate that our method is effective

at recovering the cloud properties as soon as the polarization induced by a cloud to its background stars is higher than

∼

for

the adopted survey exposure time and level of systematic uncertainty. The larger the induced polarization is, the better the method’s

performance, and the lower the number of required stars. Our method makes it possible to recover not only the mean polarization

properties but also to characterize the intrinsic scatter, thus creating new ways to characterize ISM turbulence and the magnetic

field strength. Finally, we apply our method to an existing data set of starlight polarization with known line-of-sight decomposition,

demonstrating agreement with previous results and an improved quantification of uncertainties in cloud properties.

Key words.

dust, extinction – ISM: magnetic fields – ISM: structure – polarization – methods: statistical

1. Introduction

Studies of the interstellar medium (ISM) have relied on two-

dimensional (2D) projections on the sky until recently. With the

advent of sophisticated techniques and state-of-the-art facilities,

astronomy has entered a new realm in which the third dimension

can finally be accessed with accuracy, enabling the mapping of

the ISM in three dimensions (3D). Astronomers – and the pub-

lic – will soon be able to experience the Universe in 3D flying

through real-world data sets loaded in virtual environments.

Gaia

data on stellar distances in particular (e.g., Gaia

Collaboration 2016; Gaia Collaboration 2021; Bailer-Jones et al.

2021) have allowed for the precise localization in 3D space of

more than one billion stars in our Galaxy through the accu-

rate determination of stellar parallaxes. Coupling measurements

of stellar parallaxes to reddening, Bayesian inversion methods

have already been successful at reconstructing 3D maps of the

dust density distribution (e.g., Green et al. 2019; Lallement et al.

2019; Leike & Enßlin 2019; Leike et al. 2020), leading to stun-

ning 3D images mapping the dust content of the ISM, in the

A164, page 1 of 30

Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0),

which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

A&A 670, A164 (2023)

Solar neighborhood from within the first tens of parsec and up to

much larger distances within the Galactic disk, already covering

a substantial fraction of the Galaxy (

6000

800

for

Lallement et al. 2019, 2022). Such 3D mappings of the ISM con-

tent are of great interest for several areas in astrophysics. They

shed new light on the dynamics shaping the Galaxy, breaking

degeneracies caused by 2D mapping on the 3D shapes of ISM

clouds and cloud complexes, their formation mechanisms, and

their history (e.g., Ivanova et al. 2021; Bialy et al. 2021; Rezaei

Kh. & Kainulainen 2022; Zucker et al. 2022; Lallement et al.

2022; Konstantinou et al. 2022). Ultimately, 3D images of the

dust content of the Galaxy could also help in the characteriza-

tion of Galactic foregrounds for observational cosmology and

extra-galactic astrophysics (e.g., Martínez-Solaeche et al. 2018).

Impressive as they may be, 3D reconstructions of the ISM

dust distribution are leaving out an important component of

the Galaxy: magnetic fields, which are ubiquitous in the ISM.

Magnetic fields are relevant in a variety of processes, from regu-

lating star formation (e.g., Mouschovias et al. 2006; Hennebelle

& Inutsuka 2019; Li 2021) to shaping large-scale structures in

the disk and the halo of the Galaxy (e.g., Beck 2015). Mag-

netic fields in the Galaxy also affect our ability to study the

Universe’s structure and history, from its first moments to its

later ages. Aspherical dust grains line up their shortest axis with

the ambient magnetic field (e.g., Andersson et al. 2015). As a

result, the thermal radiation emitted by those grains is polarized.

This emission constitutes the major limitation in cosmologists’

search for primordial B modes, the clear proof of primordial

gravitational waves from inflation, and cosmic birefringence in

the polarization of the cosmic microwave background (CMB;

e.g., BICEP2/Keck Collaboration & Planck Collaboration 2015;

Planck Collaboration XXX 2016; Planck Collaboration XI 2020;

Diego-Palazuelos et al. 2022). This emission also represents a

foreground in polarization studies of individual extra-galactic

objects (e.g., Pelgrims 2019).

Significant effort has been invested in the last two decades

to characterize the dust-polarized emission in order to disen-

tangle it from the cosmological signal. However, this task has

been proven to be very convoluted. Variations in dust spectral

emission distribution, either in the plane of the sky (POS) or

along the line of sight (LOS; Tassis & Pavlidou 2015; Planck

Collaboration Int. L 2017; Pelgrims et al. 2021; Ritacco et al.

2023), and unexpected signatures of the dust signal in polar-

ization power spectra (e.g., Planck Collaboration Int. XXXVIII

2016) – all rooted in the tight connection between dust clouds

and the magnetic field – add many layers of complexity. Various

sophisticated techniques are being developed to address these

problems. The most direct way of attacking them and of provid-

ing confident and accurate solutions requires 3D mapping of the

Galactic magnetic field in dusty regions (e.g., Hervías-Caimapo

& Huffenberger 2022; Pelgrims et al. 2022; Konstantinou et al.

2022; Huang 2022).

Accessing the LOS structure of the magnetic field from dust

emission alone is not feasible. Three-dimensional maps of the

dust distribution can help identify LOSs with several clouds

and place constraints on their respective significance; however,

those maps alone provide no information about the magnetic

fields permeating those clouds. While they can be combined

with maps of dust-polarized emission in a coherent analysis

to model the Galactic magnetic field (GMF) on large scales

(Pelgrims et al. 2020), they cannot provide significant informa-

tion on cloud scales, with perhaps some exceptions (e.g., Rezaei

Kh. et al. 2020).

Fortunately, there are other probes that make it possible

to infer the structure of the magnetized ISM in 3D. Among

those, the linear polarization of stars, measured from the near-

infrared (NIR) to the near-ultraviolet, is of particular interest,

and can be used to study and model the dusty and magnetized

ISM, from the smallest to the largest scales. While starlight

usually starts out unpolarized, the same aspherical dust grains

that are responsible for the polarized thermal emission induce

a polarization to it when partially absorbing it, due to dichroic

extinction (e.g., Andersson et al. 2015). Starlight polarization has

been related to the magnetic field and the ISM in the Galaxy

since its early observation (e.g., Hiltner 1949, 1951; Davis &

Greenstein 1951; Heiles 2000). In comparison to other probes

of the magnetized ISM, starlight polarization has the significant

advantage that it can provide direct 3D information as soon as

stellar distances are known. The potential of such 3D magnetic

tomography to recover information on the LOS structure of the

magnetic field has been demonstrated recently by Panopoulou

et al. (2019b) using data collected from the RoboPol polarimeter

(Ramaprakash et al. 2019), while correlation analysis of dust-

polarized emission at sub-millimeter wavelengths and starlight

polarization data has proven useful to locate the distance to

the dominant polarized dust emission component seen at high

Galactic latitude (Skalidis & Pelgrims 2019).

In recent years, several regions of the sky have been mapped

with a high density of stellar polarization measurements (

100

stars per sq.degree), including a significant portion of the inner

Galaxy (in the NIR Clemens et al. 2020), as well as more diffuse

regions of the ISM (e.g., Panopoulou et al. 2019b; Skalidis et al.

2022). These data sets have paved the way to 3D mapping mag-

netic fields in the general ISM of the Galaxy (e.g., Pavel et al.

2012), far from the dense regions of star formation that had been

traditionally studied with large stellar samples (e.g., Pereyra &

Magalhães 2004; Sugitani et al. 2011; Marchwinski et al. 2012;

Santos et al. 2014; Franco & Alves 2015; Kwon et al. 2015;

Eswaraiah et al. 2017). In the near future, the P

ASIPHAE

survey

(Tassis et al. 2018) and the SOUTH POL survey (Magalhães et al.

2005) will enable a leap forward by generating stellar polariza-

tion data for millions of stars, covering a large fraction of the sky.

In conjunction with measurements of stellar distances obtained

Gaia

, those data sets will pave the way for the characterization

of the dusty magnetized ISM in 3D. Since stellar polarization

traces the very same medium (magnetized dust) that produces

the dominant CMB polarization foreground, starlight polariza-

tion data may offer a unique independent means to model out

the polarization signal of the Galaxy, allowing the study of the

very first moments of the Universe.

The observed polarization of each single star is the inte-

grated effect of dichroic absorption from all interstellar clouds

lying between us and the star. For this reason, in order to derive

the complex 3D structure of the magnetized ISM from starlight

polarization data and stellar distances, we need to develop meth-

ods that invert this LOS integration. So far, no standard method

has been established in the literature to accomplish such a task in

an automated, Bayesian way. On the one hand, different ad hoc

methods have been considered (e.g., Andersson & Potter 2006;

Panopoulou et al. 2019b; Doi et al. 2021), but they are not eas-

ily scalable to large data sets since they are not well adapted for

automation and they do not allow for a straightforward, robust

estimation of the credible interval of the reconstruction. On the

other hand, methods developed for extinction data cannot be

used unaltered on starlight polarization data. The main reason is

that polarization is a pseudo-vector quantity. This implies that

A164, page 2 of 30

V. Pelgrims et al.: LOS tomography of the dust polarization sky

it cannot be described by a single scalar quantity – two are

needed: either the degree of polarization and polarization posi-

tion angle, or its linear Stokes parameters. Additionally, even if

contributions from individual clouds are additive (as is the case

for linear Stokes parameters in the case of low amounts of extinc-

tion), polarization increments can be either positive, negative,

or null. Because of these fundamental differences, dedicated,

specialized methods need to be developed for the problem of

starlight polarization tomography.

In this paper, we present such a specialized Bayesian method,

developed for the P

ASIPHAE

survey, implemented in Python, and

now made publicly available for use by the community

. The

inversion method developed here works on a per line-of-sight

basis. We defer to future work for information on the extension

of the method, which must take the correlation of the solutions

in the plane of the sky into account.

In Sect. 2, we present our model for the distance depen-

dence of starlight polarization along sightlines of the diffuse

ISM, and we explain how we built our data equation and derived

the corresponding likelihood. In Sect. 3, we provide details on

the implementation of our Bayesian method and validate its per-

formance by applying it to two simple examples of mock data.

In Sect. 4, we present extensive testing of the performance of the

method in the low signal-to-noise (S/N) regime and identify the

method’s limitations. We apply the method to currently available

data in Sect. 5 and compare the results from our method to the

literature. We finally summarize and discuss our work in Sect. 6.

This paper contains two appendices. In Appendix A, we explain

the creation of the mock observations used for performance test-

ing, which are based on actual star samples, realistic estimates

of the uncertainties on stellar parallax and polarization, and that

rely on a complete toy model for the 3D structure of the ISM

along sightlines. In Appendix B, we explore our toy model of

the magnetic field geometry to gain intuition on the effects of

turbulence-induced fluctuations in the ISM on the polarization

observables.

2. Model, data equation, and likelihood

In this section, we lay the foundation for a model that describes

the distance-dependence of stellar polarization toward a sightline

of the diffuse ISM. We construct a generic data equation and

build the corresponding likelihood that relates model parameters

and star data. We first discuss the case of a single cloud along

the LOS and then proceed to the generalization to cases with

multiple-clouds.

2.1. Model: Thin-layer magnetized clouds

We model the LOS polarization induced by an individual cloud

to background stars as being dominated by a single thin polar-

izing dust layer at the cloud distance (

). As already described

by many authors (e.g., Andersson et al. 2015; Hensley & Draine

2021; Draine & Hensley 2021), the polarization induced by a

dust cloud to the light of background stars depends on the dust

opacity at the frequency of observation, the polarizing efficiency

(which relates dust reddening

(

−

)

to a maximum polariza-

tion fraction) and on the apparent 3D orientation of the magnetic

field (

) that permeates the cloud. The latter is described through

the inclination angle of the magnetic field lines with respect to

the POS (

) as well as the position angle of the POS compo-

nent of the magnetic field (

). For a single star behind a cloud,

https://github.com/vpelgrims/Bisp_1/

with starlight intensity

, the vector of its relative Stokes param-

eters in the visible

(

)

(

)

equals the cloud

polarization vector

(

)

given by

max

cos

cos[2

]

sin[2

]

(1)

In this equation, in which we have neglected any possible source

of noise,

max

≈

13%

(

−

)

, (Panopoulou et al. 2019a; Planck

Collaboration XII 2020) where the reddening

(

−

)

generally

depends on the dust grain physical properties and on the col-

umn density. Using single-frequency starlight polarization only,

we have access to the position angle of the POS component of

the magnetic field (related to the electric vector position angle,

EVPA, of the stellar polarization) and to the magnitude of the

induced polarization, that is the degree of polarization:

(

)

(related to the degree of stellar polarization). The

latter is affected by the dust extinction, the dust polarizing effi-

ciency, the inclination of the magnetic field with respect to the

POS, and by possible LOS depolarization caused by turbulence

within the cloud.

If the distribution of dust and the magnetic field properties

were spatially homogeneous within a cloud, a single stellar mea-

surement would suffice to describe the polarization properties

it induces. Considering an ensemble of stars to constrain the

cloud polarization properties makes it possible to take advan-

tage of the number statistics and to sample the distance axis to

provide constraints on the cloud distance. The former is criti-

cal for the S/N regime that is expected at high Galactic latitudes

(Skalidis et al. 2018). Furthermore, observations of interstellar

dust reveal fluctuations in the dust distribution on a range of

scales (e.g., Miville-Deschênes et al. 2016) and fluctuations in

the density and the magnetic field are also expected as a result

of magneto-hydrodynamic (MHD) turbulence (e.g., Goldreich &

Sridhar 1995; Cho & Lazarian 2003; Heiles & Crutcher 2005).

Therefore, in order to obtain a realistic description of the polar-

ization properties of a cloud within a region, we consider an

ensemble of stellar measurements from LOSs within a finite cir-

cular aperture (called “beam” in the remainder of the paper)

toward the cloud. We describe the Stokes parameters induced by

the cloud to the ensemble of stars with a well-defined mean and

a measure of dispersion about that mean. We refer to this disper-

sion as intrinsic scatter, to distinguish it from other sources of

dispersion in the measurements, such as noise.

The intrinsic scatter has effects on polarization observables.

We explore these effects in Appendix B where we characterize

the variations produced by the intrinsic scatter on

and on

;

or more generally in the

(

)

plane. In particular, we show

that 3D variations of the magnetic field can generate biases and a

nonzero cross term in the polarization plane. The biases are irrel-

evant in our case since we are interested in recovering the mean

values. The cross term, however, needs to be accounted for given

that it might reach a non-negligible fraction of the variance of the

Stokes parameters

(

; Montier et al. 2015). Other sources

of variance in the polarization properties, such as fluctuations of

the dust extinction across the sky, may reduce the importance

of the off-diagonal element compared to the diagonal elements.

Nevertheless, we retain the off-diagonal terms in our analysis for

completeness.

As a first approximation, we assume that the turbulence-

induced variations generate a bivariate normal distribution about

the mean in the

(

)

plane. As a result, and in absence of

observational noise, our stochastic model for the vector of Stokes

A164, page 3 of 30

A&A 670, A164 (2023)

parameters of a star

in the background of a cloud is

int

)

(2)

where

and

now denote the mean polarization values

induced by the cloud. We denote the mean Stokes vector of the

cloud as

(

)

†

int

)

is a random realization of

a 2D bivariate normal distribution centered on

with the

2-by-2 covariance matrix,

int

. The latter encodes the variances

and covariances induced by sources of intrinsic fluctuations (e.g.,

turbulence) on the Stokes parameters. According to this generic

description, six parameters are necessary to describe the polar-

ization data of stars toward such a cloud: the distance of the cloud

(

), the mean Stokes parameters

(

)

and three numbers to

characterize the intrinsic-scatter covariance matrix

int

2.2. Data equation

To write the data equation, we need to account for the fact that

a star at distance

may either be in the foreground (

)

or in the background (

) of the cloud. In the former case

no polarization is induced by the cloud and in the latter case the

star polarization will be given by the mean polarization of the

cloud plus one random realization of the intrinsic scatter. This

piecewise-constant behavior is implemented through the use of

the Heaviside function (

(

)

otherwise).

We further add a noise term (

) to our stochastic model for

the Stokes parameters. We consider that the observational noise

which results from photon noise, instrumental polarization, and

on-sky instrumental calibration, is described by a bivariate nor-

mal distribution

obs

)

with covariance matrix,

obs

, where

the off-diagonal terms can be nonzero. Unlike the intrinsic scat-

ter, the covariance matrix corresponding to the observational

uncertainties is generally source dependent; it might depend on

the source’s brightness, for example. The variance and covari-

ance in the

(

)

plane result from both the intrinsic scatter

and the observational uncertainties. We consider them as inde-

pendent sources of Gaussian scatter in the polarization plane.

Therefore, for stars in the background of a cloud the covari-

ance matrices (

int

and

obs

) are summed. For stars that are not

background to a cloud, only the observational uncertainties are

relevant. The total covariance matrix thus takes the form:

obs

int

(

−

)

(3)

As a result, the data equation for the case of a single cloud along

the LOS is:

(

−

)

(

int

obs

)

obs

)

otherwise,

(4)

where

is the vector of the measured Stokes parameters,

is the

distance of the star, and

is the source-dependent noise term.

Denoting

the mean polarization induced by any dust cloud

between us and the star and using Eq. (3), we thus write the

vector of measured Stokes parameters as a random draw of a

bivariate normal distribution with mean

(with value either

) and 2-by-2 covariance matrix

←

(

)

exp

−

†

−

(5)

in which

det

(

)

is the determinant of

and where

we have introduced

−

with

depend-

ing on whether the star is in the foreground of background to

the cloud. In Eqs. (3) to (5), both the modeled mean and the

modeled covariance terms for the stellar polarization depend on

whether the source is in the foreground or the background of

the cloud.

2.3. Likelihood

The data equation (Eq. (4)) concerns only the Stokes parame-

ters of the stars. As explicitly written, our model for the Stokes

parameters relates to the distance of the stars which, in turn,

is also a measured quantity that comes with an uncertainty.

This source of uncertainty adds complexity to the problem that

we wish to solve, as distance uncertainties might impact the

model prediction for

, in particular for those stars that are

near a cloud. Star distances are nowadays known with great

accuracy (e.g., Bailer-Jones et al. 2021) through their parallaxes

(

) which, to good approximation, have Gaussian uncertainties

(

). For this reason we choose to work in terms of paral-

laxes rather than distances, the two being related through the

inverse relation

where parallaxes are measured in arc-

seconds and distances in parsec. To account for this extra source

of complexity, we notice that the star distance entering the data

equation above should be the true distance of the star, and there-

fore we use the true parallax of the star in the following. We

modify the argument of the Heaviside function from

(

−

)

(

−

)

, in which

denotes the true parallax of a star.

We consider that the measured parallax (

) is a random realiza-

tion of a Gaussian distribution centered on the true parallax with

uncertainty

←

,σ

)

√

πσ

exp







−

(

−

)







(6)

In this work, we assume that the measurements of the paral-

lax and of the optical polarization of stars are independent and

uncorrelated. Further, we assume that the Stokes parameters for

star polarization are functions of the true parallax of the star

through the generic data equation built in the previous section

(Eq. (4)). With these notations, the likelihood of the observa-

tional data point for star

with measured parallax

and Stokes’

vector

takes the form:

int

, π

,σ

obs

,σ

int

obs

√

πσ

−

(

−

)

−

†

−

(7)

where, in the last line, we explicitly write the parallax likeli-

hood as a 1D Gaussian with standard deviation equal to the

observational uncertainty and the polarization likelihood as a 2D

Gaussian with the total covariance matrix that accounts for both

the observational uncertainties and the contribution to intrinsic

scatter from the crossed cloud (Eq. (3)).

We are not interested in modeling the true parallax of the

star (

). Instead we wish to marginalize over it to define the

likelihood of the cloud parameters given an observation for star

. This marginalization allows us to separate the likelihood into

two parts, one corresponding to the background of the cloud and

A164, page 4 of 30

V. Pelgrims et al.: LOS tomography of the dust polarization sky

the other to its foreground as follows:

(

int

)

∞

int

, σ

obs

int

, σ

obs

}

background

∞

, σ

obs

}

foreground

int

obs

, σ

obs

∞

, σ

int

obs







erf







−

√













obs







−

erf







−

√













(8)

For a given sample of stars with polarization measurements

and with known parallaxes and uncertainties, and under the

assumption that the data are independent, the likelihood of

the cloud parameters for a given LOS is given by the product

of the likelihoods of the data points:

(

int

)

star

(

int

)

(9)

This is the total likelihood function that we need to maximize to

constrain our model parameters given the data.

2.4. Multicloud case

The generalization of the single-cloud model to the case with

multiple independent clouds along the LOS is straightforward

(we take

to be the number of clouds along the LOS). We

consider that the Stokes parameters of a star in the background

of multiple clouds result from the addition of the contributions

from individual clouds. This approximation, which is correct in

the low polarization regime (e.g., Appendix B of Patat et al.

2010), is well motivated for translucent LOSs through the diffuse

ISM given that dust clouds can be considered as weak polariz-

ers. The validity of this approximation may need to be revised

for denser regions of the ISM, such as the Galactic plane or other

LOSs through very dense molecular clouds.

In this work, we assume the linearity of the polarization sig-

nal and defer to future work the addition of more complex cases.

Hence, in the low-polarization regime, the mean polarization

vectors (

[

]

) of clouds along the LOS are additive. The same is

true for the covariance matrices from the intrinsic scatter. Here,

we have introduced the superscript

[

]

to label clouds from

[1]

[

]

, for the nearest and the farthest cloud, respectively. As for

the case of a single cloud, the total covariance (

) on the Stokes

parameters for a star

depends on the observational uncertain-

ties and on the sum of all sources of intrinsic scatter intervening

along the LOS, from the star to the observer. Thus, assuming

independent clouds along the LOS, the data equation for a

star

is:

(

)

(10)

where

≤

[

]

and

obs

≤

[

]

int

(11)

in which we implicitly assume that star

lies behind cloud

[

]

and in front of cloud

[

, if

and behind all clouds if

To determine the likelihood for the set of cloud parame-

ters given star data we marginalize over the true parallax of the

star. This step allows us to separate the likelihood into

terms, building the likelihood of a mixture model where the dif-

ferent terms correspond to a different (increasing) number of

foreground clouds:

{

[

]

[

]

[

]

int

}

obs

[1]

obs

[1]

int

...







[

]

obs

[

]

int







(12)

The coefficients

result from the integration of the probability

density function of the star parallax in the inter-cloud ranges of

parallax and take the form







erf







−

√







−

erf







−

√













(13)

Finally, from an ensemble of stars, the likelihood of the cloud

parameters for a given LOS is given by the product of the

likelihoods of the data points, as in Eq. (9).

3. The Bayesian inversion method: Implementation

and validation

In this section we first describe the implementation of our

maximum-likelihood method to decompose the polarization

properties of clouds along the LOS using stellar polarization

and distance. We validate our method using mock stellar data

sets based on a toy model of discrete clouds along the LOS and

explore the sensitivity of the log-likelihood with respect to the

different model parameters. We then provide a solution on how,

facing real observations, we can select the correct model, that

is, to choose the correct number of clouds that exist along a

sightline.

3.1. Implementation

To maximize the likelihood function and estimate the posterior

distributions of our model parameters we rely on a numerical

method. We choose to use the code

dynesty

(Speagle 2020) to

sample the parameter space using the nested sampling method

(Skilling 2004). This code has already proven to be powerful and

reliable in solving astrophysical problems similar to ours (e.g.,

A164, page 5 of 30

A&A 670, A164 (2023)

Zucker et al. 2019, 2022; Alves et al. 2020). The algorithm uses

sampling points (named ‘live points’ in

dynesty

’s definition)

to explore the parameter space in a dynamic nested sampling

scheme and estimate the posterior distributions on model param-

eters. It has two main advantages compared to other sampling

methods: first, it returns an estimate of the model evidence and

second, it includes a well-defined stopping criterion, suitable for

automation of the fitting process.

The code

dynesty

takes as input the function of the log-

likelihood that has to be maximized and the definition of

functions that implement our prior knowledge on the model

parameters. We implemented both uniform and Gaussian pri-

ors for the cloud parallaxes and cloud mean polarization. The

Gaussian priors are defined through their means and stan-

dard deviations while uniform priors are defined by their lower

and upper limits. We strongly recommend using flat priors

for the element of the covariance matrix encoding the effects

from turbulence. This is because the diagonal elements of the

matrix must be positive and can approach zero depending on

the (unknown) orientation of the magnetic field in 3D (see

Appendix B). In this case, the diagonal elements (

int

and

int

) are first proposed independently in their respective ranges

and then the off-diagonal elements are drawn such that the semi-

positive-definiteness of the covariance matrix is guaranteed, that

is,

int

is drawn from a uniform distribution in the range

(

−

int

)

; excluding the limits. In the

case of multiple clouds, the prior function makes internally the

distinction between clouds, ranked by their distances. We have

hard-coded a lower limit on the number of stars that can exist

between two clouds. We fix this limit at five.

Given a set of data points

{

(

)

}

and corre-

sponding uncertainties

(

obs

)

, the log-likelihood function

that

dynesty

has to maximize is given by

log

{

[

]

[

]

[

]

int

}

log







star

{

[

]

[

]

[

]

int

}







star

log

{

[

]

[

]

[

]

int

}

(14)

which requires the number of clouds as an entry and where the

’s are given by Eq. (12). The number of clouds populating the

LOS is a priori unknown. We discuss in Sect. 3.5 how we intend

to determine it. We implemented the likelihood functions for up

to five clouds along the LOS. Even though the generalization of

the implementation to higher number of clouds is trivial we do

not deem it necessary given that only few LOSs at intermediate

and high Galactic latitudes are expected to show a large number

of components with significant contribution to the polarization

signal (Panopoulou & Lenz 2020).

3.2. Mock data for two example LOSs

To test our method we developed a simple but complete imple-

mentation of our layer model to generate mock data sets with

realistic number of stars, stellar distance and brightness distribu-

tions, and uncertainties both on parallax and polarization. This

implementation, which includes a self-consistent prescription for

the intrinsic scatter, is presented in detail in Appendix A.3. Our

toy model for generating stellar polarization has five free param-

eters per cloud: the cloud parallax (

)

, the maximum

degree of polarization (

max

), the inclination (

reg

) and posi-

tion (

reg

) angles of

reg

and, finally, the relative amplitude of

fluctuations in magnetic field orientation (

turb

). We draw the

reader’s attention to the fact that, apart from the cloud parallax,

these parameters are not the same as the model parameters enter-

ing our data equation. Additionally, our toy-model is stochastic

due to the presence of the intrinsic scatter. Therefore, as demon-

strated in Appendix B, the mean values of the Stokes parameters

of a cloud and the values characterizing the covariance induced

by the intrinsic scatter must be read from the simulated data

before observational noise in parallax and polarization are intro-

duced. To do so, we segment the ‘clean’ data sets at the input

cloud distance(s) and we estimate the mean and covariance of

the polarization induced by the cloud to the polarization of stars

behind the cloud (but in front of the more distant cloud, if any).

We refer to these estimates as the “true” values in the remainder

of this paper.

We show in Fig. 1 two examples of simulated data for a

single-cloud case (left) and a two-cloud case (right) applied

to a sample of stars typical to intermediate to high Galactic

latitudes for a circular sky area with a diameter of

◦

(see

Appendix A.1). We show the relative Stokes parameters as a

function of distance modulus (

5 log(

)

−

, where

is the

star distance in parsec). The top row shows the simulated data

before noise in parallax and polarization is added and the bottom

row shows the same simulated data sets when noise is added.

The parameter values used to generate both mock data from

our toy model (see Appendix A.3) are reported in Table 1. The

two-cloud LOS is chosen such that (i) the far-away cloud alone

induces about half the polarization of the nearby cloud alone, (ii)

that the presence of the far-away cloud is only clear in one of the

two Stokes parameters due to the POS orientation of

permeat-

ing the second cloud, and (iii) that there are at least 20 stars in the

background of each cloud. The simulated data sets include uncer-

tainties on stellar parallaxes based on

Gaia

performance (see

Appendix A), realistic uncertainties on individual Stokes param-

eters as expected for WALOP-N (the northern instrument that

will be used for the P

ASIPHAE

survey) with 5 min of exposure

time (see Appendix A.4 and Fig. A.3), and include a prescription

for the intrinsic scatter. For completeness, we report in Table 1

the true values of the mean Stokes parameters and of the ele-

ments of the covariance matrix corresponding to the intrinsic

scatter corresponding to the two examples shown in Fig. 1. These

values correspond to the parameters entering our data equation

and, ultimately, should be retrieved from the application of an

inversion method to the data.

For the remainder of this section, we use these two exam-

ples of mock starlight polarization data to explore the sensitivity

of our log-likelihood with respect to the different sampled

parameters and to validate our implementation.

3.3. 1D conditional posterior distributions

As a first test of implementation, we perform 1D likelihood

scans through the parameter space. This exercise also allows us

to study the sensitivity of the likelihood, and therefore of the

observables, to each model parameter individually. We present

in Fig. 2 the conditional log-likelihood curves corresponding to

the scans using the one-cloud LOS mock data shown in Fig. 1

(bottom left). We show in Fig. 3 the conditional probability

distribution function (PDF) corresponding to these scans. They

are not estimates of the 1D marginalized posterior distributions

obtained from a fit since all parameters are not varied during the

scans but are fixed to their true values. An actual fit to this partic-

ular simulated data set is performed in Sect. 3.4. The validation

of our method and its implementation for several realistic cases

A164, page 6 of 30

V. Pelgrims et al.: LOS tomography of the dust polarization sky

single-cloud LOS

two-cloud LOS

(

V

V

)

Stellar Stokes [%]

(

V

V

)

Stellar Stokes [%]

Distance Modulus

(

V

V

)

Stellar Stokes [%]

Distance Modulus

(

V

V

)

Stellar Stokes [%]

Fig. 1.

Example of a single-cloud (left) and two-cloud (right) simulated data set. We show the stellar

(green circles) and

(blue diamonds)

Stokes parameters as a function of distance modulus (

). Top and bottom panels show the same data set. Top panels do not include observational

noise (they are the “true” data points) while bottom panels do include uncertainties in both parallax and polarization (shown with errorbars). The

vertical red lines indicate the input distance modulus of the clouds. The horizontal green and dashed-purple lines indicate the input (cumulative)

mean Stokes parameters (

and

, respectively) before the inclusion of intrinsic scatter and observational noise, i.e.,

in Eq. (1).

Table 1.

Setups for the simulated data set used in Sect. 3 and illustrated in Fig. 1.

Data set

Cloud parameters

“True” values

cloud #

max

turb

int

(pc)

(

)

(

◦

)

(

◦

)

−

(%) (%) (%)

(%)

Single-cloud

cloud 1

400

0.8

0.2

0.49

0.30

0.01

0.02

–0.01

Two-cloud

cloud 1

400

2.0

–15

0.2

1.56

–1.01

0.10

0.12

0.10

cloud 2

1300

1.0

0.2

0.00

0.91

0.17

0.18

0.16

Notes.

The labels of the parameters follow the notations given in the text.

will be presented in the Sect. 4 where the performance and the

limitations of the method will be assessed.

First, we note that the input parameter values always fall in

the interval where

log(

)

−

max(log(

))

−

meaning that the

input values fall inside the approximated 68% credible interval.

Second, we see that the shapes of the conditional log-likelihood

curve, and of its corresponding conditional PDF, from the explo-

ration of the cloud distance shown in Figs. 2 and 3 are somewhat

surprising while the curves obtained for the other parameters

look quite conventional. The very reasons for the unconventional

shapes of the cloud-distance curves come from the unevenly dis-

tributed constraints (the stars) on the parallax (distance) space,

their unequal uncertainties along that axis, the smearing in the

foreground and the background that the latter can generate, and,

last but not least, the unequal constraining power of each star in

the fit since polarization uncertainties are unequal. A star with

large polarization uncertainties will constrain the fit less, and

thus the position of the cloud along the LOS, as compared to

a star with small polarization uncertainties. We illustrate part of

this complexity in Fig. 4 in which we repeat the top left panel of

Fig. 3 but where the true and observed parallaxes are indicated by

vertical segments. It is clear that the likelihood of having a cloud

with any distance between two distant constraints is constant

and that the steepness of the variations depends on the parallax

uncertainties. In general, standard statistics are not appropriate

for characterizing the posterior distributions on cloud distances

A164, page 7 of 30