Press Release
Open Access
Issue
A&A
Volume 684, April 2024
Article Number A162
Number of page(s) 25
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/202349015
Published online 23 April 2024

© The Authors 2024

Licence Creative CommonsOpen 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.

1 Introduction

The polarization of starlight is a powerful probe of the magnetized interstellar medium (ISM). Starlight acquires a polarization due to dichroic absorption by aspherical interstellar dust grains, which align their minor axis with the magnetic field (e.g., Davis & Greenstein 1951; Andersson et al. 2015). The polarization position angle of starlight is parallel to the plane-of-sky (POS) component of the magnetic field, and the maximum degree of polarization is proportional to the column density of the polarizing dust through which the light beam passes. Since its discovery (Hiltner 1949; Hall 1949), the polarization of starlight has contributed significantly to the study of the magnetic field in our Galaxy and to our understanding of its role as an agent of the Galactic ecosystem, from the smallest to the largest scales (e.g., Spoelstra 1972; Ellis & Axon 1978; Goodman et al. 1990; Heiles 1996; Heyer et al. 2008; Nishiyama et al. 2010; Li et al. 2013; Berdyugin et al. 2014; Doi et al. 2023). Once stellar distances are known, starlight polarization could provide information on the properties of the magnetized ISM directly in three dimension (3D; Panopoulou et al. 2019), and with good spatial resolution given the high stellar density.

In recent years, the Gaia satellite (Gaia Collaboration 2016) has provided the data necessary for the precise localization in 3D space of more than a billion stars in our Galaxy (e.g., Bailer-Jones et al. 2021; Gaia Collaboration 2021, 2023; Lindegren et al. 2021). By combining stellar parallax and reddening data, several teams have been successful in reconstructing 3D tomography maps of the dust density distribution in large volumes centered on the Sun, up to around 3 kpc in the Galactic disk and up to around 1.2 kpc in the halo (Green et al. 2019; Lallement et al. 2019, 2022; Leike & Enßlin 2019; Leike et al. 2020; Vergely et al. 2022; Edenhofer et al. 2024), or in more focused areas. This formidable community effort has revolutionized our view of the 3D structure of dust distribution in the ISM, and it has already enabled the modeling and better understanding of some of the main structures in our cosmic neighborhood and their history (e.g., Pelgrims et al. 2020; Alves et al. 2020; Das et al. 2020; Bialy et al. 2021; Zucker et al. 2022; Großschedl et al. 2018; Marchal & Martin 2023; Ivanova et al. 2021; Tahani et al. 2022). Such a 3D mapping, with the additional knowledge of magnetic field properties, would be certain to enable breakthroughs and discoveries in several research topics, as discussed below and in Pelgrims et al. (2023).

For example, a tomographic view of the magnetized ISM would offer new avenues to address open questions such as the role of the magnetic field in star formation (see, e.g., Mouschovias et al. 2006; McKee & Ostriker 2007), and the search for the sources of ultra-high energy cosmic rays (see, e.g., Boulanger et al. 2018; Magkos & Pavlidou 2019; Tsouros et al. 2024). Likewise, by providing measurements of the polarization properties for each individual dust cloud in 3D space, such a tomography map would enable significant progress in the modeling and characterization of the dusty magnetized ISM as a contaminating foreground to observations of the polarization of the cosmic microwave background (Tassis & Pavlidou 2015; Martínez-Solaeche et al. 2018; Pelgrims et al. 2021). Hence, it would help clear the path for an unbiased study of one of the first moments of the Universe’s history. This holds true, although less directly, for the characterization of the Galactic synchrotron emission as demonstrated by Panopoulou et al. (2021). Combined with multiwavelength observations of the polarized emission of dust in the submillimeter, knowledge of the polarization properties of each individual dust cloud would also enable us to advance our modeling of astrophysical dust (e.g., Hensley & Draine 2021), to better determine its composition and understand its interaction with its cosmic environment, and to better assess its role as a building block of life.

Mapping the dust distribution in 3D over large volumes (kiloparsecs) required millions of stellar extinction and distance measurements. A similar amount of data is likely required for magnetic field mapping in such a volume. Currently, the availability of stellar polarization data is limited (Panopoulou et al. 2023), with only the inner Galaxy having millions of stars with measured near-IR polarization (GPIPS survey, Clemens et al. 2020). In the future, planned optical polarimetric surveys will deliver millions of measurements throughout the entire sky (Magalhães et al. 2005, 2012; Tassis et al. 2018; Covino et al. 2020). It is therefore timely to develop techniques to analyze these forthcoming datasets in such a way that would make it possible to extract from the data most of the information on the properties of the dusty magnetized ISM in 3D.

The observed polarization of each single star is the integrated effect of the dichroic absorption from all ISM clouds lying between us and the star. This line-of-sight (LOS) integration needs to be inverted in order to derive the complex 3D structure of the magnetized ISM from starlight polarization and distance measurements. Recently, we have developed the first standalone method to perform this LOS inversion through Bayesian modeling (Pelgrims et al. 2023). This method, which takes into account all sources of uncertainties in polarization and parallax measurements, works on a per LOS basis. It has the advantage that it can be easily automated to run on a large set of sightlines. We have extensively tested the method and its performance based on simulated data. We further applied it to existing data for two sightlines and demonstrated that this new method leads to results that are fully consistent with previously obtained results but within a robust Bayesian framework and, therefore, put on a more solid footing.

In this paper, we continue our preparation for the large datasets to come. We aim to develop a pipeline capable of inverting measurements of parallax and starlight polarization for an actual, extended volume of 3D space in order to derive the properties of the dusty magnetized ISM in 3D. We develop a first pipeline based on our LOS-inversion method and apply it to observations taken for an extended region of the sky. Our goal is to obtain, with as few assumptions as possible, a first tomography map of the POS component of the magnetic field in dusty regions from which 3D properties of the magnetized ISM can be accessed.

For the purposes of this work, we carried out a survey of starlight polarization for a continuous region of the sky covering about four square degrees. We present our survey and the resulting dataset in Sect. 2. We devote Sect. 3 to the description of our data analysis pipeline and its application to our data to obtain the first degree-scale tomography map of the magnetized ISM from starlight polarization and parallax measurements. Section 4 presents our main results and how we produce and visualize our 3D map of the POS component of the magnetic field in dusty region from the posterior distributions output by our Bayesian analysis. Our 3D map extends up to 3 kpc distance and covers a sky region of about four square degrees. We discuss our results in comparison with other probes of the (magnetized) ISM in Sect. 5 and conclude in Sect. 6.

2 Dataset

To demonstrate the feasibility of starlight-polarization-based tomography of the diffuse magnetized ISM, Panopoulou et al. (2019) obtained dense optical polarization measurements of stars for two circular beams of 9.6 arcmin radius in a sky region that they identified as being likely to exhibit complexity along the distance axis based on inspection of H I velocity data and polarized dust emission maps. These data indeed suggest possible variations of the polarization signal both in the POS and along distance as several components can be identified. As they showed through the analysis of their polarization data along distance, and as we recently confirmed (Pelgrims et al. 2023), one of the two beams likely intersects at least two dust clouds while the other one likely intersects only one. Showing complexity both along the LOS and in the POS, this region therefore seems to be well suited to develop and test tomography methods to reconstruct the magnetized ISM in 3D. Hence, we carried out a survey to expand that region of the sky with optical starlight polarization data. The final surveyed region, comprising about four square degrees, and its location in the sky is shown in Fig. 1.

thumbnail Fig. 1

Sky location of the surveyed area of about four square degrees. Left: full-sky map of the dust emission as seen by Planck at 353 GHz (Planck Collaboration XII 2020). The color represents the intensity of dust emission on logarithmic scale. The line-integral-convolution texture shows the polarization angle of the dust emission rotated by 90°. Middle: a zoom-in of the map toward the surveyed regions, which includes part of the North Celestial Pole Loop on the East of the map. Right: a closer view of the surveyed region. Black segments indicate the polarization orientation from the stars in our survey and from Panopoulou et al. (2019). The segments are scaled according to the polarization fraction. Unpolarized stars appear as dots. The blue, horizontal segment in the bottom right corner shows the scale for a 1% polarized star. Outlier candidates (see Sect. 3.3) are not shown.

2.1 Survey strategy and observation plan

We aimed at obtaining polarization measurements for a large, continuous region centered on Galactic coordinates (l, b) ≈ (103.5°, 22.25°) with complete star samples limited in magnitude. Due to limited observing time, it was infeasible to perform a deep, photometrically uniform survey for the entire area. Hence, we relied on H I observations to gauge a priori which part of the targeted region likely intersects dust clouds at large distances, in order to increase the number of stellar polarization data accordingly.

As discussed by Panopoulou et al. (2019), and also shown in Fig. 2, the averaged H I-velocity spectrum (measured in the Local Standard of Rest, LSR) in this region of the sky shows the existence of two very distinct components. The dominant component is centered on υLSR ≈ −2.5 km s−1 and we refer to it as the low velocity cloud (LVC). The second component, with a lower amplitude, is centered on υLSR ≈ −50 km s−1 and we refer to it as the intermediate velocity cloud (IVC). This velocity component may be a southern extension of the IV Arch as identified by Kuntz & Danly (1996) who determined that it is located at large distance (≳ 1 kpc). We used the HI4PI survey (HI4PI Collaboration 2016) to look at the spatial distribution of the intensity maps obtained from the integration of the velocity spectra in channels corresponding to the peaks observed in the averaged velocity spectrum.

These data are shown in Fig. 2 where we also show the outline of the region for which we obtain a tomography reconstruction in this work. As can be seen clearly in the right-hand bottom panel of Fig. 2, the H I velocity spectra show power in the IVC range in the eastern part of the region (with l ≳ 103.5°). Thus, we concluded that this region is more likely to contain dust clouds at larger distances than the other (western) part of the region. Hence, as the stellar magnitude generally increases with stellar distance, we conduct a deeper survey in that part of the sky (with a limiting magnitude R ≲ 15.5 mag, compared to R ≲ 14 mag in the remaining area); to increase the number density of data points, in particular at large distances. Despite this choice that resulted in more sparse stellar polarization measurements on the right part of the region, a hint of distant clouds is still found there (see Sects. 4 and 5). This indicates that even with a shallower polarization survey distant clouds could still be recovered - though see Sect. 5 where we discuss the reliability of cloud detection.

We have conducted our survey so that the entire region is photometrically complete in the R band up to 14 mag (shallow survey), and up to 15.5 mag for the part of the region located at l ≳ 103.5° (deep survey). The upper panel of Fig. 3 shows the locations of stars on the sky for both the shallow and the deep survey. We used the R-band magnitude and sky coordinates from the USNO-B1.0 catalog (Monet et al. 2003) to plan our observations.

2.2 Observations and data reduction

Observations of optical polarization were conducted using the 1.3-m telescope situated at Skinakas Observatory1 in Crete (1750 meters above the sea level, 24°53′57″E, 35°12′43″N). The telescope is equipped with the RoboPol polarimeter, which consists of two adjacent half-wave retarders, with their fast-axes rotated by 67.5°, relative to each other, followed by two Wollaston prisms with orthogonal fast-axes (Ramaprakash et al. 2019). This setup splits each incident ray into four rays with different polarization states on a single CCD, which provide information on the q = Q/I and u = U/I normalized Stokes parameters in the instrument’s reference frame (see Eq. (1) in King et al. 2014) with only one exposure. In this way, the use of a fixed instrument configuration eliminates random and systematic errors resulting from changes in the sky, imperfect alignment, and non-uniformity of rotating optical elements, as the instrument has no moving parts aside from the filter wheel. To enhance the signal-to-noise ratio (S/N), a special mask was placed in the center of the telescope focal plane where systematic uncertainties have been estimated to be lower than 0.1% in the degree of polarization (Skalidis et al. 2018; Ramaprakash et al. 2019).

The observations were carried out star-by-star over three observation seasons, from May 2019 to November 2022. Each star was measured in the mask. Optical polarization measurements were obtained for each star in the Johnsons-Cousins R band. For each star, the observation exposure time was estimated on-the-fly so as to guarantee that photon-noise-driven uncertainties fall below the estimated uncertainties from instrumental calibration (which turned out not always to be the case as discussed below). Zero- and highly polarized polarimetric standards were observed every observing night to monitor the instrumental polarization and polarization angle zero-point through time, and to estimate the corresponding uncertainties as described in Blinov et al. (2021, 2023). We obtained instrumental uncertainties in both q and u at the level of 0.1%.

Data reduction was performed with the standard RoboPol pipeline (King et al. 2014; Panopoulou et al. 2015; Blinov et al. 2021). For any given source, we produced the stacked image of all observations and deduced the linear Stokes q and u parameters and their photon-noise-driven uncertainties through differential aperture photometry as outlined in Ramaprakash et al. (2019). Then we removed the contribution from instrumental polarization and added in quadrature observational and instrumental uncertainties as: (1) (2)

and similarly for Stokes u, where the superscripts “measured” and “instr” refer to the values obtained from differential photometry and from estimation of the instrumental polarization, respectively. The polarization data are given in the equatorial coordinate system and follow the IAU convention for the polarization position angle (zero at north, increasing toward east). In total, we obtained reliable optical polarization data for 1530 stars, spending approximately 153 telescope hours.

thumbnail Fig. 2

H I velocity data from the HI4PI survey in the sky area toward the region surveyed in starlight polarization. The top panels show the brightness temperature as a function of the gas velocity measured in the LSR. The blue spectrum corresponds to the velocity spectrum averaged over the full region. Two dominant peaks with velocities in the LVC (≈ − 1 km s−1) and IVC (≈ − 50 km s−1) ranges are seen. The bottom panels show the column density maps resulting from the integration of the velocity spectra in the ranges depicted by shaded regions marked in the top panels in the LVC (left) and IVC (right) ranges. The maps reveal a high degree of complexity in morphological structures. The magenta outline indicates the sky region of 3.8 sq. deg. for which starlight-polarization-based tomography is obtained in this work (see Sect. 3.4).

2.3 Cross-match with Gaia and quality cuts

For the purpose of this work, we complemented our new polarization measurements with data from Panopoulou et al. (2019) for 192 stars. As we need estimates of stellar parallaxes and their corresponding uncertainties to perform the tomography decomposition along distance, we cross-matched our sample with polarization measurements with the Gaia DR3 catalog (Gaia Collaboration 2023).

We used a cross-match radius of 5 arcsec around the USNO-B1.0 coordinates of the targets. We chose to use the USNO-B1.0 coordinates as the astrometric accuracy of the RoboPol pipeline output varies with position on the field of view and can be limited to a few arcseconds in some cases as a result of distortions and the 4-spot pattern of the stars which are modeled within the software2. For most sources, there was a unique Gaia source found within the search radius, to which we assigned the match. Some sources had multiple matches with the Gaia catalog, as a result of the proximity of the target star with another star. The cross-match was not always straightforward mainly due to the lack of precise astrometry in the USNO-B1.0 catalog. This problem mostly happened when, in the USNO-B survey, adjacent stars (approximately less than 2″ spatial separation) were blended3. We visually inspected all cases where multiple matches were found by comparing the raw RoboPol images with the USNO-B1.0 and Gaia catalogs. We also checked for possible misidentifications by comparing the R- and G-magnitudes of the sources after the cross-match, correcting for a couple cases where a fainter Gaia star had been mistakenly associated with a bright USNO-B1.0 star. Finally, we note that we observed several faint stars, which exceed the photometric magnitude limits ofour survey (R > 15.5 mag), because they were close (around 1″ distance) to the target stars. These stars happened to lie within the mask of RoboPol during the observations. The resulting obtained S/N in their degree of polarization is usually low. In some cases, this proximity led to photometrically blended measurements, which we disregarded in our analysis.

