© ESO 2018

1 Introduction

Galactic dust is ubiquitous throughout the sky (e.g., Planck Collaboration XI 2014) and interacts with the large-scale Galactic magnetic field. Asymmetric dust grains tend to orient their short axis along the magnetic field lines. The most plausible mechanism of alignment is given by radiative alignment torque theory (RAT; for a recent review on grain alignment, see Andersson et al. 2015). As a result of this alignment, the dust thermal emission is polarized perpendicular to this axis (Stein 1966; Cudlip et al. 1982) at far-infrared (FIR). On the other hand, starlight that passes through a dusty region suffers dichroic extinction; this results in the starlight becoming polarized parallel to the field lines (Hiltner 1949; Hall 1949; Davis & Greenstein 1951). Therefore, the plane-of-sky component of the magnetic field can be traced through the polarization of starlight caused by dust.

The fractional linear polarization, p, is related to the dust column density, and therefore to stellar reddening, E(B–V). Observed values of p are bounded by the empirical upper limit p_max = 9(%)E(B−V) (Hiltner 1956; Serkowski et al. 1975). The majority of existing optical polarization observations have been driven by star formation studies, and consequently are agglomerated near the Galactic plane (e.g., Heiles 2000), that is, at regions with high E(B–V). There have been, however, a number of works targeting stars that have low E(B–V), either because they lie in our local neighborhood (Piirola 1977; Tinbergen 1982; Leroy 1993; Bailey et al. 2010) or because they are located at high Galactic latitudes (e.g., Appenzeller 1968; Berdyugin et al. 2014, and references therein). The samples of stars in these surveys have been selected on the basis of distance, and consist entirely of bright stars (V <13 mag for the deepest sample of Berdyugin et al. 2014, which extends out to ≤600 pc from the Sun). These works have measured the level of interstellar polarization towards individual stars that are spread out over several thousand square degrees. Though informative, these sparsely sampled (in all three dimensions) datasets form an incomplete picture of interstellar polarization at low extinctions.

The interest in understanding interstellar medium (ISM) polarization in this low dust-column regime is multifaceted. There is much to be gained in terms of understanding of the Galactic magnetic field and its effect on the diffuse ISM. In addition, this regime can offer new insights regarding the micro-physical interaction of dust with the magnetic field. A third and largely sought-after reward relates to the role of Galactic dust as a foreground to the cosmic microwave background (CMB).

Polarized thermal dust emission from our Galaxy is a major obstacle in the search for the primordial B-mode signal in the polarization of the CMB (BICEP2/Keck Collaboration 2015). This signal is predicted to arise from the effect of gravitational waves on the last scattering surface after the inflationary epoch (Seljak 1997; Seljak & Zaldarriaga 1997; Kamionkowski et al. 1997a,b; Zaldarriaga & Seljak 1997).

The problem of foreground subtraction is challenging and previously used methods have been proven largely incomplete. One common way of treating the contamination problem has been to extrapolate the signal from frequencies where dust emission dominates (~350 GHz) to frequencies where CMB emission is important (~60–150 GHz). However, this extrapolation can become problematic. For example, in the presence of two clouds along the same line-of-sight, the polarization at one frequency could be decorrelated compared to that at another frequency (Tassis & Pavlidou 2015). The conditions for this decorrelation to be significant are as follows: (1) the magnetic field of one cloudis significantly misaligned with that of the other and (2) the temperatures of the clouds are not identical.

Knowledge of the orientation of the magnetic field on the plane of the sky as a function of distance from the observer is necessary to address this effect. Thermal emission cannot provide this information because its intensity is the result of integration along the line-of-sight to infinity, and therefore distance information is lost. Stellar polarization, on the other hand, only traces the dust column out to the distance of the star. With enough measurements of stars of known distances tracing a similar sightline, one can reconstruct the plane-of-sky magnetic field orientation as a function of distance. Existing stellar polarization measurements are very sparse at high-latitudes, which are the regions targeted by CMB experiments. The optical polarizationsurvey PASIPHAE¹ is being designed to address precisely this issue.

In this work, we conduct a path-finding mini-survey for PASIPHAE in three regions with very low dust content. We wish to determine the level of polarization that can be measured in regions with very low dust emission using a flux-limited sample of stars located within a very small area (~0.05 square degrees). In contrast to previous works at high latitude, this approach allows for determination of the average polarization properties of the ISM locally.

The surveyed regions and the observing strategy are presented in Sect. 2. A description of the reduction and instrument calibration are given in Sect. 3, supplemented by Appendix A. Our results are presented in Sect. 4. Our findings are discussed in Sect. 5 and our conclusions in Sect. 6.

Fig. 1 Mollweide projection of the E(B–V) map of Lenz et al. (2017; centered at l, b = [0,0]; grid-line spacing is 30°). The red stars mark our target regions. Gray areas are not included in the map.

2 Observations

Observations were conducted with the 1.3-m telescope at Skinakas Observatory in Crete, Greece² using the RoboPol polarimeter (King et al. 2014). RoboPol is a four-channel imaging polarimeter, designed to simultaneously measure the relative Stokes parameters q and u. Each source in the 13′   × 13′ field of view (FOV) is projected on four locations on the CCD. The central part of the FOV (2′   × 2′) is shadowed by a focal plane mask whose purpose is to lower the background for the central target. For this project, we measure sources only in this masked region, to maximize measurement precision. We conducted all observations in the R-band.

As stellar polarization depends on reddening, we find the mean E(B–V) towards the observed regions. We refer to the target regions as Dark Patches (DPs). To this end, we use the Lenz et al. (2017) reddening (hereafter LHD) map (Fig. 1). This map was derived from HI emission using the HI4PI all-sky survey (HI4PI Collaboration 2016) and covers 39% of the sky with HI column densities N_HI < 4 × 10²⁰ cm⁻². The map isprovided in HEALPIX format with NSIDE = 1024 (pixel spacing 3.3′). Figure 2 shows the distribution of E(B–V) in the LHD map. The resolution of this map (16′) is comparable to the size of our regions. Therefore, we assign a single value of E(B–V) to each DP: that of the E(B–V) map at the center of the region. The reddening of the target regions is shown with vertical lines at 0.0063 mag (solid line, DP1), 0.0072 mag (dotted line, DP2) and 0.0118 mag (dash-dotted line, DP3). We compare these values with an independent estimate of E(B–V) from the map of Planck Collaboration Int. XXIX (2016). We assume R_v =3.1 (Schultz & Wiemer 1975) to convert the A_v to E(B–V). Remarkably (considering all the factors that contribute to the uncertainty), we find very small differences between these values and those of LHD. According to Planck, E(B–V) is 0.0083 mag in DP1, 0.0077 mag in DP2, and 0.0171 mag in DP3. In the E(B–V) map of Fig. 1, 24, 94 and 1751 square degrees in total correspond to regions with E(B–V) lower than that ofDP1, DP2, and DP3, respectively.

The centers in Galactic coordinates (l,b) and angular sizes of the three regions are: DP1; (124.7, 60.0), 15′ × 15′; DP2; (159.4, 49.0), 16′ × 16′; and DP3; (191.1, 48.6), 13′ × 13′. The locations of the three surveyed regions are marked with red star symbols on the LHD map in Fig. 1.

Within each target region, we constructed flux-limited samples. Stars with 10.5 mag < R < 16.5 mag were selected from the USNO-B catalog (Monet et al. 2003). We discarded stars that would suffer from confusion with nearby sources. Due to observing-time limitations, our final samples are complete to 16.47 mag for DP1, 15.7 mag for DP2, and 16.25 mag for DP3 (with an additional star of 16.6 mag).

Observations took place in the period May–August 2015 for DP1, September–November 2015 for DP3, and October–November 2017 for DP2. The total observing time for the three regions was 27.5 hr for DP1, 16 hr for DP2, and 20 hr for DP3. In total we observed 68 stars from which 24 were in DP1, 23 in DP2, and 21 in DP3. For each star the exposure time was selected so that the statistical (photon-noise) error of p would be comparable to the instrumental systematic error of 0.1% (see Sect. 3). Total exposure times for targets ranged from 2 to 165 min, with ~80% of targets having an exposure time of <50 min. Polarization standard stars were observed each night for calibration. During the 2015 observing run, two standards were observed each night on most nights, and one standard on some nights, due to visibility/time constraints. During the 2017 observing run, from two to six standards were observed during each night, with the exception of one night when only one was observed.

Fig. 2 Distribution of E(B–V) from the map of Lenz et al. (2017). Vertical lines correspond to the E(B–V) of the target regions DP1, DP2, and DP3. The truncation at 0.045 mag corresponds to the mask placed by Lenz et al. (2017) at hydrogen column densities NHI > 4 × 1020 cm−2.

3 Data reduction

The data were reduced using the RoboPol pipeline (King et al. 2014), which performs aperture photometry of each source to measure the relative Stokes parameters q and u and their (statistical) uncertainties σ_q and σ_u, respectively. These are used to calculate the fractional linear polarization, p, and the electric vector position angle (EVPA or χ) through $$p=\sqrt{q^2+u^2},{\sigma }_p=\sqrt{\frac{q^2{\sigma }_q^2+u^2{\sigma }_u^2}{q^2+u^2}},$$(1) $$\chi =\frac{1}{2}\arctan \left(\frac{u}{q}\right),{\sigma }_{\chi }=\frac{1}{2}\sqrt{\frac{u^2{\sigma }_q^2+q^2{\sigma }_u^2}{{(q^2+u^2)}^2}}.$$(2)

We use the latest version of the pipeline which selects the aperture size according to the optimization method presented in Panopoulou et al. (2015).

Two modifications have been made after the publication of that paper. First, in Panopoulou et al. (2015), the optimization was run separately on each of the four images of a target. However, this could potentially introduce artificial differences between the photometry of the stellar images, leading to low levels of spurious polarization (which could be significant for this work). To alleviate this, we use an aperture that is common for the pair of ordinary and extraordinary stellar images. In practice, we run the optimization on each stellar image separately, and then use an aperture size that is the mean of the two solutions for the pair of images used to calculate one Stokes parameter and likewise for the remaining pair of images. Second, we recalibrate the aperture optimization for sources in the RoboPol mask. The optimal aperture fora stellar image is found by constructing its growth curve (intensity as a function of aperture size, x), fitting a polynomial f(x) and solving the equation: df∕dx = λf(x). The parameter λ is calibrated using standard stars. A value of 0.02 was found for sources in the RoboPol field in Panopoulou et al. (2015). For this work, we repeated the calibration of λ for sources in the RoboPol mask and found a value of 0.01. We use the same set of processing parameters for both polarization standard stars and DP target stars.

The instrument calibration involves using polarization standard stars to (a) determine the polarization zero-point, (b) estimate theuncertainty of this determination (systematic error), and (c) identify the rotation of the coordinate system compared to the standard (IAU) reference frame. In the majority of works using RoboPol data, calibration is done using an instrument model, which is constructed by scanning standard stars across the instrument FOV (King et al. 2014). In this way, the polarization zero point is found for every point on the CCD. This approach, however, only provides an estimate of the systematic uncertainty for sources observed in the field of the instrument (Panopoulou et al. 2015). Our sources were observed in the mask, so as to minimize systematic uncertainties. For this reason, and to eliminate unknown uncertainties in the model determination, we do not make use of an instrument model for the calibration of the DP sources. Instead, we use the measurements of the standards to find the weighted mean instrumental zero-point and its corresponding uncertainty.

The set of polarization standards used for calibration are shown in Table 1, along with their literature values. Though our observations were conducted in the R-band, not all standard stars have reference values in this band. However, most stars which do not have a measurement in the R-band are polarized at a level well below the typical systematic error of our instrument (0.1% in the mask). Assuming that their polarization is interstellar, differences between the V - and R-bands will be negligible for our purposes. The only exception is the star BD+33 2642, which has a reference value of p = 0.2% in V. We present a determination of the R-band value of this standard in Appendix A. We use this value for the following calibration steps.

To find the instrument zero-point, we first subtract the literature value (q_lit) from each measurement of a standard star (q_obs): $\overline{q}=q_{\text{obs}}-q_{\text{lit}}$ and ū =u_obs − u_lit. The zero-point of the instrument q_inst (u_inst) is the weighted mean of all $\overline{q}$ (ū).

The uncertainty of the zero-point reflects the level of systematic error. For its determination, it is of critical importance to take into account all the factors that can contribute to this uncertainty. These include (a) possible intrinsic variability of the standard polarization stars, (b) errors in the determination of the literature value of a certain star (if one makes use of multiple stars for thedetermination), (c) (spatio-temporal) variability of the sky conditions, and (d) (spatio-temporal) variability of the instrument behavior.

Although the zero-point is found using the weighted mean of measurements, the standard error on the mean cannot capture all the aforementioned factors. In the limit of a very large number of measurements of standard stars, the standard error would tend to zero, even though these sources of error would still be at play. In order to properly quantify theaforementioned effects, we assign the uncertainty on the zero-point to be the standard deviation of the $\overline{q}$ measurements (and correspondingly for ū). This is a conservative approach compared to the standard error on the mean and is more likely to err on the side of caution, that is, it is likely to overestimate the systematic uncertainty.