The final cross-match was successful for 1698 stars for which we retrieved their Gaia identifier, their G-band photometric magnitude, their parallax and corresponding uncertainty, and their renormalized unit weight error (RUWE). The latter measures the quality of the single-star model to account for astromet-ric Gaia observations, and must be used as a quality criterion to guarantee the reliability of the solution, and thus of the parallax estimates. We used stars with RUWE < 1.4 (as recommended) to avoid unreliable measurements that could occur because of blended sources for example (Lindegren et al. 2018, 2021). Among the successfully cross-matched stars, 24 stars do not have parallax information and 226 stars have RUWE ≥ 1.4 or are unknown. Consequently, after applying the quality criterion we have 1448 stars with both reliable optical polarization measurements and reliable parallax estimates. This is the star sample that we use for our tomography of the magnetized ISM.

Some properties of our final dataset of 1448 stars are given in Figs. 35, which illustrate the non-homogeneous character of our sample.

The different cuts in R-band magnitude used to design our survey based on H I complexity are clearly seen in the sky distribution shown in the top panel of Fig. 3. The histogram in the bottom panel of the same figure clearly shows that fainter stars (R ≥ 14 mag) generally have larger distances (obtained simply by taking the inverse of the parallax) than brighter ones.

Consequently, the density of stars, in particular at large distances, is much larger in the eastern half of the observed region than in the western one. We also discuss this in Sect. 3.4. Figure 4 shows the distribution of parallax and parallax S/N in our sample. Most of our targets (with RUWE < 1.4) have parallax S/N higher than 10, demonstrating the reliability of our distance markers through the ISM. The distributions of the Stokes parameters and the (total) uncertainties of the entire sample with reliable estimate of the parallax are shown in Fig. 5. It is seen from the 2D distribution of the Stokes parameters that a certain number of stars have low polarization while the majority have polarization degree at the per cent level or higher4. The uncertainty on the Stokes parameters is at the level of 0.19% in both q and u for a large fraction of our sample, and therefore dominated by systematic uncertainty from the instrument calibration, while the photon noise contributes significantly for a subset of measurements. We decided to keep all measurements since we are confident in their uncertainty estimates.

thumbnail Fig. 3

Sky distribution (top) and distance (modulus) distribution obtained by inverting the parallax (bottom) of our stellar sample with reliable parallax estimates. The sample is divided according the USNO-B R-band magnitude. Blue, red and black correspond to stars with R < 14 (shallow survey), R ≥ 14 (deep survey) and unknown values due to identification mismatch. The histograms on the bottom panel are stacked.

thumbnail Fig. 4

2D distribution of the number of stars in our catalog in the parallax – S/N plane. For better visualization, a dozen stars with large parallax values and or large parallax S/N are not included in this plot.

thumbnail Fig. 5

Distribution of the polarization properties for the stellar sample with a reliable parallax estimate. Distributions of the Stokes parameters (q’s and u’s) and of the logarithm of their uncertainties (log10(σq) and log10(σu)) are shown in the top and bottom panels, respectively. The horizontal and vertical dashed lines are used for visual reference. They indicate the q = 0 and u = 0 loci in the top panel and indicate the values of 0.1% and 1% in polarization uncertainties in the bottom panel.

2.4 Coordinate system conversion

The Stokes parameters are measured in the equatorial celestial coordinate system. We construct our 3D map of the POS component of the magnetic field in the Galactic coordinate system, as this seems natural in the context of Galactic tomography. Although we could perform the tomographic decomposition in the equatorial coordinate system and then convert the resulting 3D map in the Galactic coordinate system, we prefer to first convert the starlight polarization data and then perform our analysis to obtain the tomography result directly in the Galactic coordinate system. Because of the change of coordinate system, the values of the polarization position angles, and therefore of the Stokes parameters, change to account for the change of the orientation of the meridian from one reference frame to the other at the location of the stars. The change of coordinate system thus involves a rotation of the polarization plane which is implemented using the rotation matrix (e.g., Tegmark & de Oliveira-Costa 2001) (3)

where the rotation angle (ψR) is defined locally. It depends on the celestial coordinates (right ascension and declination) of a star (α, δ) and on the celestial coordinates of the North Galactic pole (αNGP, δNGP), and is given by (e.g., Hutsemékers 1998) (4)

where we use the two arguments arctangent function to place the resulting angle in the correct trigonometric quadrant. The Stokes parameters of a star in the Galactic coordinate system are thus obtained from the Stokes parameters in equatorial coordinate system through the following: (5)

Similarly, it can be shown (e.g., see Appendix A of Planck Collaboration Int. XIX 2015) that the noise covariance matrix of the Stokes parameters in the Galactic reference frame is obtained from the rotation matrix and the noise covariance matrix in the equatorial coordinate system through (6)

where RT is the transpose of the rotation matrix. Introducing a = cos(2ψR) and b = sin(2ψR), the elements of the noise covariance matrix are thus obtained as follows: (7)

The noise covariance matrix of the Stokes parameters remains unchanged from one coordinate system to the other only when there is no noise covariance between q’s and u’s and when the noise uncertainties in both q’s and u’s are equal. In any other situation, the noise covariance matrix changes and, most notably, a non-zero off-diagonal term arises simply due to the change of coordinates. The use of the exact analytical formalism given above allows us to avoid issues related to the polarization bias in low S/N regime and to avoid the need for estimating the noise covariance matrix through Monte Carlo treatment.

Table 1

Polarization catalog (abbreviated).

2.5 Stellar polarization catalog

The polarization catalog is made publicly available through CDS. The columns contained in the catalog are described in Table 1. The coordinates given are from Gaia DR3 at epoch 2016. We include all stellar polarization measurements, including those that were identified as outliers and those that did not return a reliable distance (due to no match with Gaia or high RUWE values). We distinguish between sources used in the tomographic reconstruction in the catalog with a usage flag. The usage flag is assigned a value of 0 for stars that did not have a reliable distance estimate, 1 for stars that are used in the reconstruction, and 2 for stars identified as outliers (see Sect. 3.3).

3 Method and data analysis

In this paper we aim to obtain a 3D tomographic decomposition of the dusty magnetized ISM for an extended region of the sky, from stellar measurements of optical polarization and distance. As there is no currently available method to perform a 3D inversion of an actual volume of stellar polarization data, we design an analysis pipeline based on the LOS-inversion method that we presented in Pelgrims et al. (2023) and which is implemented in the BISP-1 (Bayesian Inference of Starlight Polarization in 1D) Python code5.

3.1 BISP-1

Our maximum-likelihood method, implemented in BISP-1, makes it possible to reconstruct the dusty magnetized ISM along a single LOS using starlight polarization and distance (parallax) only. It assumes that the effect of multiple clouds along the LOS is well approximated by the vector addition of the linear Stokes parameters induced by each cloud separately – which holds for the typically low polarization levels in the diffuse ISM (Martin 1974; Patat et al. 2010; Panopoulou et al. 2019). Relying on the nested sampling method (Skilling 2004) implemented in the code dynesty (Speagle 2020), our algorithm is able to determine the number of components (dust clouds) along the LOS and to determine the distance and polarization properties of each component using six parameters per component: cloud distance (dC), cloud mean polarization (qC, uC), and three parameters to characterize the covariance matrix Cint encoding the intrinsic scatter of stellar polarizations arising as a result of ISM turbulence. As such, the method makes it possible to recover the stellar-polarization source field which directly informs on the local orientation of the POS component of the magnetic field (the position angle) and on the local degree of polarization. The latter is related to the dust grain density, the dust grain polarization efficiency and on the angle made by the magnetic field lines with the POS, as explained in (e.g., Pelgrims et al. 2023). In this picture, it is implicitly assumed that the aspherical dust grains always align their shortest axis with the local orientation of the magnetic field, at least statistically. This assumption is expected to be true for the diffuse ISM and agrees with current state of dust alignment theory (e.g., Draine & Hensley 2021; Hensley & Draine 2021).

The model that we use to reconstruct the dusty magnetized ISM assumes that the dust clouds can be represented as thin layers distributed along the sightline. That is, we consider that the typical extent of dust clouds along the LOS is smaller than the typical separation of stars or, in practice, smaller than the distance range spanned by the number of stars needed to allow for the detection of a cloud given the amplitude of the polarization it induces and all sources of scatter in the polarization data. We expect this assumption to hold at high and intermediate Galactic latitudes and for clouds of the cold neutral medium and molecular clouds (Heiles 1976; Zucker et al. 2021; Marchal & Martin 2023). It is worth noting that the thin-layer approximation also seems to hold in the denser regions of the ISM, at least for some sightlines (Doi et al. 2023).

Our method (BISP-1) further assumes that all stellar measurements trace the dusty magnetized ISM only (as opposed to having intrinsic polarization), and it does not implement spatial variation in the POS. To circumvent these limitations while benefiting from the strengths of our LOS-inversion method, we embed it in a multistep process to form our 3D inversion pipeline, as described below.

3.2 Design of the 3D inversion pipeline

The broad outline for the workflow of our 3D inversion pipeline is illustrated in the flowchart given in Fig. 6. From the original sample of stars in the extended region, we first identify stars that likely trace the magnetized ISM (as opposed to the intrinsically polarized candidates). Then, for a set of sightlines that span the observed region, we create two sets of overlapping sub-samples centered on each LOS, as detailed in Sect. 3.4. The first is meant to capture ISM properties at small distances despite the low stellar density; the second corresponds to the highest angular resolution we can achieve. We then perform the LOS inversion for each LOS; first on the subsamples which connect different sightlines at small distances, and then at higher resolution taking into account the information obtained in the previous step. Finally, we produce the 3D map of the magnetized ISM from the posterior distributions obtained for each LOS individually. Each step is detailed in this section and in Sect. 4, and illustrated by applying the method to our dataset.

thumbnail Fig. 6

Flowchart for the 3D inversion pipeline, from the stellar catalog with reliable parallax estimates (top) to the 3D maps of differential quantities (bottom).

3.3 Identification of outliers

Some stars may exhibit intrinsic polarization, possibly due to the existence of a circumstellar disk or other asymmetries in the object (e.g., Fadeyev 2007; Cotton et al. 2016; Gontcharov & Mosenkov 2019). These stars usually show either a higher degree of polarization or an unrelated polarization position angle as compared to their neighboring stars, and sometimes they display both. As BISP-1 assumes that all starlight-polarization data points trace the dusty magnetized ISM, the first step is to remove from the original sample all stars for which their polarization is unlikely to be of only interstellar origin. This includes intrinsically polarized stars and also any other outliers.

To identify these stars, we adopt recursive sigma-clipping in groups of neighboring stars. Based on the sky coordinates and distance estimate of every star, we identify the group of N neighbors of each star in 3D space using their heliocentric Cartesian coordinates, ignoring the distance uncertainties from parallax measurements. We estimate the weighted-average Stokes parameters and associated covariance in the (q, u) plane from the neighbors, and then compute the Mahalanobis distance of the polarization of the central star (s) as compared to the 2D bivariate distribution from the neighbors as (8)

where the noise covariance matrix of the polarization of the central star is added to the covariance matrix from the bivariate distribution to compute the Mahalanobis distance. If the Mahalanobis distance exceeds some threshold, the probability that the polarization of the central star is drawn from the same parent distribution as the neighbors (assumed to be representative of the magnetized ISM) is small and, therefore, the central star must be identified as an outlier. Repeating this process for all the stars in the original sample and running the whole process until no additional outliers are identified yields a catalog of outliers and a “clean” sample of stars whose polarization likely traces the dusty magnetized ISM. Only the clean sample of stars is then considered to infer the dusty magnetized ISM. The exact list of outliers depends on our specific choice for the size of the neighbor groups and the adopted threshold in significance level. In addition, the sensitivity of the “clean sample” membership to these parameters should also ideally be tested against distance uncertainties. However, if we can recover the 3D dusty magnetized ISM from the clean sample of stars, we can test a posteriori the hypothesis that the polarization of a given star is given by the magnetized ISM only, as we do in Sect. 5.5. This has the potential to lead to a somewhat more robust list of candidate targets for being intrinsically polarized stars, or at least outliers, and can trigger follow-up observations. However, we notice that caution has to be made for the choice of the significance-level threshold. A too strong selection criterion would discard stars that merely pick up fluctuations of the magnetized ISM. This would subsequently lead to an underestimation of the turbulence-induced intrinsic scatter which we want to avoid. This point is further discussed in Sect. 5.5.

We apply the recursive sigma-clipping approach to the full sample of stars discussed in Sect. 2.3. We adopt a size of N = 30 to build the groups of neighbors and we choose to flag every star as an outlier if its polarization shows a probability of less than 1% for it to be drawn from the same parent distribution as the neighbors. In our case, the number of outliers remained constant after three iterations. Using these parameters, 18 stars out of 1448 are identified as outlier candidates. This represents a fraction of 1.2% of our full sample. In Fig. 7, we show the stellar polarization data plotted against the distance modulus (µ = 5 log10(d) − 5) and highlight the outliers. These outliers are not considered in the analysis discussed below.

3.4 Definition of subsamples

BISP-1 assumes that all the stars in a sample lie along a narrow, one dimensional beam. It consequently returns the structure of the dusty magnetized ISM averaged in the POS over the sky region spanned by the input sample, which we refer to our “beam”. The method also relies on the dust-layer model that we have introduced (Pelgrims et al. 2023) and which we expect to hold true as long as the magnetic field and dust density do not vary appreciably in our beam. The validity of these assumptions is tested by the data in the following. Depending on the geometry of the volume filled by the star sample, the averaging scale in the POS may depend on distance. For example, if a conical geometry is chosen to define the star subsamples, the averaging scale at small distances is much smaller than the one at large distances. A cylindrical geometry would keep the averaging scale constant. The number density of constraints (i.e., of stars) as a function of distance also depends on the chosen geometry of the beam, the specific size of the volume encompassed by the data, and the actual 3D spatial distribution of stars. Hence, there is a trade-off between having a sufficient number of data points to constrain the model and the achieved resolution. The resolution also needs to be sufficiently good to minimize the POS variations of ISM properties in the beam and thus ensure that the intrinsic scatter, which we fit for, is not dominated by POS variations of the mean ISM properties.

We find that running the BISP-1 decomposition solely on subsamples defined according to conical beams is not appropriate. The reason is that the spatial distribution of nearby stars is sparse and that, unless a very large opening angle is chosen, nearby clouds may be missed or their parameters loosely determined due to the absence of a sufficient number of constraints in the required range of distances. If a large opening angle is chosen in order to be able to capture those nearby clouds, then the signal of any faraway structure would be averaged out and therefore likely missed. Alternatively, a cylindrical geometry for our beam would lead to a constant resolution in the transverse direction to the distance axis. However, at high and intermediate Galactic latitudes the number density of stars decreases substantially after approximately 1 kpc. Faraway clouds would then be missed, due to the sparsity of data points, unless a large cylinder radius is chosen. A cylindrical geometry would also imply that any detail at small distances would be averaged out over large angular scales. Thus, a cylindrical geometry alone is also not appropriate to define our beam samples.