Figure 3 shows the $\overline{q}-\overline{u}$ plane of theobserved standards. Measurements of standards during the observation time span of DP1 are illustrated in the left panel, in the middle panel for DP2, and in the right panel for DP3. It is clear thatfor all observing epochs, the instrument biases the observations towards more positive q and more negative u values. The instrument zero-point is marked with a red star. In Table 2, we present the instrument zero-point (q_inst, u_inst) for all the regions, and its corresponding uncertainty. We find the systematic uncertainty to be at the level of 0.1%.

In addition to the zero-point shift, instrumental effects may also result in a rotation of the q− u plane compared tothe standard (EVPA zero at north, increasing towards east). In practice, to calculate the instrumental rotation we select polarized standards and correct their values for the zero-point shift found previously. Then, we find the average EVPA from these corrected q, u (χ_pol,mean) and subtract from it the literature value of the EVPA: χ_pol,mean − χ_lit.

This rotation is very small for RoboPol, with measurements of polarized standards in 2017 and 2015 placing it at 0.5° ± 0.1°, much smaller than statistical uncertainties for all values of EVPA quoted in this work.

Fig. 3 The $\overline{q}-\overline{u}$ plane of standard polarization stars. The red star illustrated in each case is the weighted mean of measurements and indicates the instrumental zero-point. The labels correspond to the names of the standard stars.

Fig. 4 Measurements of Stokes parameters, q and u, of stars in the DP fields, after instrumental calibration (Sect. 3) shown with blue points. Error bars correspond to the combined statistical and systematic uncertainties. The red star in each panel corresponds to the weighted mean of q and u measurements in each DP (qmean, umean in Table 3). Empty circles denote the outliers defined in Sect. 4.2.

4 Results

4.1 Measurements on the q−u plane

After reducing the data with the RoboPol pipeline, we correct each measurement of a DP target for the instrumental zero-point, and propagate the statistical and systematic uncertainty to the final result. We plot the corrected measurements on the q− u plane in Fig. 4. As is clear from Eqs. (1) and (2), p measures the offset from the origin, while χ measures the (half) angle with respect to the line of positive q, in the counter-clockwise direction.

In each DP, most of the measurements cluster near the origin. Only two stars in DP1, three in DP2, and two in DP3 have signal-to-noise ratio (S/Ns) in p higher than 3. We want to investigate if this is due to the fact that the observed stars are nearby, meaning the polarization does not trace the full extinction across the line-of-sight. We use Gaia Data Release 2 (Gaia Collaboration 2018) to derive the distances of the observed stars. In DP1, there are parallaxes for 22 stars out of 24, in DP2 for 22 out of 23, and in DP3 parallaxes exist for all the stars we observed. However, inverting parallaxes is not a reliable method for inferring distances (Luri et al. 2018). For this reason we use the distances published by Bailer-Jones et al. (2018)³.

In Fig. 5, we present the degree of polarization (p) versus distance for all the stars with associated distances. The majority of measurements are non-detections (S∕N in p less than 3) and we present their upper limits (blackarrows) within 3σ. For stars with S∕N > 3, we present their observed uncertainties. For several stars, the distance uncertainty is negligible and is not visible. We compare this with the line-of-sight distribution of E(B–V) for each DP using the latest version of the three-dimensional (3D) Galactic E(B–V) map⁴ of Green et al. (2018). This map’s beam size ranges from 3.4′ (for high extinction regions) to 13.7′ (for low extinction regions). We query the map for the E(B–V) as a function of distance across the line-of-sight at which each DP is centered. In DP1, DP2, and DP3 the plateau of maximum E(B–V) is reached at distances of 398, 501, and 794 pc, respectively (shown with vertical dotted lines in Fig. 5). Altogether, in all DPs, the majority of stars are far enough to trace the full extinction across the line-of-sight. Therefore, our measurements are tracing the line-of-sight-averaged polarization.

Our approach of observing multiple stars within a small area enables us to infer the mean interstellar polarization towards these regions, even though we do not have a significant detection of p for individual stars. This information is encoded in the observed anisotropy towards a certain direction on the q− u plane. Starting with DP3, the clustering of q, u measurements towards the first quadrant indicates a non-zero mean polarization; the location of this clustering is related to the meandirection of the local polarization of this region. As we move from DP3 to DP1, the anisotropy of measurements around the origin becomes less pronounced. In DP1, most of the measurements are distributed roughly isotropically around zero, appearing consistent with a non-polarized region at the accuracy level of our instrument. DP2 measurements show some clustering towards the first and fourth quadrant, although the q and u measurements are not as anisotropic as in DP3. In Sect. 4.3, we calculate the mean p and χ of each region using the weighted mean q and u (marked with a red star in Fig. 4).

4.2 Search for indicators of intrinsic polarization

While most stars are clustered on the q−u plane, there are some prominent outliers; these are the measurements for which p is further than 3 σ_p from the mean p. More specifically, in DP1 there are two outliers in the first and fourth quadrants of the q− u plane. In DP2 there are three outliers: one in the second quadrant, located opposite the majority of the measurements, one in the fourth quadrant, and the last one in the first quadrant. These measurements are denoted with empty circles in Figs. 4 and 5. Inspecting the latter, we cannot attribute the high p we measured to the fact the stars are far away, except for the outlier in DP1 located at a distance of 2100 pc.

One possible reason for the existence of these outliers could be that they are intrinsically polarized. If this is the case, their measurements should be excluded from our analysis of the properties of the mean interstellar polarization in the surveyed regions. We therefore searched for complementary information on the sources that could help us judge whether they are potentially intrinsically polarized.

One type of source that exhibits intrinsic polarization is an active galactic nucleus (AGN; Angel & Stockman 1980). AGNs are easily distinguished from stars by their non-black-body multiwavelength emission. In our search, we used VOSA⁵ (Bayo et al. 2008), a tool which uses historical multi-band photometric data in order to construct the spectral energy distribution (SED) of a source. We fit the simplest stellar spectral model of VOSA to the outlier sources and find that all are consistent with a black body spectrum. We also computed their effective temperatures, as these can indicate if a star is young (and therefore likely to have a polarization-inducing circumstellar disk). However, all the temperatures found are typical of main sequence stars.

Our second approach to search for intrinsic variability utilizes data from the second data release of the Catalina Sky Surveys⁶ (Drake et al. 2009). For each of our targets, we inspected the seven-year photometric light curves and found no sign of intrinsic (photometric) variability. In addition, we checked the observed B-V colors of the stars in our sample to see if there could be any Be star candidates. We used the SDSS G-R color and converted to B-V according to Jester et al. (2005). We found no stars with negative B-V values (that would be consistent with O or B types). Finally, according to the Besancon stellar population synthesis model of the Galaxy (Czekaj et al. 2014), no stars of type O-A are found within our survey magnitude range in the observed regions. We therefore proceed using the entire sample of observed targets.

Fig. 5 Degree of polarization vs. distance for all the stars with measured parallaxes by Gaia in the DP fields. Downward-pointing arrows represent the 3σ upper limit in polarization fraction when p∕σp < 3. Some distance error bars are too small to be seen. The dotted vertical lines correspond to the distance at which maximum E(B–V) is reached along the line-of-sight of each region according to Green et al. (2018). Empty circles denote the outliers defined in Sect. 4.2.

4.3 Mean polarization

In order to measure the mean fractional linear polarization, p_mean, and the mean EVPA, χ_mean, for each DP, we computed the weighted mean of the q and u measurements of each region. We obtain two single values (q_mean and u_mean) for each region, and we apply Eqs. (1) and (2) to derive p_mean and χ_mean. We follow this approach for two reasons. First, it minimizes the contribution to the mean of the aforementioned outliers which might not be consistent with ISM polarization. Second, this approach avoids the bias of individual measurements. The majority of our measurements have low p∕σ_p and are therefore biased towards higher values of p (e.g., Simmons & Stewart 1985; Vaillancourt 2006). Therefore, if we were to compute the p_mean of each region by averaging individual stellar p, the p_mean value would be overestimated.

In Table 3, we present the mean values of the Stokes parameters and p and χ in each region. The p_mean in DP1 lies within 1.4σ of the origin, while that of DP2 is 3.1σ from zero, and that of DP3 is 4.7σ from zero. Therefore, we have measured significant polarization in DP3 and DP2 but not in DP1. We consider the p_mean of DP2 to be alower bound on the detectable mean polarization that can be achieved using our methodology.

The isotropy of q, u measurements of DP1 about the origin is consistent with the non-detection of polarization. As we move from DP2 to DP3, the anisotropy of q, u measurements becomes more prominent. We elaborate more on the anisotropy of the measurements with regard to the significance of the p_mean in the DPs in Sect. 4.5.1.

4.4 Intrinsic spread in the distribution of Stokes parameters

Due to the distribution of stars in three dimensions, different stars could be tracing different materials along the line-of-sight. As a result, deviations from the mean polarization (or equivalently from q_mean and u_mean) can arise. This means that there is an underlying intrinsic distribution of q and u, respectively, which needs to be fully characterized. Except from q_mean and u_mean, our strategy of measuring a large number of stars within a small region of the sky allows us to constrain the intrinsic spread of the distributions of q and u. We proceed to obtain such constraints only for the two regions where we have detected significant p_mean (DP2, DP3).

The spreads (sample standard deviations) of the observed q and u are s_q,obs = 0.31% and s_u,obs = 0.13% for DP2, and s_q,obs = 0.22% and s_u,obs = 0.12% for DP3. The errors in each individual stellar measurement are comparable to the observed standard deviations of the distributions. This means that the spreads of these distributions are determined by the uncertainty of our measurements and not by a physical process.

We use a Bayesian approach to estimate the intrinsic spread of the distributions. First we make the assumption that the intrinsic q, u of the stars in a region are normally distributed about q_mean, u_mean with equal spreads s_q,0 = s_u,0 = s₀. We also use the fact that measurements of the Stokes parameters of an individual star have Gaussian uncertainties. The likelihood of observing N stars with q_obs,i, u_obs,i (i =1, ..., N) and measurement uncertainties σ_q,obs,i and σ_u,obs,i if the intrinsic spread is s₀ is $$\begin{array}{lll}L\hfill & =\hfill & ({\displaystyle \prod _{i=1}^N\frac{1}{2\pi \sqrt{(s_0^2+{\sigma }_{q,\text{obs},i}^2)(s_0^2+{\sigma }_{u,\text{obs},i}^2)}}})\hfill \\ \hfill & \hfill & \times \text{exp}[-\frac{1}{2}({\sum }_{i=1}^N\frac{{(q_{\text{obs},\text{i}}-q_{\text{mean}})}^2}{{\sigma }_{q,\text{obs},i^2+s_0^2}+\frac{{(u_{\text{obs},i}-u_{\text{mean}})}^2}{{\sigma }_{u,\text{obs},i}^2+s_0^2}})].\hfill \end{array}$$

The analytical proof of the likelihood function can be found in the appendix of Venters & Pavlidou (2007).

To obtain an estimate of the true s₀ we search for the value of s₀ that maximizes the likelihood function in Eq. (4). Figure 6 shows the likelihood L as a function of s₀ for DP2 (left) and DP3 (right). The functions are normalized so that the area under each curve is equal to 1. A vertical dotted line shows the value of s₀ which maximizes the likelihood. In the case of DP2, this corresponds to s₀ = 0.07%, while for DP3 the likelihood peaks at s₀ = 0%. Both likelihood functions in Fig. 6 are bounded by zero and it is therefore not possible to obtain symmetric bounded confidence intervals on s₀. We can, however, place upper limits on s₀ for both regions: the 99% upper confidence interval is s₀ ≤ 0.187% for DP2 and s₀ ≤ 0.127% for DP3.

For the purpose of this work, deriving upper limits on s₀ is sufficient. In DP3, we find the maximum likelihood s₀ to be zero. This is reasonable since the majority of measurements are within ~ 1σ of the q_mean and u_mean. The situation in DP2 is altered mainly because of the presence of the outliers.

Fig. 6 Normalized likelihood L as a functionof the intrinsic spread s0 of the distributions of Stokes parameters for DP2 (left panel) and DP3 (right panel). A vertical line shows the maximum likelihood s0. The 99% confidence intervals for DP2 and DP3 are s0 ≤ 0.187% and s0 ≤ 0.127%, respectively.

4.5 Confidence of the measured mean fractional linear polarizations

Due to the low signal in the DPs, we must assess the confidence that we can place on the presented mean interstellar polarization towards the three regions. We do this by quantifying the likelihood of a false detection, that is, the likelihood that the signal could result simply from uncertainties in our analysis. We use two observables to test this null hypothesis: the anisotropy of measurements on the q−u plane (Sect. 4.5.1) and the ratio of the mean fractional linear polarization over its uncertainty, p_mean∕σ_p,mean (Sect. 4.5.2).

4.5.1 Significance of q−u plane anisotropy