To guarantee a good angular resolution at all distances, and to avoid missing clouds or only placing loose constraints on their distance and polarization properties, we adopt a two-step hierarchical decomposition process. The first step is performed on samples defined according to a hybrid geometry for our beam: the beam follows a cylindrical geometry at low distances and a conical geometry at large distances. For the second step, the star samples are defined in a beam with conical geometry only. The idea is to perform the BISP-1 decomposition in the second step using priors defined from the posterior distributions obtained in the first step. In this way, constraints on distances and polarization properties of nearby clouds are obtained from the decomposition on hybrid-beam samples (with larger number of stars at lower distances) while good angular resolution is achieved, even at low distances, from the conical-beam samples. We use the same opening angle for the conical part of the hybrid beam and the purely conical beam so that they match at large distances.

In Fig. 8, we show the distance-dependent angular radius of our beams along with the corresponding physical scale in the POS for the specific choice of opening angle and cylinder radius that we use. Our hybrid beam centered on a given LOS has an angular size that depends on distance. Any star is considered as part of a given subsample if its angular separation to the LOS is lower than the maximum between the opening angle of the cone and the distance-dependent angular size of the cylinder evaluated at the star distance. For the conical beam, the angular separation cutoff is constant.

To cover an extended region of the sky, such as the one we observed, we adopt a moving-window strategy. That is, we sample the observed region with a large number of sightlines. For each LOS, we define a subsample of stars according to our choice for the beam geometries. The star samples of neighboring sight-lines are overlapping and, in the first step, the overlap extends to larger angular scale for nearby stars than that for distant stars. This strategy ensures a continuous scan of the observed region and also implies that the results of the LOS decomposition of neighboring sightlines will not be independent.

For our stellar polarization measurements spanning a region of about four square degrees, we sample the sky according to a HEALPix tessellation (Górski et al. 2005; Zonca et al. 2019) with the resolution parameter Nside = 512. Each pixel center defines a LOS which we take as the symmetry axis of our beam. The angular separation of neighboring sightlines is about 6.9 arcmin. To define our hybrid beam, we adopt a value of 2.5 pc for the cylinder radius and an angular radius of 13.74 arcmin for the cone. The conical beam has the same value of 13.74 arcmin for the angular radius. Each conical beam spans a sky area of about 0.16 square degrees. We choose the angular radius of the conical beam so that (i) at least 20 stars are contained within each LOS and (ii) every pixel has at least two neighbors to ensure continuity at the edges of the map. The sky region is covered by an ensemble of 287 sightlines. Panels a–c of Fig. 9 show the number of stars per bins of distance modulus for all the 287 sightlines covering the regions of the hybrid beams, of the conical beams, and of the cylindrical beams, respectively. We see that the density of stars is generally very low at distances smaller than 300 pc for the conical beams and at distances larger than 1 kpc for the cylindrical beams. This justifies the use of the hybrid beams in extracting information of nearby and distant ISM clouds. With our choice of beam parameters, the geometry of the hybrid beam transitions from cylindrical to conical at a distance of about 630 pc (µ ≈ 9). Figure 10 shows the number of stars per beam samples for the conical beam geometry. The higher density for points at l ≳ 103.5° is due to the aforementioned survey strategy (Sect. 2.1).

thumbnail Fig. 7

Stellar relative Stokes parameters (q, u) versus the distance modulus for the full sample. The qV and uV (in the Equatorial coordinate system) are shown by green circles and blue diamonds, respectively. The vertical error bars indicate the uncertainties (noise and systematic) on the Stokes parameters and the horizontal error bars represent the uncertainties on the star distance modulus converted from the 1σ parallax uncertainties. Outliers identified from the iterative sigma-clipping approach in groups of nearest neighbors are highlighted with red-filled circles for q and purple-filled squares for u. To facilitate the visual clarity of the plot we restrict the range of (q, u) and µ values. Ten stars have µ < 5 and 30 stars have µ > 14.

thumbnail Fig. 8

Beam sizes as a function of distance from the observer. Top: angular aperture radius of our hybrid beam (continuous gray) resulting from the combination of a conical beam (red dashed) and cylindrical beam (dot-dashed). Bottom: extent of the beams in the transverse direction to the LOS. Same line convention as in the top panel. In this example the cylinder radius is fixed at 2.5 pc and the angular radius of the cone at 13.74 arcmin. For the hybrid beams, the cylindrical geometry prevails at distances smaller than ≈630 pc, while the conical geometry prevails at larger distances.

3.5 LOS decomposition of the dusty magnetized ISM

Once the star subsamples are defined for the entire observed region, we independently apply the BISP-1 code to each sub-sample, first on the hybrid-beam samples and then on the conical-beam samples.

3.5.1 Step 1 – LOS decomposition in hybrid beams

For each LOS, we test the layer-model for one to four clouds along the distance with uniform priors on all the model parameters. The limits of the priors on cloud parallaxes are set so that the minimum allowed distance is 20 pc and the maximum allowed distance of the farthest cloud is the minimum between 3.5 kpc and the maximum distance of the stars in the analyzed sample. The upper distance limit may thus vary from one LOS to another. The value of 3.5 kpc (corresponding to a distance modulus of µ ≈ 12.7) originates from the fact that beyond this distance our data sample generally become very scarce as also seen in Fig. 3. In addition, we make sure that there are at least five stars between clouds. This is required by the BISP-1 code (Pelgrims et al. 2023).

The uniform priors for the cloud polarization parameters are set as follows. For the mean polarization parameters (qC and uC), we compute the maximum of the absolute values of both the stellar Stokes parameters in the star samples. This value is used to define the limits of the top-hat priors on both qC and uC. This definition of the prior limits, while not fully general, is valid in our case as we know that the extinction is dominated by nearby components as shown for example by comparing Planck dust column-density map and 3D star extinction maps (e.g., O’Callaghan et al. 2023). The diagonal elements of the intrinsic-scatter covariance matrix are positive definite and we require Cint,qq, Cint,uu ∈ [0, 10−4]. This is a very loose range for possible values as it allows for a spread in stellar Stokes parameter due to turbulence as high as 1% in degree of polarization. The off-diagonal element Cint,qu is initially bounded by ±10−4 but is further constrained to verify |Cint,qu| < (Cint,qq Cint,uu)1/2 inside BISP -1 to ensure that the covariance matrix is invertible.

We use BISP -1 to run the nested sampling experiment using 1000 live points and sample the parameter space until an uncertainty of around 0.01 is achieved on the log of the model evidence. Typically, this requires 15 000–80 000 nested sampling iterations, corresponding to several hundred million calls of the log-likelihood function of each model.

For each of the tested models, BISP-1 leads to estimates of the log of the evidence, the maximum log-likelihood value and to the estimated posterior distributions of all model parameters. As already demonstrated and pointed out in Pelgrims et al. (2023), the posterior distributions of the cloud parallaxes have generally complex shapes, mainly due to the sparse and uneven distribution of stars along distance, and might also be piled up on the far edge of the prior domains. The latter case happens whenever the data are not enough to determine the cloud properties or if there is no cloud to be found.

To deal with this peculiarity, we follow the same idea as in Pelgrims et al. (2023) and rely on the following automated analysis of the marginalized posterior distribution on the cloud parallax. The idea is that the value at maximum-likelihood must belong to (one of) the main mode(s) of the marginalized posterior distribution on the cloud parallax and that this mode must not be piled-up on the lower limit of the prior. If both criteria are verified, then we qualify the fit as valid. We then select from the 6×NC dimensional posterior distribution (where Nc is the number of layer in the model) all samples that belong to this mode for the remainder of the analysis. In practice, we analyze the marginalized posterior distribution using the peak-finder algorithm f ind_peaks of the SciPy Python library (Virtanen et al. 2020) which identifies all local maxima through simple comparison of neighboring values. We thus find local maxima of the marginalized posterior distribution and the range of parallax values corresponding to the extent of the corresponding peaks. We compute the fraction of the posterior distribution which corresponds to each peak and we consider that it is (one of) the dominant peak(s) if this fraction exceeds the threshold of 30%. Our results do not depend strongly on this choice.

We consider all solutions for which the posterior distribution on the cloud parallax of the farthest cloud passes this selection criterion to be valid. Any solution which does not satisfy this criterion is discarded in this first step. Finally, for each tested model j, we compute the Akaike Information Criterion (AIC) as (9)

based on the estimated maximum likelihood and the number of parameters in the model (M). We then compare the model performances by computing the probability (10)

that, among the tested models {m}, each model j is actually the one that minimizes the loss of information (Boisbunon et al. 2014), as in Pelgrims et al. (2023). The best model is the one with Pj|{m} = 1. We note that, according to the selection procedure explained above, we are giving some chances to solutions with large number of clouds (up to four) to be considered as the best model while the data itself might not be enough to place strong constraints on the farthest cloud. However, we checked all the solutions and found no evidence of spurious detections that might have resulted from bad fits.

A map of the number of clouds per LOS determined in this first step from the analysis of the hybrid-beam samples is shown in Fig. 11. For most of the surveyed area, we find evidence of two clouds per LOS. The maximum number of clouds is three, and a few sightlines show only one cloud.

thumbnail Fig. 9

Distribution of the number of stars in beam samples. (a) Number of stars per bin of distance modulus for the 287 sightlines sampling the observed regions when the beam geometry is hybrid. The colors indicate the number of sightlines for which a specific number of stars in a given distance bin is observed. The red-continuous line indicates the median number of stars per bin of distance for the entire set of sightlines. (b) Same as for (a) but for the conical-beam geometry. (c) Same as for (a) but for the cylindrical-beam geometry.

thumbnail Fig. 10

Number of stars per conical beam. The blue circle in the top right corner indicates the size of our conical beam with an angular radius of 13.74 arcmin. The magenta contour surrounds all HEALPix pixels whose center coincides with the center of the beams defining our star subsamples.

thumbnail Fig. 11

Map of the number of dust clouds identified at the end of Step 1 from the analysis of the hybrid-beam samples obtained for each beam centered on the pixel centers.

3.5.2 From posteriors to priors

We now aim at informing the maximum-likelihood analysis of data in conical beams from the results of the analysis of hybrid beam samples. To do so, we would ideally want to resample the posterior distributions obtained in Step 1, at least for the cloud-parallax parameters. This is however currently not possible using BISP -1 which does not make it possible to input prior samples as it relies on the Python nested sampling software dynesty which does not have this feature. For the purpose of the this work, we thus work around this limitation as explained below and leave the development of improved solutions for the future. We choose to ignore the possible correlations between model parameters and only extract (and propagate) information on cloud parallax and cloud mean-polarization properties from the marginalized posterior distributions of these parameters. By ignoring possible correlations, we know we are losing part of the information gained in Step 1.

For each of the valid models selected in Step 1 and their corresponding estimated (marginalized) posteriors, we define priors on the parameters of the models to be constrained by the data in Step 2. For the cases discussed above where the marginalized posterior distribution on the cloud parallax of the farthest cloud is multimodal, only the “valid” subset of the posterior samples is used to define the priors.

Neither a Gaussian nor a top-hat distribution generally represent the posterior distributions of the cloud parallax well; however, we find that in our case, they are sufficiently effective to impose constraints on the cloud’s parallax. We choose to construct Gaussian or top-hat priors on cloud parallaxes from the estimated posteriors as follows. We denote the posterior sample of parallax for a given cloud. For a given cloud, the estimated mean and standard deviation of the Gaussian prior, and the estimated minimum and maximum limits of the top-hat prior are obtained from percentiles of the posterior sample of the cloud parallax as (11)

where denotes the X-th percentile of the sample distribution. The prior distributions on the cloud parallax are then defined as follows: (12) (13)

for the Gaussian and the top-hat respectively.

For our data, we find that, even for the cases where the three-layer model is favored, only the two nearest clouds are located in the distance range where the beam geometry is cylindrical. Therefore, we do not need to update the priors of the most distant cloud parameters and simply left it as is for Step 1. That is, we inform the modeling of the data in conical beams (Step 2) from the fit in hybrid beams only for clouds that are found in the distance range where the hybrid and conical beams differ. If and if there are at least five stars with parallaxes lower than , then we adopt the Gaussian prior. In all other cases, we adopt a top-hat prior for this parameter. This is to avoid unphysical negative parallax values and to allow for the possibility of a farther away cloud.

Defining the priors on the cloud mean polarization from the corresponding posteriors is easier than for the parallax as the marginalized posteriors are generally close to Gaussian. To build the Gaussian priors on the mean Stokes parameters of the cloud, we consider the means of the posterior distribution samples and both the standard deviations of those samples (taken individually) and the mean of the posterior distribution of the corresponding element of the intrinsic-scatter covariance matrix. If are the estimated mean and standard deviation of the posterior distributions on qc (the mean q Stokes parameter of the cloud), and if is the estimated mean of the posterior on the qq element of the covariance matrix of the intrinsic scatter, then the Gaussian prior on qc is defined with a mean of and a standard deviation of . The same applies for uc· We include the terms from the intrinsic scatter in the definition of the priors on the mean Stokes parameters to account for the fact that, in the conical beams (higher angular resolution at lower distance), the “mean” polarization may, in general, pick a local fluctuation of the turbulent magnetized ISM defined at larger angular scales in the hybrid beams. Accordingly, we do not modify the priors for the elements of the intrinsic-scatter covariance matrix because, without further assumptions, we do not know how the intrinsic scatter evolves as a function of physical and angular scales. Consequently, we keep the same uniform priors defined above throughout the analysis.

thumbnail Fig. 12

Map of the number of dust clouds along the sightlines identified at the end of Step 2 in the conical beam centered on the pixel centers. The blue circle in the top-right corner indicates the extent of the conical beam.

3.5.3 Step 2 – LOS decomposition in conical beams

For each LOS, we have obtained a LOS decomposition in the hybrid beam, we have selected the best model, and we have defined priors from the estimated posterior distributions of this model. In this second step, we use BISP-1 to decompose the starlight polarization along distance for the star samples defined in the conical beams using the information gained above. For each LOS, we test the models with two, three, and four layers and use the informed priors for the two nearest clouds only. For the cases where the one-layer model was favored in Step 1, we test this model and a higher number of layers using the informed priors only on the nearest layer. In all cases we force the additional cloud(s) to be located at larger distances.

We use BISP-1 to run the nested sampling experiment using 1000 live points and sample the parameter space until an uncertainty of around 0.01 is achieved on the log of the model evidence. The required number of nested-sampling iterations was similar or higher than in Step 1.

When all models have been evaluated, we proceed as in Step 1 to inspect the solutions based on the marginalized posterior on the cloud parallax and keep only the valid reconstructions, and to select the best-model which minimizes the loss of information based on the AIC criterion. Figure 12 shows the number of clouds along each of the sightlines sampling the observed regions, as in Fig. 11. Again, the number of clouds ranges from one to three. A large fraction of the sightlines in this sky area intersects two clouds along the LOS. This map is further discussed in the Sect. 4.