In DP1, there are seven measurements located in the first quadrant, eight in the second, six in the third, and three in the fourth. In DP2, there are nine measurements in the first, seven in the second, zero in the third, and seven in the fourth quadrant. For DP3, the numbers are ten for the first, three for the second, zero for the third, and eight for the fourth quadrant. To quantify how anisotropically the measurements are distributed around zero, we define an isotropy parameter κ as the ratio of the number of q, u measurements detected in the quadrant with the least measurements over the number of points detected in the quadrant with the most measurements. If a quadrant is empty, κ = 0, while if all the measurements are isotropically distributed, κ = 1. For DP3 and DP2, κ = 0 because the third quadrant in each region has no measurements. In DP1, the fourth quadrant has the least number of points, three, while the second has the most, eight, and therefore κ = 0.375.

In order to investigate how probable it is to reproduce the observed anisotropies from unpolarized stars, we performed Monte Carlo simulations to obtain the probability distribution of κ under the null hypothesis that all stars are unpolarized, (q_true, u_true) = (0,0). The q_i and u_i measurements follow Gaussian distributions centered on q_true,i and $u_{\text{true},i}$ with standard deviations σ_q,i and σ_u,i, respectively, with i = 1,…, N where N is the number of stars. Assuming all stars are unpolarized (q_true,i = 0 and u_true,i = 0), we created mock observations by drawing random values for q, u from zero-centered Gaussians with standard deviations equal to the observational uncertainties of each star. We produced a sample of N mock (q_mock,i, u_mock,i) sets that match the number of stars measured in each region, and we computed the parameter κ. We repeated this process 10⁶ times, produced the distribution of mock κ, and then compared it with the κ value obtained from the real data. For DP1, we find that 40% of the mock κ are smaller than the observed one. For DP2 and DP3, this probability is 0.96% and 0.95%, respectively.

We conclude that the anisotropy of observations on the q− u plane in DP1 is consistent with the anisotropy produced by stars with zero ISM polarization (unpolarized). In the other DPs, the observed anisotropy could not be produced by unpolarized stars, with confidence more than 99%.

4.5.2 Significance of p_mean∕σ_p,mean

In this section, we calculate the probability of measuring the observed weighted mean p S/N (p_mean∕σ_p,mean), if all stars in each DP field were unpolarized, given the uncertainties in our analysis. These uncertainties include the photon noise of individual DP target star measurements, but also uncertainties in the measurements of the standard stars that are used for our calibration (zero-point offset correction). We perform the calculation as follows.

For each standard star used in the calibration of a DP, we draw a mock observation from a Gaussian centered on the existing measurement with a standard deviation equal to the statistical error of the measurement. We then compute the weighted mean of this mock set of standard observations (as in Sect. 3). This is the zero-point, around which any unpolarized star should lie. We now generate mock observations of the stars in the DP, assuming they are unpolarized. In practice, for each star we draw a value (q and u) from a Gaussian centered on the zero-point that we had just calculated, with a standard deviation equal to the observed photon-noise error of the star. We now have mock observations of hypothetical zero-polarized stars in the DP. Next, we correct each mock star measurement for the instrumental polarization using the generated zero-offset and standard deviation, as in Sect. 3. We then calculate the weighted mean and error of these corrected mock measurements. This process is repeated 10⁴ times.

We construct the distribution of p_mean∕σ_p,mean from this test for each DP in Fig. 7. The distributions for all three DPs are very similar. There is a peak at 0.8 and a long tail that extends out to 3.1 (3.6 for DP3). The observed p_mean∕σ_p,mean in the DPs (from Table 3) are shown with vertical lines. The p_mean∕σ_p,mean of DP1 falls well within the spread of the distribution, showing that it can be completely explained by the uncertainties present in our analysis. That of DP2 falls at 3.13, slightly higher than the maximum value of the 10⁴ DP2 zero-polarization realizations. The p_mean S/N of DP3 is much larger than the corresponding maximum value of 10⁴ DP3 zero-polarization realizations.

Both tests in this section are in agreement that DP2 and DP3 have yielded significant detections of the mean polarization, in contrast to DP1.

4.6 Characteristics of HI emission towards the DPs

Having assessed the confidence of our measurements of the mean polarization in the three regions, we can now study the source of the signal; that is, the diffuse atomic medium, as traced by HI line emission.

We use spectral data cubes from the first data release of the Effelsberg-Bonn HI survey (Winkel et al. 2016), which have a beam size of 10′. For each DP, we locate the pixel in the HI data that corresponds to the center of the observed field. We then average the HI spectra of this pixel and its eight nearest neighbors (which yields an averaging area ~ 10 ′ in width). Figure 8 shows the (averaged) HI spectra for all DPs (blue for DP1, orange for DP2, green for DP3). A vertical dotted line shows the velocity, v_max, at which each spectrum intensity is maximum. For DP1 and DP2, v_max is –52 and –9 km s⁻¹, respectively.For DP3, we show both the primary and the secondary peaks which correspond to velocities of –13 and − 1 km s⁻¹, respectively.

Inspecting the spectrum of DP3, we distinguish two prominent peaks at velocities that imply that the gas is local. These two intense components likely account for the main contribution to the polarization we observe. In the spectrum of DP2, there is a single component that dominates the signal, and it is also local. The marginally detected polarization of DP2 can be attributed mostly to this component. Meanwhile, the non-detection of polarization in DP1 could be the result of the much lower total intensity of the two components seen in the spectrum. Another alternative is that the component at high (absolute) velocities, which would be the main contributor to the column, is at a distance such that the bulk of the stars in our sample are foreground to it. Finally, there is also the possibility that the low polarization is a result of depolarization occurring along the line-of-sight to DP1, due to the presence of these two components.

By integrating over (approximately) the width of each velocity component, we can inspect the local morphology of the HI emission on the plane of the sky. Recent works have found the structure of the HI gas to be well correlated with the orientation of the magnetic field (using starlight polarization and HI data, Clark et al. 2014; or dust thermal emission and its polarization, Planck Collaboration Int. XXXII 2016). We investigate the relation of HI structure to the mean polarization in our regions.

Figures 9 and 10 show maps of the HI intensity integrated over the velocity range [v −Δv, v + Δv], where Δv has been set equal to 6 km s⁻¹ and v is the velocity at which the spectra peak (including the secondary peak of DP3 at –1km s⁻¹). The maps are centered on DP2 (Fig. 9) and DP3 (Fig. 10). DP1 is not presented because its p_mean is consistent with zero. For each region, the black (or white) solid line segment forms an angle χ_mean with respect to the north (increasing towards the east according to the IAU convention for the EVPA) and shows the mean polarization orientation. Its length is proportional to p_mean (Table 3). The dotted lines around the solid segment indicate the error in the mean polarization angle and the green boxes mark the surveyed regions. A polarization segment indicating the polarization scale is plotted in the top left of each image.

Inspecting Fig. 9, we find that the single component of the spectrum in DP2 shows a filamentary morphology.The orientation of the observed structure is not aligned with that of the mean polarization, at the resolution of 10′. The HI emission around DP3 is much more complex. The integrated intensity of the velocity component seen at –13 km s⁻¹ is shown in the top panel of Fig. 10. DP3 appears to fall on the edge of a low-aspect-ratio (“blobby”) structure, whose mean orientation is not clear. There does seem to be an asymmetry of the emission indicating an orientation from the south-east towards the north-west. The mean EVPA does not coincide with the axis of the asymmetry. The second velocity component, centered around –1 km s⁻¹, does not allow for a clear determination of a mean orientation in the HI morphology (bottom panel, Fig. 10). At higher resolution, such a comparison may be facilitated; for example, if the high-aspect-ratio structures identified by Clark et al. (2014) in GALFA-HI data (4′) are present in this region. Finally, it is possible that the polarization of stars in DP3 is affected by both velocity components, meaning that an alignment of the polarization in each component with the corresponding HI structures cannot be excluded.

Fig. 7 Distribution of the ratio of weighted mean p over its error for random realizations of hypothetically unpolarized DP stars. The vertical lines show the actual ratio for our DP measurements. Colors: gray – DP1, red – DP2, blue – DP3.

Fig. 8 HI line intensity as a function of radial velocity for regions of 10′ in size, centered on the DPs. Vertical dotted lines correspond to the velocities at which the intensity peaks.

Fig. 9 Mean polarization orientation over-plotted on an image of HI integrated intensity centered on DP2. Velocities are integrated around –9 km s−1. The green box marks the surveyed region. The segment in the top-left corner is for scale in p.

Fig. 10 As inFig. 9, but for the two components in the spectrum of DP3. Top panel: velocity range centered on –13 km s−1. Bottom panel: velocity range centered on –1 km s−1.

4.7 Comparison of p_mean with E(B–V)

Measurements of interstellar p are bound by an upper envelope in the p −E(B−V) plane. The first work to define such an envelope empirically was that of Hiltner (1956), using observations of 1259 O and B stars with E(B−V) ≥ 0.13 mag (we convert his presented total extinctions (A_V) to reddenings assuming a ratio of total-to-selective extinction R_V = A_V∕E(B−V) = 3.1). He found that the majority of measurements are bound by the relation $$p_{\text{max}}\approx 9E(\text{B-V})(\%/\text{mag}),$$(5)

where we have used R_V = 3.1 and the conversion from p in magnitudes to the fractional polarization (Whittet 1992). Following works confirmed the existence of this envelope (Serkowski et al. 1975, E(B−V) > 0.1 mag). In a much later work, Fosalba et al. (2002) used stars with 0.01 mag < E(B-V) < 1 mag to fit an expression for the mean p as a function of E(B–V): $$\langle p\rangle =3.5*E{(\text{B-V})}^{0.8}(\%/\text{mag}).$$(6)

We note that $\langle p\rangle$ is the mean p of a sample of stars within a given range in E(B–V), and is therefore distinct from our determination of the mean fractional linear polarization within the DPs (p_mean).

We wishto compare the p_mean found in Sect. 4 to relations (5) and (6) above. In order to make a fair comparison, we must take into account that the observed p_mean is a biased estimator of the true p_mean (e.g., Simmons & Stewart 1985). We calculatethe debiased p_mean using the formula of Vaillancourt (2006): ${\widehat{p}}_{\text{mean}}=\sqrt{p_{\text{mean}}^2-{\sigma }_{p,\text{mean}}^2}$, for DP2 and DP3 where p_mean > 3σ_p,mean. The debiased estimate of p_mean in DP1 is zero.

A directestimate of each star’s reddening cannot be obtained with the available information on our sample. As an estimate of the mean reddening, we use the total E(B–V) from the LHD map towards each DP. This most likely overestimates the reddening for the starsthat are not tracing the full line-of-sight. As there are systematic uncertainties associated with the conversion from HI column density to E(B–V), we adopt the following approach to estimate this uncertainty for the DPs. Schlegel et al. (1998) derived LHD E(B–V) from the relation of dust-emission-based E(B–V) to HI column, NHI (Fig. 1, right, in LHD). For a given NHI, there is a range of E(B–V) values observed, likely resulting from the combination of different effects (including variations in the gas-to-dust ratio and effects of cosmic infrared background fluctuations etc.). For each DP, we have found the NHI and the spread of the corresponding E(B–V) ± 0.001 values(from the entire LHD sky footprint). We assign this spread as the error to the LHD E(B–V) in each DP.

In Fig. 11, we compare our measurements to Eq. (5) and (6). The figure shows both p_mean and the debiased ${\widehat{p}}_{\text{mean}}$ as open andfilled circles, respectively. We find that both of our significant detections lie higher than both the $\langle p\rangle$ and p_max curves. The ${\widehat{p}}_{\text{mean}}$ of DP2 lies 1.1 σ from both relations. The ${\widehat{p}}_{\text{mean}}$ of DP3 lies 2.2 σ from Eq. (5) and 2.3 σ from Eq. (6). Since our estimate of E(B–V) is an upper limit, it is likely that the points will be shifted towards the left, and the inconsistency with the p_max envelope will be augmented.

Polarization measurements at higher extinctions have consistently shown agreement with Eq. (5). Throughout the literature, only a small fraction of stellar p measurements lie above the relation (e.g., Fig. 1(a) from Andersson et al. 2015 and Fig. 15 in Panopoulou et al. 2015). At E(B–V) > 0.03 mag, the measurements of Santos et al. (2011) are also largely consistent with this upper envelope. This is the first comparison with the aforementioned relations for such low extinctions.

Interestingly, the relation for the upper envelope (5) coincides with that for $\langle p\rangle$ for E(B–V) < 0.01 mag. This implies that at least one of the two relations will lead to an erroneous estimate of p for a given E(B–V). Both relations have been calculated at higher extinctions and have been extrapolated to these low E(B–V). However, in the case of the p_max curve, the extrapolation has been made from E(B–V) that are an order of magnitude higher than the ones studied here. The Fosalba et al. (2002) dataset used points down to E(B–V) ~ 0.01 mag and is therefore more reliable for these low extinctions.