The solutions of neighboring sightlines are not independent. First, the samples of stars from neighboring sightlines overlap, even for the conical beams, because we explicitly decided to oversample the observed region with a large number of sight-lines and the beam size is more than twice as large as the angular separation between adjacent sightlines. Second, the solutions obtained at the end of the second step are constrained by the larger angular scales at short distances since the priors are defined from the results obtained from the hybrid beam samples. As a result, we obtain the posterior distributions for our decomposition of the POS component of the magnetic field in dusty regions for a simply connected 3D volume without gaps.

3.6 Validation

To validate the result of our tomographic reconstruction, we rely on inspection of the significance of the polarization residuals of all the stars on which the reconstruction is based. We consider that we can be confident in our reconstruction if only a few data points show significant deviations from our model prediction, and if these data points are not clustered in the 3D volume. If, on the contrary, too large a fraction of stars show significant residuals, or if they are clustered in 3D space, this would indicate flaws or limitations in our tomographic reconstruction. In particular, the clustering of significant residuals would indicate the presence of features in the data that the model is not able to account for within uncertainties.

To estimate the significance of the residuals for all the stars taken individually, we compare the individual stellar observational data to the modeled data from the tomography results. We proceed as follows. Since we do not have reconstruction for each LOS toward each star individually but only for each LOS toward the center of HEALPix pixels at the center of our sample beams, we first identify the pixel in which the star falls in. We then read the posterior samples of the best-model decomposition obtained earlier for the corresponding beam.

By resampling the posterior distribution, we can then evaluate the expected polarization of the star at its given distance. According to our layer model (Pelgrims et al. 2023), the ISM contribution to a star polarization at distance d toward a given LOS is thus described by the stochastic model: (14)

with and where denotes a bivariate normal distribution with mean and covariance matrix Cint. The values of the cloud’s mean polarization and covariance matrix are obtained from the posterior samples and added, cloud wise, up to the distance of the star. The noise covariance matrix of the polarization measurements of the star can then be added to the intrinsic-scatter covariance matrix and the observation (s = (q u)) compared to the value predicted from our 3D reconstruction using the Mahalanobis distance: (15)

For each individual star, a single value of is obtained for each sample of the posteriors. By resampling the posteriors of the model parameters and resampling the parallax distribution of the star, we obtain a distribution of . This distribution informs us on the likelihood that the observed polarization is due to the dusty magnetized ISM given our 3D reconstruction. This estimate accounts for both the turbulence-induced scatter and observational noise in polarization and parallax. Large values of indicate significant residuals. In practice, a star is identified as an outlier (i.e., a star whose polarization is poorly predicted by our LOS model and thus showing significant residuals) if the median of its Mahalanobis-distance distribution is larger than a given threshold value. In 2D, the square of the Mahalanobis distances is expected to follow a χ2 distribution with two degrees of freedom. Accordingly, we can compute the threshold value corresponding to a given probability (Pth) to observe by chance a Mahalanobis distance greater than that. In 2D the threshold value is obtained as , where log is the natural logarithm.

Having obtained the Mahalanobis-distance distributions for all the stars in our sample, we look at the locations of stars with significant residuals in the 3D space. We show in Fig. 13 the sky locations of the stars for which the median of their Mahalanobis-distance distributions exceed threshold values corresponding to the probabilities of 5% and 1% of observing by chance greater values than that. The stars with significant residuals corresponding to lower probability of being compatible with the model are highlighted with blue and green circles, respectively.

The fractions of our star sample which show p-values lower than the 5% and 1% threshold are 3% and 0.3%, respectively, and no star shows a p-value lower than 0.02%. These values are further discussed in Sect. 5.5. We see from Fig. 13 that the stars with significant residuals do not cluster in particular places of the sky. Visually it seems that there is an excess of significant residuals in the eastern half of the observed region but this is due to the larger stellar density in that region. However, we observed (not shown on the figure) mild preference for the stars with significant residuals to be located at large distance. This is likely an effect of the low number density of stars at large distance and of the limited angular size of our sample beam as, for example, two of the four stars with their p-value lower than 1% are the farthest stars in their beam. Though, as the residuals are not very significant and that this trend is mild, we consider that this effect has no substantial effect on the present results.

This validation test makes it possible to verify the reliability of the assumptions underlying our modeling by looking at the spatial distribution of the polarization residuals. On the one hand, the violation of the thin-layer assumption would lead to systematic increases of the polarization residuals close to the reconstructed cloud distances. And, on the other hand, any significant variation of the polarization signal in the POS within the beam would lead to gradients in the residuals. We searched for such possible trends and could not find any. This suggests that our working assumptions are appropriate to model our dataset.

thumbnail Fig. 13

Sky map of residual significance. All the stars in the sample from which we performed our tomographic inversion are represented. Transparent gray dots (gray crosses) show stars that do (do not) fall in an HEALPix pixel for which we have tomography data. Blue and green circles show the stars for which the median of their distribution of Mahalanobis distances exceed a threshold values corresponding to the p-value indicated in the legend (see text). The lower the p-value, the more significant the residuals.

4 Results

The output of the 3D-inversion pipeline developed in the previous section consists of an ensemble of LOS decomposition of stellar polarization for non-independent samples. For each beam, we identified the number of components (shown in Fig. 12) and obtained the posterior distributions on their distances (cloud parallaxes, ) and polarization properties (qc, uc, Cint). From these, we can now infer the properties of the POS component of the magnetic field in dusty regions. We examine the results at the mean of the posterior distribution in Sect. 4.1 to infer the main features of the dusty magnetized ISM in the observed 3D volume. Then we build 3D maps of the posterior distributions in Sect. 4.2 and visualize the main output in Sect. 4.3. The discussion, interpretation, and validation of the results are provided in Sect. 5 along with caveats of the method.

thumbnail Fig. 14

Histogram of estimated mean cloud distances . Distribution of mean posterior cloud distances for the best fit models determined for all sightlines. The vertical lines indicate the limits to define the distance bins used in Sect. 4.1.

4.1 Basic exploration of the output

From the posterior distributions, we extract the estimated mean values for the cloud parallax and mean polarization for each LOS. In Fig. 14 we show the histogram of the (estimated) mean cloud distance for all clouds and all analyzed sightlines. This histogram shows two main separated peaks centered on 62 pc and 380 pc. A number of clouds are also found at distances larger than 1 kpc. Relying on this observation, we divide the distance axis in three bins for nearby (dc ≤ 265 pc), intermediate (dc ∈ [265, 650] pc), and distant (dc > 650 pc) clouds. For each distance bin, we generate maps of the cloud mean distance and mean polarization as shown in Figs. 15 and 16, respectively. The polarization maps are obtained by introducing the mean degree of polarization and mean polarization angle from the mean Stokes parameters as where we make use of the two-argument arctangent function to handle the π-ambiguity of the arctangent. In these maps, empty pixels of the observed regions mark the absence of detected clouds with mean distance in the specific range of distances. We find that all sight-lines intersect a cloud in the intermediate range of distances. A large fraction of the sightlines intersect both the nearby- and intermediate-distance components. As seen in Figs. 15 and 16, these components show a high degree of regularity in their distances and mean polarization properties, suggesting that each of these two components is related to physical entities (real interstellar clouds). The distances of detected clouds in the large distance bin are more scattered but we notice several clusterings of cloud detections in the 3D space for components with similar polarization properties.

Over the observed region of the sky, the mean estimated distance to the nearby cloud ranges from 30 pc to 170 pc and has a median of 58 pc. Its mean estimated degree of polarization ranges from 0.06% to 0.33% with a median value of 0.19%. The estimated POS orientation of the magnetic field of this component varies smoothly across the region. It ranges from 49.4° to 126.9° with a mean of about 80.5°.

The mean estimated distance to the intermediate cloud ranges from 354 pc to 448 pc with a median of 374 pc. Its mean estimated degree of polarization shows spatial variation over the observed region. It is also much higher than the nearby component, spanning the range from 0.68% to 2.03% with a median value of 1.45%. This is the dominant polarizing “screen” in the region, the effect of which was already clearly seen in Fig. 7. The mean estimated polarization angle of this component is nearly uniform throughout the observed region. It ranges from 28.4° to 51.9° and has a mean value of 39.5°. The standard deviation of all the polarization angles is 4.7°. This small dispersion indicates that the POS projection of the ordered component of the magnetic field in this cloud is nearly uniform despite the inhomogeneities seen in the mean degree of polarization.

At larger distances, we find several distinct clouds intersected by different sightlines. The distances of the distant clouds in neighboring pixels agree within uncertainties. Most noticeably, we find a cloud toward (l, b) ≈ (103.8°, 22.4°) with mean distance ranging from about 1600 pc to 1850 pc and which has a mean degree of polarization of about 0.26% and a mean polarization angle of about 136°. This is nearly perpendicular (≈85°) to the mean polarization orientation of the dominant polarization screen for the same sky pixels. This small region of the sky is slightly to the northwest of the “2-cloud LOS” studied in Panopoulou et al. (2019), Clark & Hensley (2019), Pelgrims et al. (2023) at (l, b) = (104.1°, 22.3°). In the east-southeast of this “2-cloud LOS", toward (l, b) ≈ (104.3°, 22.2°), we detect a cloud with mean distances ranging from 1700 pc to 2300 pc which has a mean degree of polarization of about 0.28% and a mean polarization angle of 66.2°, about 23° away from the mean polarization angle of the dominant polarization screen in the same sky pixels. In the southwestern part of the surveyed region, toward (l, b) ≈ (102.4°, 21.9°), we detect a cloud at distance between 1270 pc and 1500 pc. It has a mean degree of polarization of about 0.58% and a mean polarization angle of 16.5°, about 20° away from the mean polarization angle of the dominant polarization screen in the same sky pixels.

We also detect other groups of pixels with cloud detection close to the edges of the observed region. One at (l, b) ≈ (104.4°, 22.8°) with a somewhat large scatter in mean distance, ranging from 1760 pc to 2400 pc, and another one in the upper right corner of the observed region, with a mean distance of about 1170 pc. Other distant clouds are detected sparsely in the region. We have checked the tomographic decompositions of all these LOS individually, and none of these cloud detections appear to be the result of a statistical flaw in our analysis. We checked the posterior distributions of the different tested models to make sure that the automated model selection worked as expected. We additionally checked the residuals and visually compared the data and models in the (q, u) − µ space. In all cases, the best model was correctly identified.

We provide more discussion on our results and their reliability in Sect. 5. There we also relate them to known structures of the ISM that have already been studied in the literature, based on stellar extinction and H I data, in particular.

As a last step in the basic exploration of the output, we look at the covariance matrices that encode the intrinsic scatter of polarization properties within clouds (Cint). We evaluate the determinant of Cint for each sample of the posterior distributions, for each dust layer of the best model, and for each LOS individually. This quantity is related to the level of turbulence-induced intrinsic scatter within a cloud (see Appendix B of Pelgrims et al. 2023). We then estimate the median of the distributions of the determinants for each layer and each LOS. We finally build the distributions of the medians splitting the sample in terms of the cloud distances, as before. We show these distributions in Fig. 17 in the form of violin diagrams. We note that, while all the sightlines likely trace the same clouds in the nearby and intermediate distance ranges, this is unlikely to be the case for the distant distance range, as argued before. This is likely the reason for the apparent wider distribution of in this distance range than for the others. From Fig. 17, it is seen that, at least for the clouds at nearby and intermediate distances, the intrinsic scatter is detected above observational noise. This indicates that, in principle, subsequent analyses of the covariance matrices could lead to a detailed characterization of fluctuations in the magnetized ISM. We will explore such an avenue in future work.

thumbnail Fig. 15

Maps of the estimated mean cloud distances in the three bins of distances identified from Fig. 14, from far away (top) to nearby (bottom). The color scales span the full range of distance in each bin. Empty pixels of the observed regions mark the absence of clouds in the specific distance range. The magenta contour and the blue circle are as in Fig. 10.

thumbnail Fig. 16

Same as for Fig. 15 but for the cloud mean polarization. The color scales is identical for the three maps and informs on the mean degree of polarization . The red segments indicate the orientation of the mean polarization in the clouds which directly traces the orientation of the mean POS-component of the magnetic field.

thumbnail Fig. 17

Violin diagrams of the distributions of the square root of obtained from the posterior distributions and for the three distance ranges defined from Fig. 14. In each diagram, the horizontal segments indicate the 16th, 50th and 84th percentiles of the distribution.

4.2 3D-map making

In this subsection, we build 3D maps of the dusty magnetized ISM from the posterior distributions of all model parameters obtained for each LOS of the observed region. For each LOS, we construct the probability density function (PDF) of having a cloud at a given distance as follows. We stack the marginalized posterior distributions of the cloud parallax of each component, and obtain the distance distribution P(dc) by inverting every parallax value. The PDF is estimated from this distribution using a moving window of length ∆d as: (16)

By construction PDF(d) is normalized to the number of components intersected by the LOS. Then, using the estimated PDF(d), we construct the distance profiles of the different polarization properties by marginalizing the posterior distributions of each polarization parameter in distance bins. We thus estimate at any distance (d) the differential of any of the Stokes parameters or of the elements of the intrinsic-scatter covariance matrix as: (17)

where sc is any of qc, uc, Cint,qq, Cint,uu or Cint,qu and dc is the distance given by the inverse of the parallax. P(sc | d) is the conditional probability of having a value sc given a cloud distance dc at value d: (18)

where P(d, sC) is the 2D-marginalized posterior distribution between distance and the chosen polarization parameter sC, obtained by mapping the posterior in distance space. The units of the differentials of the Stokes parameters, δq and δu, are in polarization fraction per parsec, and the units of the differential of the intrinsic scatter covariance matrix, (δCint,xy, where the subscripts x and y denote either q or u) are polarization fraction per parsec squared. In the remainder of this paper, we focus on the mean properties of the dusty magnetized ISM that we can infer from stellar polarization and refer to future work for the study and characterization of its fluctuations.

Due to the limited number of samples used in estimating the posterior distributions, spurious noise is observed in the distance profiles. To reduce this, we smooth the estimated distance profiles along LOS using a Gaussian kernel. The posteriors are not sampled uniformly along LOS. This is because the sampling happens in the parallax space and also because the density of stars varies significantly with distance. It is not possible to choose a constant kernel value that would both smooth the noise at large distances and not severely dilute the signal at small distances. Therefore, we choose to smooth the profiles with a Gaussian kernel that varies with distance as (19)

where σk(d) and d are given in parsec. This choice is rather arbitrary, but it effectively smooths the noise at all distances. The choice for σk is such that it is close to 10 pc at a distance of 50 pc, it increases smoothly around 500 pc, and it reaches a value of about 50 pc for all distances larger than 1 kpc. The radial profiles and subsequent visualization are not strongly dependent on this choice. However, we recommend using the posterior distributions directly rather than the profiles for any subsequent quantitative analysis.

We show a set of such radial profiles of the Stokes parameters for 15 sightlines in Fig. 18 (middle and bottom panels). Nine sightlines are randomly chosen in the observed regions to which we add six sightlines that intersect three components as identified in Fig. 11 (those with b ≤ 22.5° and l ≥ 103.5°). Similar profiles can be constructed for the three parameters δCint,xy. The top panel of Fig. 18 shows the PDF of cloud distances for the same set of sightlines. To construct these profiles, we used a moving-window width of ∆d = 10 pc and sample the distance axis every parsec up to a distance of 3 kpc.

By repeating this process for all sightlines, we finally obtained the values of the differentials at any grid point sampling the 3D space corresponding to the observed sky region. These are the 3D maps that describe the dusty magnetized ISM in the observed sky region. Obtaining these maps is the main result of this paper. We present some visualizations of the maps in the next subsection. The 3D maps come naturally in a spherical coordinate system centered on the observer. From the 3D maps of the differential of the Stokes parameters, 3D maps of the polarization degree (δp) and of the polarization angle (δψ) can be derived. The 3D map of the polarization degree can be interpreted as a map of the density of polarizing material through space, weighted by the geometrical factor from the local inclination of the magnetic field on the LOS and by the local polarization efficiency. The 3D map of the polarization angle informs us on the orientation of the POS component (as seen from us, the observer) of the magnetic field locally. However, because the degree of polarization and the polarization angle are not additive quantities, as opposed to the Stokes parameters (and so of their differentials), the maps of δp and δψ must not be integrated along distance. In the first and second rows of Fig. 19, we show the cumulatives of the differentials δq and δu as a function of distance for the same sightlines as presented in Fig. 18. The cumulatives at a distance d correspond to the Stokes parameters that would be observed for a test star at that distance if we neglect the effects from turbulence. In the third and bottom rows we show the derived (mean) degree of polarization (p(d)) and (mean) polarization angle (ψ(d)) that would be observed for a star at the specific distance. These quantities are obtained at any distance d from the cumulatives of δq and δu, not from the cumulative of δp and δψ. The depolarization effect (decrease in p) due to the superposition of clouds with misaligned POS component of the magnetic field is clearly observed for some of the sightlines.

thumbnail Fig. 18

Variation along the distance of the normalized probability density distribution of the cloud distances (top) and of the differentials of the Stokes parameters δq and δu (middle and bottom) obtained from the marginalization of the full posterior distributions in distance bins for 15 sightlines as explained in the text. The moving-window length of 10 pc is adopted. The curves are smoothed with a Gaussian kernel that varies with distance as explained in the text. The same color corresponds to the same LOS across panels.

thumbnail Fig. 19

Polarization of a test star as a function of its distance for the same set of sightlines as shown in Fig. 18. The first and second rows show the cumulatives of the radial profiles of the differentials δq and δu as a function of distance. The third and fourth rows show the degree of polarization and polarization angle observed for the test star, as derived from the cumulatives of derived δq and δu. For p < 0.05%, ψ is masked (shown by thin gray line).

4.3 Visualization of the 3D maps

The portion of the Galactic space covered by our observations extends over ≈3.8 square degrees across the sky and extends up to 3.5 kpc. As a result, the geometry of the 3D maps that we construct is very elongated in radial distance and resembles a “pencil beam” having the geometry of a rectangular pyramid. This makes it difficult to visualize the results in 3D. To circumvent this limitation, we transformed the spherical coordinate system using the distance modulus as the radial coordinate as (20)

where µ(d) = 5(log10(d) − 1) is the distance modulus with d given in parsec and (l′, b′) are the angular coordinates (longitude and latitude) defined such that (l′, b′) = (0, 0) points toward the center of the observed region. In our case, we thus have Xµ(d), Yµ(d) ľ, and Z ≈ µ(d) b'. This coordinate system is ill-defined at small distances. However, we are not affected by this issue since all stars have a distance larger than 20 pc and that our 3D maps do not extend at distance modulus smaller than µ = 2.

We thus project the 3D maps of the differentials constructed above in this coordinate system. The sampling on angular coordinates is fixed by the centers of the HEALPix pixels used to define our samples from which the LOS decompositions were performed. The sampling on radial distance is such that we have values for the differentials at every parsec for the range µ ∈ [2, 12.7]. As defined, the X-axis runs through the center of the observed region at (l, b) = (103.3°, 22.3°), Y increases with increasing longitude and Z with latitude.

Despite the use of this coordinate system, the volume is still more extended toward the X-axis. For visualization purposes, we thus shrink that axis by a factor of 10 as compared to the others. We then construct a regular Cartesian grid made of 2563 voxels, with limits such that the volume X/10 × Y × Z fits in. Due to the rectangular pyramid geometry of the inverted dataset, most of the volume is empty (no data). We use linear interpolation to obtain the values of the differentials at every voxel position from the 3D maps of the differentials in the modified spherical coordinate system. The data cubes can now be visualized.

4.3.1 Plane projections

We start by producing plane projections of the data cubes. The data cubes are integrated along the Y-axis and the Z-axis of the Cartesian grid to produce the vertical and horizontal plane projections, respectively. Figure 20 shows the results of such projection for the differentials of the (mean) Stokes parameters. In this figure, and as explained in its caption, we tweak the color scales in order to visualize the very faint features at large X values. The faint signal of the distant components seen in Fig. 18 at µ ≈ 10 are indeed dimmed because of their small extent in space and by the integration over the full length of the data cube axes. The nearby and dominant components already seen in Fig. 18 are striking on these plane projections. The dominant component appears nearly planar at constant X (distance) and with very coherent polarization properties. This is reminiscent of what we already discussed in Sect. 4.1. The change of colors along distance indicates that the POS component of the magnetic field in the nearby, intermediate, and large distance ranges are misaligned.

Plane-projection maps reveal features that are somewhat elongated along the distance axis. This is clearly visible at small distances. Most of it is due to the use of the moving window of 10 pc and the smoothing that we use to construct the distance profiles of the differentials and due to the choice of the logarithmic scale (distance modulus) for visualization. However, it is worth emphasizing that the constructed 3D maps result from the marginalization of the posterior distributions along the distance axis. Therefore, some real extensions along the LOS exist in the maps and are related to the level of constraints we can impose on cloud distances given the stellar data. They do not inform the actual extension of a dust cloud along the LOS but rather reflect our uncertainties on cloud distances. These features are similar to the “finger-of-god effect” commonly seen in 3D dust map reconstructions. We recall that the model that we use to reconstruct the dusty magnetized ISM from stellar data in polarization and distance assumes that the clouds are thin - they have no dimension along the LOS.

thumbnail Fig. 20

Plane projections of δq (left) and δu (right) in the vertical (top) and horizontal (bottom) planes. The vertical and horizontal planes are defined by averaging the data cubes (rectangular pyramid) in the Cartesian grid along the Y-axis and the Z-axis, respectively. The observer (us) is on the left and distance increases to the right. In order to visualize faint features, the color values (c) are obtained from the projected values (v) as · The units of the color scale are therefore . The color scales are symmetrical about zero and the range of δu is twice as large as that of δq, reflecting the difference in magnitude of the two quantities that is observed in Fig. 18.

4.3.2 Visualization in 3D

Visualizing a pseudo-vector field in 3D presents more challenges than visualizing a scalar field, such as the dust density distribution for example. While preparing this paper, we have started addressing these challenges building up on the framework underlying ASTERION, the tool developed in Konstantinou et al. (2022) to simulate magnetized dust clouds in 3D. ASTERION relies on real-time 3D visualization techniques with virtual reality capabilities. It is implemented in the realtime-engine called Unreal Engine 4 (UE4)6 and allows for the rendering of details of the magnetized ISM and enables the user to fly through the simulated environment as is done in video games.

Extending the capabilities of this software, we are now able to visualize the main properties of our tomography map. Yet, the current version of the software does not render the contribution from the intrinsic scatter nor does it visualize the uncertainties on our 3D reconstruction. Only the solution corresponding to the means of our posterior distributions (or at the maximum likelihood values) on the cloud distances and their mean polarization can be visualized. Specifically, the software makes it possible to place the magnetized dust cloud in 3D and to visualize the “differentials” of δpC and δψC using colors, transparency and a 3D version of the line integral convolution technique (Cabral & Leedom 1993). The visualization of the 3D map obtained in this paper can be accessed online7.

5 Discussion

In this section, we compare the results of our 3D reconstruction of the POS component of the magnetized ISM from starlight polarization and distance data to other datasets that inform on the complexity of the ISM. We also discuss the limitations of our results and caveats of our method.

5.1 Comparison with 3D dust extinction data

We start with a comparison of our results with a 3D dust extinction map that was obtained from Gaia parallax and from both spectroscopic and photometric extinction measurements, including from Gaia and 2MASS, and which informs on the dust density distribution in 3D space in a Cartesian box of 3 × 3 × 0.8 kpc3 centered on the Sun (Vergely et al. 2022). Our starlight-polarization-based tomography map is independent from the 3D extinction map as it relies on different observables (aside from the parallax) and as the stellar polarization sample is (much) smaller than the stellar extinction sample. Both 3D maps inform us on different properties of the ISM (extinction versus polarization). In addition, the 3D map of dust extinction and our map of stellar polarization have different spatial resolutions both along the LOS and in the POS. The spatial resolution of the 3D dust map that we consider is 10 pc (and sampled by 3D voxels with side length of 5 pc). Our map has a varying spatial resolution, first because the LOS inversion is not bound to a particular sampling (or gridding) of the distance, and second because the signal along each LOS has been smoothed with a distance-dependent Gaussian kernel. Despite these differences, we wish to verify whether our approach identifies clouds at distances similar to the cloud locations in the 3D dust map. We thus construct 1D profiles of extinction as a function of distance to compare the locations of peaks in the dust distribution along the LOS with the locations of clouds in our tomography map.

For this first qualitative comparison between dust extinction and dust polarization tomography data, we sample both 3D maps with a set of sightlines that go through the surveyed sky area and set by the HEALPix tessellation with Nside = 128. For each LOS we extract the distance profiles of differential extinction and of the “differential” of the degree of polarization (δpC), computed locally from the differentials of the Stokes parameters (δqC and δuC). We present the profiles in Fig. 21 where we also indicate with a vertical strip the estimate of the distance to the inner surface of the Local Bubble as derived in Pelgrims et al. (2020)8 from the 3D extinction map of Lallement et al. (2019). Because of the finite size of the 3D map of Vergely et al. (2022), the profiles do not extend over the whole distance range for which we recover information from stellar polarization. Comparing the profiles from dust extinction and dust polarization tomography data, we notice that both tracers indicate the presence of a very nearby component at µ < 5 and a dominant component at µ ≈ 7.9.

It is remarkable that the stellar polarization data makes it possible to recover the very nearby component at µ ≲ 5 although it appears to be very shallow in the extinction profiles. The difference between extinction and polarization data in relative amplitudes between the nearby component and the dominant one could point to differences in magnetic field inclination or dust polarization properties. However, differences in 3D-map resolutions could artificially dilute more the signal of the nearby component in 3D dust density map than in polarization. Thus, this comparison will require confirmation from a dedicated analysis which goes beyond the scope of this paper. Meanwhile, it should be noted that our polarization tomography data independently confirms the presence of a very close dust cloud identified in the 3D dust extinction map allowing us to argue that these small features in the 3D dust map are real (at least in the surveyed region). Evidence for nearby ISM clouds within the Local Bubble are also provided by pulsar-scintillation studies (Ocker et al. 2024 and references therein) and stellar line absorption (Peek et al. 2011).

The agreement between distances to the dominant peak seen in the and δpC profiles at µ ≈ 7.9 is remarkable. However, while this component appears as an isolated an narrow peak in the δpC profiles, it appears broader toward lower distances, possibly featuring a second peak, in the profiles. According to Pelgrims et al. (2020), the nearby peak would correspond to the wall of the Local Bubble. The main peak (centered in µ ≈ 7.9) could correspond to the diffuse Cepheus Flare as several molecular clouds in this part of the sky show similar distances (see Schlafly et al. 2014).

However, it is unclear whether the broadening of the peak in the profiles actually indicates the presence of two close dust components, or if this shape simply stems from the fact that there is too little stellar extinction data right in front of the dominant component to effectively constrain the lower distance limit of a single cloud centered at the main peak. We recall that A' profiles, like our distance profiles in the cloud parameters, reflect the shape of the posterior distribution, which is closely related to the uneven distribution of stars along the distance axis. Additionally, we note that we have identified the presence of a cloud with distance µ ≈ 7, several degrees away in the southeastern direction, outside of the surveyed region, and that the outskirt of this cloud is intersected by our region. This is observed both in 3D maps from Vergely et al. (2022) and Edenhofer et al. (2024). It is therefore possible that the peak broadening in the profiles is due to the limited resolution of the 3D dust density maps, or that the broadening indicates a real dust overdensity. In the latter scenario, we understand why we are unable to find evidence for the existence of this cloud using starlight polarization by the following. If there are two dust components in the range µ ∈ [6, 9], they are relatively close in distance. Thus the number of stars with polarization measurements that would be useful to identify a cloud and constrain its polarization properties in between the two peaks is very low. We have checked that, indeed, we have few polarization measurements in that distance range in our sample. Furthermore, the amplitude of the jump in degree of polarization that is induced by the dominant screen is very large. This hampers the possible identification of a counterpart of the wall of the Local Bubble in its vicinity as it would require a large number of data points. This illustrates one of the limitations in our reconstruction which likely comes from the limited depth of our survey. An alternative explanation would be that the nearest of these two clouds simply does not induce polarization. This could happen if this cloud is devoid of polarizing dust grains or if both its column density is very low and its permeating magnetic field lines are nearly perpendicular to the POS.

One major difference between our approach and the approach of Lallement et al. (2019) and Vergely et al. (2022) is that they use a fixed spatial kernel in the 3D space to infer the dust density in a grid while the distance to the cloud is a free parameter in our model. Together with the thin-layer assumption in our model, this might explain the tighter constraints on cloud distances that we obtained (378 pc ± 14 pc), compared to the fixed maximum resolution of 10 pc in the 3D dust density map of (Vergely et al. 2022).

Overall, the agreement between our tomography map from stellar polarization with the tomography map of dust density from stellar extinction is remarkable given the very different nature of the two datasets. We find good agreement in both the number of components and their distances to the Sun. This indicates that a systematic and dedicated analysis will help study spatial variations of the ISM properties such as the polarization efficiency of the dust grains and the inclination of the magnetic field lines with respect to the LOS. However, our qualitative comparison taught us that comparing tomography data in polarization and density obtained from different approaches (as it is the case here) and with different resolutions could lead to confusion. To benefit from the comparison of dust extinction and polarization data, it is essential to perform detailed analysis treating the different datasets in a self-consistent manner. We leave the task of devising such a framework for future work.

thumbnail Fig. 21

Qualitative comparison of distance profiles between dust polarization and dust extinction tomography results. Top: profiles of the “differential” of the degree of polarization (δpC, computed locally from δqC and δuC) as a function of distance modulus as inferred from our polarization 3D map. Bottom: differential extinction as a function of distance modulus as inferred from the 3D extinction map of Vergely et al. (2022). No extinction data is available at distance larger than about 1 kpc. The surveyed sky area is sampled according to an HEALPix map with Nside = 128. Each LOS is represented with a different color, the same color is used across panels. The green vertical strip indicates the range of distances to the inner surface of the Local Bubble in this sky area (Pelgrims et al. 2020).

5.2 Detection of distant clouds and HI data