This result implies that p_max has previously been underestimated by an unknown factor, which may have implications on existing dust models (e.g., Planck Collaboration Int. XXI 2015). Even at extinctions of 0.1 mag, significant depolarization exists along the line-of-sight.If we consider a model of the 3D dust distribution as discrete polarizing “screens” (or clouds), our findings would suggest thatthe number of “screens” along the line-of-sight is much smaller than that at higher extinctions. In fact, the HI emission spectra of the DP fields in Sect. 4.6 would suggest that 1–2 distinct components exist at these sightlines. We reserve a more detailed investigation of the polarizing efficiency of the low-dust-extinction sky for a future work.

Fig. 11 Fractional linear polarization vs. reddening. Open circles show the pmean of the DPs, while filled circles show the debiased ${\widehat{p}}_{\text{mean}}$ (Table 3). The E(B–V) is found from the LHD map as explained in the text. The solid black line indicates the upper envelope from Eq. (5) (Hiltner 1956). The expression for the mean p, Eq. (6) (Fosalba et al. 2002), is traced by the dashed gray line. The gray regions indicate the interval in which E(B–V) values lie within 99% confidence.

5 Discussion

We set out to perform our stellar polarization mini-surveys in order to identify requirements for future optopolarimetric experiments targeting the high-latitude sky. With our setup of ~20 stars in ~ 0.05 square degrees, and a systematic uncertainty of 0.1% in p, we find the vast majority of individual stellar measurements to be non-detections. However, our strategy enables us to measure the mean fractional linear polarization within our target regions, with good enough precision so as to obtain significant detections in two out of three regions. With our systematic uncertainty and survey depth R ≤ 16 mag, we obtain a marginal detection of p_mean in DP2 (~0.11% ± 0.04%, Table 3, or p_mean∕σ_p,mean ~ 3). The level of the signal is low, but recoverable.

Our survey strategy is qualitatively different from that of existing works at high latitudes, which have selected bright stars sampling sightlines that are distant on the plane of the sky. On the other hand, we have targeted a flux-limited sample in each of three very small areas of the sky. It is of interest to compare the results of the two different approaches.

The vast majority of existing high-latitude measurements belong to the surveys of (Berdyugin et al. 2001, 2014), and Berdyugin & Teerikorpi (2001, 2002). These works have cataloged the polarization of over 2800 bright stars (V ≤ 13 mag) with known distances (from the HIPPARCOS mission). Their sample is located at Galactic latitudes b > 30°, b < −60°. The bulk of stellar p from the aforementioned catalogs are clustered around 0.1%. We find that the p_mean measured in this work are within the range of observed polarizations of the Berdyugin sample.

6 Conclusions

In this work we performed stellar optical polarization surveys in three ~ 15′ × 15′ regions of the high-galactic latitude sky. Our aim was to determine the level of interstellar fractional linear polarization (p) that can be recovered in regions of very low dust emission using the RoboPol polarimeter. Our surveys are photometrically complete down to R ~ 16 mag, and are therefore unique in depth for the high-Galactic-latitude sky. Our findings can be summarized as follows.

The determination of the systematic uncertainty is critical for our goal of measuring interstellar p in regions of very low dust emission. We have taken care to provide a reliable estimate of the uncertainty of our instrument for every observing run. We find the uncertainty to be at the level of 0.1% in p (in the focal plane mask where our observations were conducted; Sect. 3). Furthermore, we have provided a measurement of the intrinsic polarization of the standard star BD+33 2642 in the R-band (Appendix A).

We havedetected significant p for only seven stars in our samples. Most of them are outliers in the Stokes q and u plane when compared to the location of the bulk of the measurements, but we could not identify indications of variability for these (or any other) stars in our sample (Sect. 4). Even though most measurements have yielded only upper limits on p, we are able to locate the mean interstellar p with high significance in the region with highest dust content (reddening), DP3, at the level of p_mean = (0.208 ± 0.044)%. In the region with intermediate reddening, DP2, we have obtained a marginal detection at p_mean = (0.113 ± 0.036)%. In the region with least dust content, DP1, the p_mean = (0.054 ± 0.038)%, has not been significantly detected (Sect. 4). In addition, we have estimated the intrinsic spread of the Stokes q and u distributions, under the assumption that the distributions are Gaussian and have equal spread. We place upper limits on the intrinsic spread of 0.187(%) for DP2 and 0.127(%) for DP3 (Sect. 4). Through statistical tests we have assessed the confidence of the measured p_mean and have found that the signal in DP2 and DP3 cannot be produced by uncertainties in our analysis (Sect. 4.5). Using HI line emission from the EBHIS survey we identify two dominant components in DP3, one in DP2 and two much fainter components in DP1 (Sect. 4.6). The morphology of the component in DP2 is filamentary, and is not aligned with the mean EVPA at the resolution of 10′. In DP3 the morphology of both components is more complex and cannot be easily compared to the mean polarization angle (Sect. 4.6).

Acknowledgements

We would like to thank the reviewer for insightful comments. In addition, we thank K. Kovlakas, T. Kougentakis, N. Mandarakas, G. Paterakis, and A. Steiakaki for technical support and as well as N. D. Kylafis, B. Hensley, I. K. Wehus and D. Lenz for their helpful comments on the paper. RS would like to thank Dr. K. Christidis for fruitful discussions. GVP acknowledges support from the National Science Foundation, under grant No. AST-1611547. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 771282. This publication makes use of VOSA, developed under the Spanish Virtual Observatory project supported from the Spanish MICINN through grant No. AyA2011-24052. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC; https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Appendix A Reference value of the polarization standard star BD+33 2642 in the R-band

During the course of the 2017 RoboPol observing season, a large sample of polarization standard stars was observed with high cadence. The wealth of data, combined with the stability of the instrument, enables us to characterize the properties of individual standard stars. We select the most well-sampled and stable stars from our set of calibrators (to minimize any possibility of intrinsic variability) and use them to derive the Stokes parameters of BD+33 2642 in the R-band. Apart from the standard stars in Table 1, we also include the star HD154892 (p = 0.05 ± 0.03, in the B-band; Turnshek et al. 1990)in this analysis.

Figure A.1 shows the measurements of $\overline{q},\overline{u}$ (literature-subtracted measurements of standards, following Sect. 3) for this set of stars observed throughout the 2017 season (May–November). All observations were conducted in the R-band and reduced as in Sect. 3. For BD+33 2642, we use its literature value for the V -band (Table 1). It is clear that this star exhibits an offset from the bulk of standards in the $\overline{q}$,ū plane. Calculating the weighted mean of the $\overline{q}$ of BD+33 2642 (marked with a purple star in Fig. A.1), we find that it lies ~ 3σ away from the weighted mean $\overline{q}$ of the other star measurements (blue star). We show the weighted mean from Table 2 for DP2 for comparison (black ×). This was found using only standards from the nights when DP2 targets were observed (black ×). The two zero-point estimates are within 1σ of each other, demonstrating the stability of the instrument.