In Fig. 21, we observe farther away clouds in the δpC profiles where no information exists in profiles, although indications for faraway clouds might be guessed in the far edge of the extinction profiles (at µ ≈ 10). To corroborate our findings of faraway clouds (≳ 1 kpc), we first turn to the inspection of H I velocity spectra. As discussed in Sect. 2, we expect to observe complex H I spectra with features at intermediate velocities if faraway dust clouds are present along these sightlines.

For this purpose, we extract the H I velocity spectra from the HI4PI data cube (HI4PI Collaboration 2016) in pixels corresponding to the sightlines for which we detected a cloud with mean distance in the large distance bin (dC > 650 pc) in Sect. 4.1. These sightlines can be seen in the top panel of Fig. 16, for example. We sorted the sightlines according to their Galactic longitude and present the spectra in Fig. 22 from blue (low longitude) to red (high longitude). We exclude the range of high-velocity clouds (vLSR ≲ −90 km s−1) and velocities above 40 km s−1 where there is no signal.

The H I velocity spectra clearly show complexity with several components. All sightlines at large Galactic longitude show strong IVC components (|vLSR| ≳ 35 km s−1). These IVC components do not vanish at lower longitude but are shallower. The fact that there is power in the velocity spectra for |vLSR| > 35 km s−1 for all these sightlines makes us confident in our detection of faraway clouds. However, the evidence for IVC components at low Galactic longitudes is milder than for high longitudes; the column density of these clouds is much lower in this part of the surveyed region, as already inferred from Fig. 2. The spectra corresponding to the upper right part of the region, with l < 103° and b > 22.1°, are the only ones not showing local maxima at vLSR < −20 km s−1. At this stage we do not know if this is an indication for spurious detection of faraway clouds in our starlight-polarization based tomography. We discuss further the reliability of our cloud detections in Sect. 5.3.

In their analysis of a small subset of the polarization data studied in this paper, Panopoulou et al. (2019) studied starlight polarization within a beam of about 9.6 arcmin radius toward (l, b) = (104.1°, 22.3°) and made the identification of a faraway dust cloud with an IVC with velocity of about −50 km s−1. Having enlarged the surveyed area with stellar polarization, we see that IVCs are detected in this part of the sky using starlight polarization, although not specifically in the pixel closest to the LOS studied by Panopoulou et al. (2019), but in neighboring pixels (eastward and westward). It is not clear from our polarization tomography data whether these faraway clouds form a continuous entity in 3D space or if they are independent. Assuming a common distance of 1700 pc, the detected clouds are at least 15 pc apart. The orientations of the POS component of the magnetic field that we find are also different in the three main clusters of pixels with faraway-cloud detection. A different POS magnetic field orientation is however compatible with the finding of Clark & Hensley (2019) who used the orientation of H I fibers to infer this quantity. This is illustrated in Fig. 23 where we overlaid our results from starlight-polarization-based tomography for the clouds in the large distance bin (as in the top panel of Fig. 16) to the H I-orientation data from Clark & Hensley (2019). We also obtain consistent results if we instead use the H I orientation maps from Halal et al. (2023). The H I data is integrated over the velocity range vLSR ∈ [−60.45, −44.99] km s−1. The POS component of the magnetic field implied by the orientation of H I fibers is shown using line-integral-convolution (LIC) texture. The alignment between H I structures and the magnetic field is driven by the orientation of anisotropic, cold gas structures (Clark et al. 2019; Peek & Clark 2019; Murray et al. 2020; Kalberla et al. 2020; Kalberla 2023), so the LIC pattern may not trace the magnetic field well away from regions of prominent H I emission. An IVC component, with varying POS magnetic field component, is prominent in the upper left corner of the surveyed area. Visual comparison of the orientation of the POS component of the magnetic field obtained from starlight-polarization tomography and H I data reveals a good qualitative agreement where faraway clouds are detected, in particular toward the IVC. However, we notice that we do not detect faraway components toward all sightlines sampling the prominent IVC, especially in the closest pixel toward the “2-cloud” region of Panopoulou et al. (2019), and that our cloud detections in the western part of the region do not have counterparts in this velocity range. We further discuss these differences in the next subsection.

Overall, our comparison between starlight-polarization and H I data points toward the great potential of combining the different tracers of the magnetized ISM to study and characterize its properties in detail. Going far beyond the scope of this paper, we will explore the potential of such a synergy in future work. Meanwhile, all of the above indicates that our inversion method for obtaining, from stellar polarization and parallax alone, a tomo-graphic view of the POS component of the magnetic field in dusty regions leads to reliable results.

thumbnail Fig. 22

H I velocity spectra for sightlines with far away clouds. The spectra are sorted according to the Galactic longitudes of their sightlines and shown with color ranging from blue to red. For visual perception, the velocity spectra are smoothed with a Gaussian kernel with width of 1.9 km s−1.

thumbnail Fig. 23

H I-orientation data integrated in the velocity range vLSR ∈ [−60.45, −44.99] km s−1 toward our surveyed area. The background color shows the H I intensity on a linear stretch, and the texture represents the POS magnetic field orientation inferred from HI fibers. The red segments indicate the mean POS magnetic field orientation obtained from starlight polarization tomography for clouds with distances larger than 1 kpc, as in Fig. 16. The white-dotted circle shows the “2-cloud region” studied in Panopoulou et al. (2019). The magenta contour and the blue circle are as in Fig. 10.

5.3 Reliability of cloud detection and limitations

To present our results, in Sect. 4, and to compare them with other probes of the ISM, in the two subsections above, we focus on the tomographic solutions corresponding to the best model (number of clouds) obtained per LOS. However, it is interesting to look at the second-best model (defined below) and to compare the performance between the best and second-best models that we obtained. In fact this helps quantify the robustness of a solution (a given model) and, also, the reliability of a cloud detection.

In Sect. 3.5.3, at the end of Step 2 in the inversion process, we obtained for each LOS, the estimated maximum-likelihood values for all tested models (different number of layers). Among the tested models, the best model was identified comparing the model performances based on their AIC values (see Eqs. (9) and (10)). By ranking the AIC values, we can also identify the second-best model. The probability that this model is actually the model that minimizes the loss of information as compared to the best model (i.e., its Pj|{m} value) quantifies by how much it is outperformed by the best model.

In the top panel of Fig. 24, we show the map of the number of clouds per LOS corresponding to the second-best model. In the bottom panel of the same figure, we show the map of the second-best model probabilities. In these maps, white pixels (enclosed in the magenta outline) indicate that no second-best model can be identified. That is, there is no “valid” reconstruction with a different number of clouds along the LOS. In practice, this happens for several sightlines where the best model is a 2-layer model, and for which the 3-layer and 4-layer models led to bad posterior distributions (see Sect. 3.5.1). This means that either there is no additional cloud farther away than the second cloud, or the data is not good enough to make it possible to detect it solely based on stellar data. We recall that when a 2-layer model is selected as the best model in Step 1, the 1-layer model is not considered in Step 2.

The comparison of the number of clouds per LOS from the best model (Fig. 12) and from the second-best model (top panel of Fig. 24), taking into account its probability, suggests that for most of the sightlines where the 3-layer (2-layer) model is favored by the data, we cannot safely ignore the 2-layer (3-layer) model. Assuming the presence of faraway clouds in the eastern (left) part of the region, as motivated by H I data, this indicates that the data is not enough to allow for a strong and robust cloud detection at large distances, suggesting that we reach the limit of the cloud-detection capability given the data. In the western (right) part of the region, the probability for the 2-layer model is generally high for the 3-layer sightlines, suggesting either unreliable (spurious) detection of faraway cloud or, again, that we reach the limit of the cloud-detection capability.

The cloud-detection capability is primarily limited by the number of measurements and their uncertainties as compared to the strength of the polarization signal induced by a cloud to its background stars (see Sect. 4 in Pelgrims et al. 2023). This shows that the reliability of our tomographic results at large distances is limited by the current depth of our survey and that more polarization data is needed to enable a stronger model separation. However, given the complementary nature of stellar extinction, H I and starlight polarization data, which we also highlighted above, there is the possibility to incorporate such external data in the inversion process. In principle this would help differentiate between competing models, even in the absence of additional stellar polarization data. We will undertake this endeavor in future work.

The depth of the survey, and therefore the number density of the measurements, also sets a limit on the maximum angular resolution at which we can achieve the LOS inversion of starlight polarization data. The angular resolution is directly linked to one of the main limitations of using BISP-1. The method does not explicitly take into account POS variations of the polarization signal within the beam other than through the intrinsic scatter term. We use a fixed geometry to define our beam and consider a top-hat acceptance window to include stars in our sample. As a result, while modeling the ISM along distance, stars at the edge of the beam (on the POS) contribute the same as stars at the center of the beam. This is fine as long as the polarization signal does not vary much on the POS and can be described by our model. However, artifacts (such as overestima-tion of the intrinsic-scatter covariance matrix and biases in the mean polarization properties) may be expected as soon as the polarization data cannot be well described by a bivariate normal distribution. This may happen in the presence of substantial POS variations at the angular scale comparable to the size of our beam or if only part of the beam intersects a cloud. Such shortcomings in polarization data modeling can naturally hamper cloud detection.

This limitation likely explains differences between the starlight-polarization and H I tomographic data at large distances seen in Fig. 24, in particular toward the prominent IVC. Specifically, H I data suggests substantial variation of the signal in our beam of 13.76 arcmin radius, and in particular, an abrupt change of the POS component of the magnetic field toward the “2-cloud region” of Panopoulou et al. (2019). We notice that part of these limitations could be overcome by introducing weights on the polarization data while computing the log-likelihood to account for the angular distance of each star with respect to the center of the beam, that is, of the LOS. We will test this idea in future work.

thumbnail Fig. 24

Maps of the second-best model: (top) number of clouds and (middle) probability for the second-best model to be the model that minimizes the loss of information against the best model. The bottom panel combines the information from the two upper panels. It indicates the number of clouds (color) whose opacity is given by the probability. Transparent pixels correspond to low probabilities for second-best models.

5.4 Caveats of the inversion method

The main caveats of our approach to obtain a 3D map of the POS component of the magnetic field in dusty regions come from the use of our Bayesian method (BISP-1) which works along the LOS. This method requires the definition of beams within which the starlight polarization data is “averaged” on the POS and modeled as a bivariate normal distribution with a mean and covariance matrix according to the dust-layer model we rely on.

The first limitation of the designed approach is that it uses a fixed beam geometry and does not explicitly correlate solutions for different beams. We address this problem by choosing to oversample the sky with non-independent beam samples and by applying our decomposition method to each of them. However, while the solutions for overlapping beam samples are not independent, the correlation is not quantified nor controlled through, for example, the use of density power spectrum as it is the case for 3D dust density mapping (e.g., Green et al. 2019; Lallement et al. 2019, 2022; Leike & Enßlin 2019; Vergely et al. 2022; Edenhofer et al. 2024). Introducing such a correlation for polarization tomography is not a trivial task and would require assumptions that we wish to avoid at this stage, in particular because of the intertwined nature of matter and magnetic field, and of the inherent degeneracy between the density of polarizing dust and the inclination of magnetic field with respect to the LOS. A careful analysis of our tomography results in combination with simulation-based studies will help shed light on how to implement such correlations for the case of dust polarization. We will address this question in future work.

A second, important limitation of BISP-1 is that it does not explicitly account for POS variations of the polarization signal. For the diffuse sightlines targeted in this work, significant POS variations within the beam do not seem to be present, as our modeling provides a good description of the data. As mentioned in Sect. 3.6, we searched for possible systematic variations of the polarization residuals within our beam and could not find any. This suggests that any POS variation of the polarization signal in our beam are successfully characterized by the intrinsic-scatter covariance matrix. However, we expect more significant POS variations to arise toward denser regions, for example toward nearby molecular clouds, where the column density varies by an order of magnitude within tens of arcsec-onds (e.g., filaments measured by Herschel, André et al. 2010). Some nearby molecular clouds have been targeted with deep polarimetric surveys (e.g., Pereyra & Magalhães 2004; Santos et al. 2017), approaching a density of 1000 measurements per square degree. In such regions, a careful examination of the choice of beam size and the POS variations of column density would be necessary when using BISP-1 to decompose polarization along the LOS, as envisioned by e.g. Soler et al. (2016). Given the existing stellar polarization data in the literature, it appears that the dataset presented here is particularly favorable for applying BISP-1 in a moving-window scan scheme due to its combination of high number density of stars and the fact that it probes diffuse sightlines, conditions that will be encountered by the PASIPHAE survey targeting high and intermediate Galactic latitudes.

5.5 Astrophysical use of the output

Mapping continuously the stellar-polarization source field in 3D opens the way to tackle several science objectives that were thus far out of reach or left to the study of specific clouds or LOS. A direct use of our results, which does not require postprocessing of the tomography map, is the production of a list of intrinsically polarized-star candidates. Other possible uses, which will however require further postprocessing and specific analysis, are mentioned in Sect. 6.

Using the same procedure as in Sect. 3.6, we proceed to the estimation of the significance of the residuals in polarization for all the stars making our polarization samples, with successful Gαiα cross-match and which satisfy the quality criterion on parallax estimate (RUWE ≤ 1.4, see Sect. 2.3). This sample is made of 1448 stars among which 18 were identified as possible outliers according to the recursive sigma-clipping approach employed in Sect. 3.3. Among the 1448 stars, 1392 fall in an HEALPix pixel for which we have a model for the magnetized ISM along distance. For each of them we obtain a distribution of the Mahalanobis distance values informing us on the likelihood that the measured polarization is compatible with our picture of the dusty magnetized ISM, taking into account all sources of uncertainties and scatter in model and observations. We show the histogram of the median of the distributions for the full sample in Fig. 25 where we separate the stars flagged as outliers in Sect. 3.3 from the others. This figure confirms that most of the outliers discarded from the tomography analysis indeed show polarization properties that are not compatible with the picture of the dusty magnetized ISM that we reconstructed. In some sense, this also further validates our 3D reconstruction.

Using our tomography map to estimate the likelihood that the polarization of any star falling in the reconstructed 3D volume is solely due to the dusty magnetized ISM provides us with a robust way to identify outlier candidates and therefore intrinsically polarized-star candidates. We thus add the values of the median of the distribution that we obtained for each star to our published catalog (see Table 1). The higher this value, the more likely the target is to be an outlier. This may be used to plan follow-up observations to study these sources in more detail. Only additional study will confirm whether these targets are real outlying data points or if they merely pick up fluctuations of the magnetized ISM that are unaccounted for in our reconstruction. Our analysis also confirms that the fraction of intrinsically polarized stars in the ISM is rather low, at least at these Galactic latitudes. Only 14 stars out of 1448 have a p-value lower than 0.2% for their polarization to be induced by the ISM only. That is, only about 1 % of our sample may be made of intrinsically polarized stars whereas we did not try to minimize this fraction while planning our observation, for example based on stellar type.

We must note a possible intrinsic degeneracy in our procedure. As mentioned in Sect. 3.3, during the sigma-clipping procedure, it is possible to inadvertently exclude from the analysis data points that merely capture fluctuations in the magnetized ISM. In this case, a too low level of intrinsic scatter would be obtained from the fit and, accordingly, any discarded points would show a large value in the a posteriori test described above. The list of outlier candidates thus depends on the choice of the hyper parameters of the sigma-clipping procedure, and more specifically on the choice of the used significance threshold. This illustrates the relevance of follow-up observations and analysis of the outlier candidates. However, if the turbulence-induced intrinsic scatter had to be systematically underestimated in our reconstruction, we would obtain a distribution of values that would be statistically shifted toward large values and that would not correspond to a bivariate normal distribution of the polarization residuals. This is not what we obtain and what is also shown in Fig. 25. On the contrary, we obtain a distribution of values that is slightly shifted to lower values than expected and which, therefore, suggests a small overestimation of the covariance matrices. We understand the latter as coming from unaccounted for variation of the polarization signal in the POS within our beam.

thumbnail Fig. 25

Histogram of the median of the per star distributions. The histograms for the “ISM probes” and “outliers” identified in Sect. 3.3 are separated (in blue and red) and stacked on top of each other. Notice the logarithmic scale of the vertical axis. The dashed vertical lines correspond to the p-value thresholds of 5%, 1%, and 0.2%. Most of the outliers show significant residuals.

6 Summary and concluding remarks

In this work, we performed a survey of optical starlight polarization for a continuous region covering about four square degrees centered on (l, b) = (122°, 33°), and designed a pipeline to obtain the first 3D map of the dusty magnetized ISM based on a Bayesian analysis of the measurements in stellar polarization and distances only. Obtaining this map, which corresponds to an extended volume in 3D space, is the main result of this paper. Our reconstruction corresponds to a sky area of about 3.8 square degrees and extends up to 3 kpc from the Sun. We found that the 3D volume covered by our tomography data is populated by several clouds and that a large fraction of the sightlines in the surveyed regions intersect at least two clouds, one being very close to the Sun with a distance of about 62 pc and a dominant polarizing screen at about 375 pc. Distant clouds are also detected up to a distance of about 2 kpc. We were able to corroborate our findings using a 3D dust extinction map and H I-velocity spectra. We are thus confident that our inversion pipeline works and that stellar data in polarization and distance alone are a powerful probe of the magnetized ISM. Specifically, for the diffuse ISM where dust grains are expected to align their shortest axes with the ambient magnetic field lines, starlight polarization allows us to determine locally, in 3D space, the orientation of the POS component of the magnetic field (the position angle) and the amplitude of the starlight-polarization source field (the degree of polarization) which depends on the local dust density, dust grain polarization efficiency, and on the inclination of the magnetic field lines with respect to the sightlines.

We obtained our polarization tomography map by adopting a moving-window strategy to scan the surveyed region with non-independent beams and making use of the LOS-inversion method implemented in BISP-1. This allowed us to invert the data and reconstruct the dusty magnetized ISM in 3D while keeping the number of assumptions to its minimum. Namely, we relied on the thin-dust layer model developed in Pelgrims et al. (2023) and assumed that its assumptions are valid for the adopted beam size of 13.74 arcmin radius. We expect that the analysis of the obtained 3D map will enable the unbiased study and characterization of the properties of the dusty magnetized ISM in 3D, such as through 3D correlation functions. This will enable us in the future to develop 3D inversion methods capable of taking into account variations and correlations in space, both along the distance and in the POS.

Finally, we expect that polarization tomography maps, as the one we obtained in this work and that provides local measurements of the POS component of the magnetic field and individual cloud polarization properties, will enable breakthroughs in the modeling of the Galactic magnetic field, in the modeling and characterization of the dusty magnetized ISM as a contaminant foreground in observations of the cosmic microwave background polarization, and in the modeling of astrophysical dust. These research goals, along with the estimation of the strength of the magnetic field through the quantification of the variation of the POS component of the magnetic field (e.g., Skalidis & Tassis 2021; Skalidis et al. 2021), will necessitate dedicated analyses and postprocessing of our 3D maps, such as to obtain proper boundaries of dust clouds. In future works we will explore these research directions based on the 3D map that we have presented here. This will set the stage for future analyses that will greatly benefit from the most awaited polarization data from the PASIPHAE survey (Tassis et al. 2018).

Acknowledgements

We thank an anonymous referee for a thorough and constructive report that helped us improve the clarity of the paper. The PASIPHAE program is supported by grants from the European Research Council (ERC) under grant agreements No. 7712821 and No. 772253; by the National Science Foundation (NSF) award AST-2109127; by the National Research Foundation of South Africa under the National Equipment Programme; by the Stavros Niarchos Foundation under grant PASIPHAE; and by the Infosys Foundation. This work was also partly supported by the ERC grant agreements No. 819478. V.P. acknowledges funding from a Marie Curie Action of the European Union (grant agreement No. 101107047). V.P. also thanks Philipp Frank and Sebastian Hutschenreuter for the fruitful and inspiring discussions related to this work during the program “Toward a Comprehensive Model of the Galactic Magnetic Field” at Nordita in April 2023, which was partly supported by NordForsk and Royal Astronomical Society. V.Pa. acknowledges support by the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the “First Call for H.F.R.I. Research Projects to support Faculty members and Researchers and the procurement of high-cost research equipment grant” (Project 1552 CIRCE). V.Pa. also acknowledges support from the Foundation of Research and Technology – Hellas Synergy Grants Program through project MagMASim, jointly implemented by the Institute of Astrophysics and the Institute of Applied and Computational Mathematics. K.T. and A.P. acknowledge support from the Foundation of Research and Technology – Hellas Synergy Grants Program through project POLAR, jointly implemented by the Institute of Astrophysics and the Institute of Computer Science. T.G. is grateful to the Inter-University Centre for Astronomy and Astrophysics (IUCAA), Pune, India for providing the Asso-ciateship programme under which a part of this work was carried out. This work was supported by NSF grant AST-2109127 and was carried out in part at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. This research has used data, tools or materials developed as part of the EXPLORE project that has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 101004214.

References

  1. Alves, J., Zucker, C., Goodman, A. A., et al. 2020, Nature, 578, 237 [Google Scholar]
  2. Andersson, B. G., Lazarian, A., & Vaillancourt, J. E. 2015, ARA&A, 53, 501 [Google Scholar]
  3. André, P., Men’shchikov, A., Bontemps, S., et al. 2010, A&A, 518, A102 [Google Scholar]
  4. Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147 [Google Scholar]
  5. Berdyugin, A., Piirola, V., & Teerikorpi, P. 2014, A&A, 561, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Bialy, S., Zucker, C., Goodman, A., et al. 2021, ApJ, 919, L5 [NASA ADS] [CrossRef] [Google Scholar]
  7. Blinov, D., Kiehlmann, S., Pavlidou, V., et al. 2021, MNRAS, 501, 3715 [NASA ADS] [CrossRef] [Google Scholar]
  8. Blinov, D., Maharana, S., Bouzelou, F., et al. 2023, A&A, 677, A144 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Boisbunon, A., Canu, S., Fourdrinier, D., Strawderman, W., & Wells, M. T. 2014, Int. Stat. Rev., 82, 422 [CrossRef] [Google Scholar]
  10. Boulanger, F., Enßlin, T., Fletcher, A., et al. 2018, JCAP, 2018, 049 [CrossRef] [Google Scholar]
  11. Cabral, B., & Leedom, L. 1993, in SIGGRAPH'93 Proc. 20th Annual Conf., Computer Graphics and Interactive Techniques (New York, USA: Association for Computing Machinery), 263 [Google Scholar]
  12. Clark, S. E., & Hensley, B. S. 2019, ApJ, 887, 136 [NASA ADS] [CrossRef] [Google Scholar]
  13. Clark, S. E., Peek, J. E. G., & Miville-Deschênes, M. A. 2019, ApJ, 874, 171 [NASA ADS] [CrossRef] [Google Scholar]
  14. Clemens, D. P., Cashman, L. R., Cerny, C., et al. 2020, ApJS, 249, 23 [NASA ADS] [CrossRef] [Google Scholar]
  15. Cotton, D. V., Bailey, J., Kedziora-Chudczer, L., et al. 2016, MNRAS, 455, 1607 [NASA ADS] [CrossRef] [Google Scholar]
  16. Covino, S., Smette, A., & Snik, F. 2020, https://doi.org/10.5281/zenodo.3906379 [Google Scholar]
  17. Das, K. K., Zucker, C., Speagle, J. S., et al. 2020, MNRAS, 498, 5863 [Google Scholar]
  18. Davis, Leverett, J., & Greenstein, J. L. 1951, ApJ, 114, 206 [NASA ADS] [CrossRef] [Google Scholar]
  19. Doi, Y., Nakamura, K., Kawabata, K. S., et al. 2023, ApJ, accepted [arXiv:2311.13054] [Google Scholar]
  20. Draine, B. T., & Hensley, B. S. 2021, ApJ, 919, 65 [NASA ADS] [CrossRef] [Google Scholar]
  21. Edenhofer, G., Zucker, C., Frank, P., et al. 2024, A&A, in press, https://doi.org/10.1051/0004-6361/202347628 [Google Scholar]
  22. Ellis, R. S., & Axon, D. J. 1978, Ap&SS, 54, 425 [NASA ADS] [CrossRef] [Google Scholar]
  23. Fadeyev, Y. A. 2007, Astron. Lett., 33, 103 [NASA ADS] [CrossRef] [Google Scholar]
  24. Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Gaia Collaboration (Vallenari, A, et al.) 2023, A&A, 674, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Gontcharov, G. A., & Mosenkov, A. V. 2019, MNRAS, 483, 299 [Google Scholar]
  28. Goodman, A. A., Bastien, P., Myers, P. C., & Menard, F. 1990, ApJ, 359, 363 [NASA ADS] [CrossRef] [Google Scholar]
  29. Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [Google Scholar]
  30. Green, G. M., Schlafly, E., Zucker, C., Speagle, J. S., & Finkbeiner, D. 2019, ApJ, 887, 93 [NASA ADS] [CrossRef] [Google Scholar]
  31. Großschedl, J. E., Alves, J., Meingast, S., et al. 2018, A&A, 619, A106 [Google Scholar]
  32. Halal, G., Clark, S. E., Cukierman, A., Beck, D., & Kuo, C.-L. 2023, arXiv e-prints [arXiv:2306.10107] [Google Scholar]
  33. Hall, J. S. 1949, Science, 109, 166 [NASA ADS] [CrossRef] [Google Scholar]
  34. Heiles, C. 1976, ARA&A, 14, 1 [NASA ADS] [CrossRef] [Google Scholar]
  35. Heiles, C. 1996, ApJ, 462, 316 [NASA ADS] [CrossRef] [Google Scholar]
  36. Hensley, B. S., & Draine, B. T. 2021, ApJ, 906, 73 [NASA ADS] [CrossRef] [Google Scholar]
  37. Heyer, M., Gong, H., Ostriker, E., & Brunt, C. 2008, ApJ, 680, 420 [Google Scholar]
  38. HI4PI Collaboration (Ben Bekhti, N., et al.) 2016, A&A, 594, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Hiltner, W. A. 1949, Science, 109, 165 [NASA ADS] [CrossRef] [Google Scholar]
  40. Hutsemékers, D. 1998, A&A, 332, 410 [Google Scholar]
  41. Ivanova, A., Lallement, R., Vergely, J. L., & Hottier, C. 2021, A&A, 652, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Kalberla, P. M. W. 2023, arXiv e-prints [arXiv:2401.00190] [Google Scholar]
  43. Kalberla, P. M. W., Kerp, J., & Haud, U. 2020, A&A, 639, A26 [EDP Sciences] [Google Scholar]
  44. King, O. G., Blinov, D., Ramaprakash, A. N., et al. 2014, MNRAS, 442, 1706 [NASA ADS] [CrossRef] [Google Scholar]
  45. Konstantinou, A., Pelgrims, V., Fuchs, F., & Tassis, K. 2022, A&A, 663, A175 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Kuntz, K. D., & Danly, L. 1996, ApJ, 457, 703 [NASA ADS] [CrossRef] [Google Scholar]
  47. Lallement, R., Babusiaux, C., Vergely, J. L., et al. 2019, A&A, 625, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Lallement, R., Vergely, J. L., Babusiaux, C., & Cox, N. L. J. 2022, A&A, 661, A147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Leike, R. H., & Enßlin, T. A. 2019, A&A, 631, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Leike, R. H., Glatzle, M., & Enßlin, T. A. 2020, A&A, 639, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Li, H.-b., Fang, M., Henning, T., & Kainulainen, J. 2013, MNRAS, 436, 3707 [NASA ADS] [CrossRef] [Google Scholar]
  52. Lindegren, L., Hernandez, J., Bombrun, A., et al. 2018, A&A, 616, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Lindegren, L., Klioner, S. A., Hernandez, J., et al. 2021, A&A, 649, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Magalhães, A. M., Pereyra, A., Melgarejo, R., et al. 2005, Astronomical Polarimetry: Current Status and Future Directions, eds. A. Adamson, C. Aspin, C. Davis, & T. Fujiyoshi, ASP Conf. Ser., 343, 305 [Google Scholar]
  55. Magalhães, A. M., de Oliveira, C. M., Carciofi, A., et al. 2012, Stellar Polarime-try: from Birth to Death, eds. J. L. Hoffman, J. Bjorkman, & B. Whitney, AIP Conf. Ser., 1429, 244 [Google Scholar]
  56. Magkos, G., & Pavlidou, V. 2019, JCAP, 2019, 004 [CrossRef] [Google Scholar]
  57. Maharana, S., Kypriotakis, J. A., Ramaprakash, A. N., et al. 2020, SPIE Conf. Ser., 11447, 114475E [Google Scholar]
  58. Marchal, A., & Martin, P. G. 2023, ApJ, 942, 70 [NASA ADS] [CrossRef] [Google Scholar]
  59. Martin, P. G. 1974, ApJ, 187, 461 [NASA ADS] [CrossRef] [Google Scholar]
  60. Martínez-Solaeche, G., Karakci, A., & Delabrouille, J. 2018, MNRAS, 476, 1310 [Google Scholar]
  61. McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565 [Google Scholar]
  62. Monet, D. G., Levine, S. E., Canzian, B., et al. 2003, AJ, 125, 984 [Google Scholar]
  63. Mouschovias, T. C., Tassis, K., & Kunz, M. W. 2006, ApJ, 646, 1043 [Google Scholar]
  64. Murray, C. E., Peek, J. E. G., & Kim, C.-G. 2020, ApJ, 899, 15 [CrossRef] [Google Scholar]
  65. Nishiyama, S., Hatano, H., Tamura, M., et al. 2010, ApJ, 722, L23 [NASA ADS] [CrossRef] [Google Scholar]
  66. O’Callaghan, M., Gilmore, G., & Mandel, K. S. 2023, arXiv e-prints [arXiv: 2311.03199] [Google Scholar]
  67. Ocker, S. K., Cordes, J. M., Chatterjee, S., et al. 2024, MNRAS, 527, 7568 [Google Scholar]
  68. Panopoulou, G., Tassis, K., Blinov, D., et al. 2015, MNRAS, 452, 715 [Google Scholar]
  69. Panopoulou, G. V., Tassis, K., Skalidis, R., et al. 2019, ApJ, 872, 56 [NASA ADS] [CrossRef] [Google Scholar]
  70. Panopoulou, G. V., Dickinson, C., Readhead, A. C. S., Pearson, T. J., & Peel, M. W. 2021, ApJ, 922, 210 [CrossRef] [Google Scholar]
  71. Panopoulou, G. V., Markopoulioti, L., Bouzelou, F., et al. 2023, AASJ, submitted [arXiv:2307.05752] [Google Scholar]
  72. Patat, F., Maund, J. R., Benetti, S., et al. 2010, A&A, 510, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  73. Peek, J. E. G., & Clark, S. E. 2019, ApJ, 886, L13 [NASA ADS] [CrossRef] [Google Scholar]
  74. Peek, J. E. G., Heiles, C., Peek, K. M. G., Meyer, D. M., & Lauroesch, J. T. 2011, ApJ, 735, 129 [NASA ADS] [CrossRef] [Google Scholar]
  75. Pelgrims, V., Ferrière, K., Boulanger, F., Lallement, R., & Montier, L. 2020, A&A, 636, A17 [EDP Sciences] [Google Scholar]
  76. Pelgrims, V., Clark, S. E., Hensley, B. S., et al. 2021, A&A, 647, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  77. Pelgrims, V., Panopoulou, G. V., Tassis, K., et al. 2023, A&A, 670, A164 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  78. Pereyra, A., & Magalhães, A. M. 2004, ApJ, 603, 584 [CrossRef] [Google Scholar]
  79. Planck Collaboration Int. XIX. 2015, A&A, 576, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  80. Planck Collaboration XII. 2020, A&A, 641, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Ramaprakash, A. N., Rajarshi, C. V., Das, H. K., et al. 2019, MNRAS, 485, 2355 [NASA ADS] [CrossRef] [Google Scholar]
  82. Santos, F. P., Ade, P. A. R., Angilè, F. E., et al. 2017, ApJ, 837, 161 [NASA ADS] [CrossRef] [Google Scholar]
  83. Schlafly, E. F., Green, G., Finkbeiner, D. P., et al. 2014, ApJ, 786, 29 [NASA ADS] [CrossRef] [Google Scholar]
  84. Skalidis, R., & Tassis, K. 2021, A&A, 647, A186 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  85. Skalidis, R., Panopoulou, G. V., Tassis, K., et al. 2018, A&A, 616, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  86. Skalidis, R., Sternberg, J., Beattie, J. R., Pavlidou, V., & Tassis, K. 2021, A&A, 656, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. Skilling, J. 2004, Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, eds. R. Fischer, R. Preuss, & U. V. Toussaint, AIP Conf. Ser., 735, 395 [NASA ADS] [CrossRef] [Google Scholar]
  88. Soler, J. D., Alves, F., Boulanger, F., et al. 2016, A&A, 596, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  89. Speagle, J. S. 2020, MNRAS, 493, 3132 [Google Scholar]
  90. Spoelstra, T. A. T. 1972, A&A, 21, 61 [NASA ADS] [Google Scholar]
  91. Tahani, M., Lupypciw, W., Glover, J., et al. 2022, A&A, 660, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  92. Tassis, K., & Pavlidou, V. 2015, MNRAS, 451, L90 [NASA ADS] [CrossRef] [Google Scholar]
  93. Tassis, K., Ramaprakash, A. N., Readhead, A. C. S., et al. 2018, arXiv e-prints [arXiv:1810.05652] [Google Scholar]
  94. Tegmark, M., & de Oliveira-Costa, A. 2001, Phys. Rev. D, 64, 063001 [Google Scholar]
  95. Tsouros, A., Edenhofer, G., Enßlin, T., Mastorakis, M., & Pavlidou, V. 2024, A&A, 681, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  96. Vergely, J. L., Lallement, R., & Cox, N. L. J. 2022, A&A, 664, A174 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  97. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
  98. Zonca, A., Singer, L., Lenz, D., et al. 2019, J. Open Source Softw., 4, 1298 [Google Scholar]
  99. Zucker, C., Goodman, A., Alves, J., et al. 2021, ApJ, 919, 35 [NASA ADS] [CrossRef] [Google Scholar]
  100. Zucker, C., Goodman, A. A., Alves, J., et al. 2022, Nature, 601, 334 [NASA ADS] [CrossRef] [Google Scholar]