The only plausible explanation for the observed offset of BD+33 2642 is that its polarization properties in the R-band are different than those in the V -band. We calculate its “true” Stokes parameters in the R-band in the following way. We find the weighted mean $\langle \overline{q}\rangle ,\langle \overline{u}\rangle$ of the measurements of all other standard stars. We also find the weighted mean of the measurements of BD+33 2642 (not corrected for the literature value) $\langle q_{B33}\rangle ,\langle u_{B33}\rangle$. An estimate of the “true” Stokes parameters is $$q_{\text{true}}=\langle q_{B33}\rangle -\langle \overline{q}\rangle ,u_{\text{true}}=\langle u_{B33}\rangle -\langle \overline{u}\rangle ,$$(A.1)

and their uncertainty: $${\sigma }_{q,\text{true}}=\sqrt{{\sigma }_{\langle q\rangle ,B33}^2+{\sigma }_{\langle \overline{q}\rangle }^2},{\sigma }_{u,\text{true}}=\sqrt{{\sigma }_{\langle u\rangle ,B33}^2+{\sigma }_{\langle \overline{u}\rangle }^2},$$(A.2)

where ${\sigma }_{\langle q\rangle ,B33}$ and ${\sigma }_{\langle \overline{q}\rangle }$ (and similarly for u) are the weighted standard deviations of the valuesused for calculating $\langle q_{B33}\rangle$ and $\langle \overline{q}\rangle$, respectively.As in Sect. 3, we use the weighted standard deviation of measurements as an estimate for the error on the weighted mean, as a conservative choice. In this way, we take into account any possible intrinsic variability of the sources.

The resulting values for the polarization of BD+33 2642 are shown in Table A.1. The star has been observed in previous years less times than in 2017. We have not used measurements from previous year for this determination. However, we have checked that repeating the same process with the data from seasons 2015 and 2016 does not produce “true” polarization parameters that are inconsistent with the determination presented here.

Fig. A.1 Measurements of the residual q,u of standard stars after subtraction of their literature values (Table 1 circles: each color represents a different star). Blue star: weighted mean of $\overline{q},\overline{u}$ excluding BD+33 2642. Purple star: weighted mean of BD+33 2642 measurements corrected using the literature value in the V -band (Table A.1). The black “×” marks the weighted mean found for DP2 in Sect. 3.

We note that in the analysis of Sect. 3, the R-band polarization of BD+33 2642 is used only for the calibration of DP1, which was observed in 2015. Therefore, we are using completely independent measurements to derive its “true” polarization.

As a final remark, we find the polarization angle in the R-band to be significantly different from that in the V -band. This may indicate that the polarization properties of BD+33 2642 have changed since the original measurements of Schmidt et al. (1992). This calls for verification with the use of different instruments.

References

All Tables

All Figures

	Fig. 1 Mollweide projection of the E(B–V) map of Lenz et al. (2017; centered at l, b = [0,0]; grid-line spacing is 30°). The red stars mark our target regions. Gray areas are not included in the map.
In the text

	Fig. 2 Distribution of E(B–V) from the map of Lenz et al. (2017). Vertical lines correspond to the E(B–V) of the target regions DP1, DP2, and DP3. The truncation at 0.045 mag corresponds to the mask placed by Lenz et al. (2017) at hydrogen column densities NHI > 4 × 1020 cm−2.
In the text

	Fig. 3 The $\overline{q}-\overline{u}$ plane of standard polarization stars. The red star illustrated in each case is the weighted mean of measurements and indicates the instrumental zero-point. The labels correspond to the names of the standard stars.
In the text

Fig. 4 Measurements of Stokes parameters, q and u, of stars in the DP fields, after instrumental calibration (Sect. 3) shown with blue points. Error bars correspond to the combined statistical and systematic uncertainties. The red star in each panel corresponds to the weighted mean of q and u measurements in each DP (qmean, umean in Table 3). Empty circles denote the outliers defined in Sect. 4.2.

In the text

Fig. 5 Degree of polarization vs. distance for all the stars with measured parallaxes by Gaia in the DP fields. Downward-pointing arrows represent the 3σ upper limit in polarization fraction when p∕σp < 3. Some distance error bars are too small to be seen. The dotted vertical lines correspond to the distance at which maximum E(B–V) is reached along the line-of-sight of each region according to Green et al. (2018). Empty circles denote the outliers defined in Sect. 4.2.

In the text

	Fig. 6 Normalized likelihood L as a functionof the intrinsic spread s0 of the distributions of Stokes parameters for DP2 (left panel) and DP3 (right panel). A vertical line shows the maximum likelihood s0. The 99% confidence intervals for DP2 and DP3 are s0 ≤ 0.187% and s0 ≤ 0.127%, respectively.
In the text

	Fig. 7 Distribution of the ratio of weighted mean p over its error for random realizations of hypothetically unpolarized DP stars. The vertical lines show the actual ratio for our DP measurements. Colors: gray – DP1, red – DP2, blue – DP3.
In the text

	Fig. 8 HI line intensity as a function of radial velocity for regions of 10′ in size, centered on the DPs. Vertical dotted lines correspond to the velocities at which the intensity peaks.
In the text

	Fig. 9 Mean polarization orientation over-plotted on an image of HI integrated intensity centered on DP2. Velocities are integrated around –9 km s−1. The green box marks the surveyed region. The segment in the top-left corner is for scale in p.
In the text

	Fig. 10 As inFig. 9, but for the two components in the spectrum of DP3. Top panel: velocity range centered on –13 km s−1. Bottom panel: velocity range centered on –1 km s−1.
In the text

Fig. 11 Fractional linear polarization vs. reddening. Open circles show the pmean of the DPs, while filled circles show the debiased ${\widehat{p}}_{\text{mean}}$ (Table 3). The E(B–V) is found from the LHD map as explained in the text. The solid black line indicates the upper envelope from Eq. (5) (Hiltner 1956). The expression for the mean p, Eq. (6) (Fosalba et al. 2002), is traced by the dashed gray line. The gray regions indicate the interval in which E(B–V) values lie within 99% confidence.

In the text

Fig. A.1 Measurements of the residual q,u of standard stars after subtraction of their literature values (Table 1 circles: each color represents a different star). Blue star: weighted mean of $\overline{q},\overline{u}$ excluding BD+33 2642. Purple star: weighted mean of BD+33 2642 measurements corrected using the literature value in the V -band (Table A.1). The black “×” marks the weighted mean found for DP2 in Sect. 3.

In the text