2

The limited astrometric accuracy of RoboPol is primarily due to uncertainties in the joint modeling of the 4-spot pattern and field distortions. New generation wide-field polarimeters such as the Wide Area Linear Optical Polarimeters (e.g., Maharana et al. 2020) will enable improved astrometric accuracy with wider coverage and the absence of the 4-spot pattern.

3

In the USNO-B1.0 catalog, blended objects have a single identifier and the magnitude of the star represents the total photometry of both stars.

4

The polarization as a function of distance is shown in Sect. 3.3 and we study it starting from Sect. 3.5.

6

UE4 is a complete suite of development tools that allows for the visualization and immersive virtual worlds, multiplatform deployment, asset and plugin marketplace, among other features (https://www.unrealengine.com/en-US/unreal).

All Tables

Table 1

Polarization catalog (abbreviated).

All Figures

thumbnail Fig. 1

Sky location of the surveyed area of about four square degrees. Left: full-sky map of the dust emission as seen by Planck at 353 GHz (Planck Collaboration XII 2020). The color represents the intensity of dust emission on logarithmic scale. The line-integral-convolution texture shows the polarization angle of the dust emission rotated by 90°. Middle: a zoom-in of the map toward the surveyed regions, which includes part of the North Celestial Pole Loop on the East of the map. Right: a closer view of the surveyed region. Black segments indicate the polarization orientation from the stars in our survey and from Panopoulou et al. (2019). The segments are scaled according to the polarization fraction. Unpolarized stars appear as dots. The blue, horizontal segment in the bottom right corner shows the scale for a 1% polarized star. Outlier candidates (see Sect. 3.3) are not shown.

In the text
thumbnail Fig. 2

H I velocity data from the HI4PI survey in the sky area toward the region surveyed in starlight polarization. The top panels show the brightness temperature as a function of the gas velocity measured in the LSR. The blue spectrum corresponds to the velocity spectrum averaged over the full region. Two dominant peaks with velocities in the LVC (≈ − 1 km s−1) and IVC (≈ − 50 km s−1) ranges are seen. The bottom panels show the column density maps resulting from the integration of the velocity spectra in the ranges depicted by shaded regions marked in the top panels in the LVC (left) and IVC (right) ranges. The maps reveal a high degree of complexity in morphological structures. The magenta outline indicates the sky region of 3.8 sq. deg. for which starlight-polarization-based tomography is obtained in this work (see Sect. 3.4).

In the text
thumbnail Fig. 3

Sky distribution (top) and distance (modulus) distribution obtained by inverting the parallax (bottom) of our stellar sample with reliable parallax estimates. The sample is divided according the USNO-B R-band magnitude. Blue, red and black correspond to stars with R < 14 (shallow survey), R ≥ 14 (deep survey) and unknown values due to identification mismatch. The histograms on the bottom panel are stacked.

In the text
thumbnail Fig. 4

2D distribution of the number of stars in our catalog in the parallax – S/N plane. For better visualization, a dozen stars with large parallax values and or large parallax S/N are not included in this plot.

In the text
thumbnail Fig. 5

Distribution of the polarization properties for the stellar sample with a reliable parallax estimate. Distributions of the Stokes parameters (q’s and u’s) and of the logarithm of their uncertainties (log10(σq) and log10(σu)) are shown in the top and bottom panels, respectively. The horizontal and vertical dashed lines are used for visual reference. They indicate the q = 0 and u = 0 loci in the top panel and indicate the values of 0.1% and 1% in polarization uncertainties in the bottom panel.

In the text
thumbnail Fig. 6

Flowchart for the 3D inversion pipeline, from the stellar catalog with reliable parallax estimates (top) to the 3D maps of differential quantities (bottom).

In the text
thumbnail Fig. 7

Stellar relative Stokes parameters (q, u) versus the distance modulus for the full sample. The qV and uV (in the Equatorial coordinate system) are shown by green circles and blue diamonds, respectively. The vertical error bars indicate the uncertainties (noise and systematic) on the Stokes parameters and the horizontal error bars represent the uncertainties on the star distance modulus converted from the 1σ parallax uncertainties. Outliers identified from the iterative sigma-clipping approach in groups of nearest neighbors are highlighted with red-filled circles for q and purple-filled squares for u. To facilitate the visual clarity of the plot we restrict the range of (q, u) and µ values. Ten stars have µ < 5 and 30 stars have µ > 14.

In the text
thumbnail Fig. 8

Beam sizes as a function of distance from the observer. Top: angular aperture radius of our hybrid beam (continuous gray) resulting from the combination of a conical beam (red dashed) and cylindrical beam (dot-dashed). Bottom: extent of the beams in the transverse direction to the LOS. Same line convention as in the top panel. In this example the cylinder radius is fixed at 2.5 pc and the angular radius of the cone at 13.74 arcmin. For the hybrid beams, the cylindrical geometry prevails at distances smaller than ≈630 pc, while the conical geometry prevails at larger distances.

In the text
thumbnail Fig. 9

Distribution of the number of stars in beam samples. (a) Number of stars per bin of distance modulus for the 287 sightlines sampling the observed regions when the beam geometry is hybrid. The colors indicate the number of sightlines for which a specific number of stars in a given distance bin is observed. The red-continuous line indicates the median number of stars per bin of distance for the entire set of sightlines. (b) Same as for (a) but for the conical-beam geometry. (c) Same as for (a) but for the cylindrical-beam geometry.

In the text
thumbnail Fig. 10

Number of stars per conical beam. The blue circle in the top right corner indicates the size of our conical beam with an angular radius of 13.74 arcmin. The magenta contour surrounds all HEALPix pixels whose center coincides with the center of the beams defining our star subsamples.

In the text
thumbnail Fig. 11

Map of the number of dust clouds identified at the end of Step 1 from the analysis of the hybrid-beam samples obtained for each beam centered on the pixel centers.

In the text
thumbnail Fig. 12

Map of the number of dust clouds along the sightlines identified at the end of Step 2 in the conical beam centered on the pixel centers. The blue circle in the top-right corner indicates the extent of the conical beam.

In the text
thumbnail Fig. 13

Sky map of residual significance. All the stars in the sample from which we performed our tomographic inversion are represented. Transparent gray dots (gray crosses) show stars that do (do not) fall in an HEALPix pixel for which we have tomography data. Blue and green circles show the stars for which the median of their distribution of Mahalanobis distances exceed a threshold values corresponding to the p-value indicated in the legend (see text). The lower the p-value, the more significant the residuals.

In the text
thumbnail Fig. 14

Histogram of estimated mean cloud distances . Distribution of mean posterior cloud distances for the best fit models determined for all sightlines. The vertical lines indicate the limits to define the distance bins used in Sect. 4.1.

In the text
thumbnail Fig. 15

Maps of the estimated mean cloud distances in the three bins of distances identified from Fig. 14, from far away (top) to nearby (bottom). The color scales span the full range of distance in each bin. Empty pixels of the observed regions mark the absence of clouds in the specific distance range. The magenta contour and the blue circle are as in Fig. 10.

In the text
thumbnail Fig. 16

Same as for Fig. 15 but for the cloud mean polarization. The color scales is identical for the three maps and informs on the mean degree of polarization . The red segments indicate the orientation of the mean polarization in the clouds which directly traces the orientation of the mean POS-component of the magnetic field.

In the text
thumbnail Fig. 17

Violin diagrams of the distributions of the square root of obtained from the posterior distributions and for the three distance ranges defined from Fig. 14. In each diagram, the horizontal segments indicate the 16th, 50th and 84th percentiles of the distribution.

In the text
thumbnail Fig. 18

Variation along the distance of the normalized probability density distribution of the cloud distances (top) and of the differentials of the Stokes parameters δq and δu (middle and bottom) obtained from the marginalization of the full posterior distributions in distance bins for 15 sightlines as explained in the text. The moving-window length of 10 pc is adopted. The curves are smoothed with a Gaussian kernel that varies with distance as explained in the text. The same color corresponds to the same LOS across panels.

In the text
thumbnail Fig. 19

Polarization of a test star as a function of its distance for the same set of sightlines as shown in Fig. 18. The first and second rows show the cumulatives of the radial profiles of the differentials δq and δu as a function of distance. The third and fourth rows show the degree of polarization and polarization angle observed for the test star, as derived from the cumulatives of derived δq and δu. For p < 0.05%, ψ is masked (shown by thin gray line).

In the text
thumbnail Fig. 20

Plane projections of δq (left) and δu (right) in the vertical (top) and horizontal (bottom) planes. The vertical and horizontal planes are defined by averaging the data cubes (rectangular pyramid) in the Cartesian grid along the Y-axis and the Z-axis, respectively. The observer (us) is on the left and distance increases to the right. In order to visualize faint features, the color values (c) are obtained from the projected values (v) as · The units of the color scale are therefore . The color scales are symmetrical about zero and the range of δu is twice as large as that of δq, reflecting the difference in magnitude of the two quantities that is observed in Fig. 18.

In the text
thumbnail Fig. 21

Qualitative comparison of distance profiles between dust polarization and dust extinction tomography results. Top: profiles of the “differential” of the degree of polarization (δpC, computed locally from δqC and δuC) as a function of distance modulus as inferred from our polarization 3D map. Bottom: differential extinction as a function of distance modulus as inferred from the 3D extinction map of Vergely et al. (2022). No extinction data is available at distance larger than about 1 kpc. The surveyed sky area is sampled according to an HEALPix map with Nside = 128. Each LOS is represented with a different color, the same color is used across panels. The green vertical strip indicates the range of distances to the inner surface of the Local Bubble in this sky area (Pelgrims et al. 2020).

In the text
thumbnail Fig. 22

H I velocity spectra for sightlines with far away clouds. The spectra are sorted according to the Galactic longitudes of their sightlines and shown with color ranging from blue to red. For visual perception, the velocity spectra are smoothed with a Gaussian kernel with width of 1.9 km s−1.

In the text
thumbnail Fig. 23

H I-orientation data integrated in the velocity range vLSR ∈ [−60.45, −44.99] km s−1 toward our surveyed area. The background color shows the H I intensity on a linear stretch, and the texture represents the POS magnetic field orientation inferred from HI fibers. The red segments indicate the mean POS magnetic field orientation obtained from starlight polarization tomography for clouds with distances larger than 1 kpc, as in Fig. 16. The white-dotted circle shows the “2-cloud region” studied in Panopoulou et al. (2019). The magenta contour and the blue circle are as in Fig. 10.

In the text
thumbnail Fig. 24

Maps of the second-best model: (top) number of clouds and (middle) probability for the second-best model to be the model that minimizes the loss of information against the best model. The bottom panel combines the information from the two upper panels. It indicates the number of clouds (color) whose opacity is given by the probability. Transparent pixels correspond to low probabilities for second-best models.

In the text
thumbnail Fig. 25

Histogram of the median of the per star distributions. The histograms for the “ISM probes” and “outliers” identified in Sect. 3.3 are separated (in blue and red) and stacked on top of each other. Notice the logarithmic scale of the vertical axis. The dashed vertical lines correspond to the p-value thresholds of 5%, 1%, and 0.2%. Most of the outliers show significant residuals.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.