© ESO 2021

1 Introduction

Cosmic microwave background (CMB) polarization experiments have reached sufficient sensitivity to demonstrate that, even in the most diffuse regions of the sky, cosmological signals of interest lie below the polarized emission from Galactic foregrounds (BICEP2 Collaboration & Keck Array Collaboration 2018; Planck Collaboration IV 2020). In particular, the B-mode signature from primordial gravitational waves (Kamionkowski & Kovetz 2016), quantified by the tensor-to-scalar ratio r, is now constrained to be at least approximately ten times fainter than B-mode emission from Galactic dust at 150 GHz, even in the diffuse BICEP/Keck region (BICEP2 Collaboration & Keck Array Collaboration 2018). Next-generation experiments like the Simons Observatory (Ade et al. 2019), CMB-S4 (Abazajian et al. 2016), and LiteBIRD (Suzuki et al. 2018) seek constraints on r that improve on current upper limits by an order of magnitude or more and will therefore require foreground mitigation at the level of one percent or better.

One of the most challenging aspect of modeling dust foregrounds is that the spectral energy distribution (SED) of dust emission is not uniform across the sky. Variations in dust temperature and opacity law are now well attested across the Galaxy (e.g., Finkbeiner et al. 1999; Planck Collaboration XI 2014; Meisner & Finkbeiner 2015; Planck Collaboration IV 2020; Irfan et al. 2019), with evidence for correlations with gas velocity (Planck Collaboration XXIV 2011; Planck Collaboration XI 2014), strength of the ambient radiation field (Planck Collaboration XXIX 2016; Fanciullo et al. 2015), and location in the Galactic disk (Schlafly et al. 2016).

Such variations greatly restrict the ability to use maps of dust emission at one frequency to constrain dust emission at another frequency; that is, two maps at different frequencies differ by more than just an overall multiplicative factor (frequency decorrelation). The three-dimensional (3D) structure of the interstellar medium (ISM) adds to the complexity of this problem (Tassis & Pavlidou 2015). If a single line of sight (LOS) intercepts multiple dust clouds with different SEDs and magnetic field orientations, then the frequency scaling of each of the Stokes Q and U parameters may be different even in a single pixel (LOS frequency decorrelation). Frequency decorrelation has already been identified as a critical uncertainty in current r constraints and will be even more important at higher sensitivities (BICEP2 Collaboration & Keck Array Collaboration 2018; CMB-S4 Collaboration 2020).

Frequency decorrelation is often quantified at the power spectrum level through the ratio of the BB cross-spectrum of two frequencies at some multipole ℓ to the geometric mean of their auto-spectra (Planck Collaboration L 2017). Computing over large areas of the Planck polarization maps at 353 and 217 GHz, the channels with the greatest sensitivity to polarized dust emission, has yielded limits of only (Sheehy & Slosar 2018; Planck Collaboration XI 2020). While this limit suggests frequency decorrelation may not be a limiting concern if r ≳ 0.01, Planck Collaboration XI (2020) caution that the level of decorrelation may be variable across the sky with some limited sky regions potentially having much greater values.

Line-of-sight frequency decorrelation can have a particularly pernicious effect on parametric component-separation methods working at the map level, especially if the SEDs of Stokes Q and U are not modeled with independent parameters (Poh & Dodelson 2017; Ghosh et al. 2017; Puglisi et al. 2017; Hensley & Bull 2018; Martínez-Solaeche et al. 2018; CMB-S4 Collaboration 2020). New techniques employing moment decomposition (Chluba et al. 2017) have shown promise for mitigating LOS averaging of dust SEDs in polarization at the expense of additional parameters (Mangilli et al. 2019; Remazeilles et al. 2020). Distortions of the SED from effects like LOS frequency decorrelation are also important for power-spectrum-based modeling of foregrounds. In particular, Mangilli et al. (2019) showed that ignoring effects like LOS frequency decorrelation can bias r determinations at consequential levels for next-generation experiments even if a frequency decorrelation parameter is used when fitting an ensemble of power spectra.

In this work, we focus on LOS frequency decorrelation, adopting a different approach, based on the fact that regions of the sky where the effect is expected to be important can be astrophysically identified using ancillary ISM data. Specifically, we use HI emission data to identify sightlines that are potentially most susceptible to this effect. We combine information on the discrete number of HI clouds on each sightline (Panopoulou & Lenz 2020) with an estimate of the magnetic field orientation in each cloud inferred from the morphology of linear HI structures (Clark & Hensley 2019). This entirely HI-based sample selection is “agnostic” to the Planck dust polarization data. We then compare the difference in polarization angles at 353 and 217 GHz along sightlines with and without an expected LOS frequency decorrelation effect, finding that the HI data indeed identify sightlines with more significant electric vector position angle (EVPA) rotation. This is the first detection of LOS frequency decorrelation with Planck data and illustrates the power of ancillary data such as HI and stellar polarizations in identifying regions of the sky where the effect is most pronounced.

This paper is organized as follows. In Sect. 2 we briefly review the phenomenology of frequency decorrelation of polarization. In Sect. 3 we describe the data sets that are used in the analysis. Section 4 presents the sample selection, the statistical tools that are used, and our handling of biases and systematic errors. Section 5 presents our results. We discuss the robustness of our findings, and present further supporting observational evidence in Sect. 6. An estimate of the required SED variation to reproduce the observed magnitude of LOS frequency decorrelation is presented in Sect. 7. We discuss our findings in Sect. 8 and conclude in Sect. 9.

This paper demonstrates that the effect of LOS frequency decorrelation exists at the pixel level and can be measured in the high-frequency polarization data from Planck. It does not address whether the amplitude of the effect is large enough at the sky-map level to affect any particular experimental search for primordial B-modes.

2 Phenomenology of LOS frequency decorrelation

We seek to detect LOS frequency decorrelation between Planck polarization data at 353 and 217 GHz, frequencies dominated by Galactic thermal dust emission and the CMB. Given that the polarized intensity of the CMB and of thermal dust emission feature different SEDs and that they are uncorrelated, their relative contribution to the observed polarization signal depends on the frequency. A change with frequency of the polarization position angle is therefore expected even if the polarization pattern of emission from dust remains constant across frequencies. Additionally, statistical and systematic errors induce scatter in polarization position angles at each frequency. Therefore, a measured difference in polarization direction (EVPA) between frequencies cannot be immediately attributed to a LOS frequency decorrelation induced by multiple dust polarized-emission components. Similarly, when the EVPA difference between frequencies is computed for a large statistical sample of different LOSs, EVPA differences form a distribution with a finite spread. The three sources of EVPA differences mentioned above (noise, relative contributions of the CMB and the dust, and SED difference between dust components) each contribute to the width of the EVPA difference distribution. We wish to detect a signal that can be directly attributed to frequency decorrelation of the dust polarized emission, in turn originating in the 3D structure of interstellar clouds and their magnetic field. We therefore have to construct a sample of LOSs where dust decorrelation is expected to be significant, and then test whether the EVPA differences between frequencies are larger for that sample than for LOSs where we expect that dust decorrelation is subdominant to effects from the CMB and noise.

The LOS frequency decorrelation of dust polarized emission is more likely to be observed for a given LOS if the following three conditions are met (Tassis & Pavlidou 2015): (i) at least two clouds are present along the LOS and both have a measurable emission contribution; (ii) the mean plane-of-sky magnetic field orientations of the clouds differ by an angle ≳ 60°; and (iii) the SEDs of the clouds are different. The first two conditions imply an emission with polarized intensity weaker than the sum of polarized intensities from individual clouds (LOS depolarization), and a modified polarization angle as compared to the emission from the dominant cloud. The third condition causes the polarization angle to be frequency dependent and is met if the dust clouds have different temperature and/or different polarization spectral index, for example if the dust grain properties differ between clouds.

In this work we rely on the fact that HI column density correlates well with dust in the diffuse ISM (e.g. Boulanger et al. 1996; Planck Collaboration XI 2014; Lenz et al. 2017) and use recent HI datasets to infer whether or not the aforementioned conditions are met.

3 Data sets

In order to identify LOSs where the LOS frequency decorrelation effect is most likely to be significant, we use two types of information that can be extracted from HI observations. The first is the number of clouds along the LOS, obtained via a decomposition of HI spectra by Panopoulou & Lenz (2020). We use publicly available¹ results from this analysisto find sky pixels for which multiple clouds contribute to the dust emission signal in intensity. The second is the plane-of-sky magnetic field orientation as a function of velocity, estimated via the morphology of HI emission by Clark & Hensley (2019). We use publicly available² results from this analysis to further constrain our pixel selection to LOSs that contain clouds with significantly misaligned magnetic fields, that is, the magnetic fields of the clouds form an angle with an absolute value of between 60° and 90°. These HI datasets allow us to define samples of sky pixels with which to study the submillimeter polarized emission as measured by Planck. We concentrate on the high-frequency Planck data, at 217 and 353 GHz, where thermal dust emission is known to dominate the measured polarization signal. In this section we describe the datasets that we use and the post-processing that we apply.

3.1 HI velocity components along the line of sight

If multiple components of dust lie along the LOS, and have different bulk kinematic properties, then the emission spectrum of the HI line willshow multiple peaks at different velocities with respect to the observer. This property of HI emission was used by Panopoulou & Lenz (2020) to measure the number of clouds along the LOS. The authors developed a method to identify the number of peaks in HI spectra and applied it to data from the HI4PI survey (HI4PI Collaboration 2016) over the high-Galactic-latitude sky. The analyzed area covers the parts of the sky where HI column density is well correlated with far-infrared dust emission, as defined by Lenz et al. (2017).

Panopoulou & Lenz (2020) decomposed each HI spectrum into a set of Gaussian components. The Gaussian parameters were grouped within HEALPix pixels of N_side = 128 (termed ‘superpixels’), in order to construct a probability distribution function (PDF) of the centroid velocity of the components. The PDFs were smoothed at a velocity resolution of 5 km s⁻¹. Within each superpixel, clouds were identified as kinematically distinct peaks in the PDF of Gaussian centroid velocity. The Gaussian components belonging to each peak were used to construct a velocity spectrum for each cloud. The published data products include: (a) the column density of each cloud, N_HI, and (b) the first and second moments of the spectrum of each cloud (v₀, σ₀, respectively).

In sightlines with multiple components, not all components will contribute equally to the column density (and similarly to the total dust intensity). Panopoulou & Lenz (2020) introduced a measure of the number of clouds per LOS that takes into account the column densities of clouds, defined as: (1)

where is the column density of the ith cloud in the superpixel and is the column density of the cloud with the highest N_HI in the superpixel. If the column density of a single cloud dominates the total column density of a superpixel, then . If there are two clouds with equal column density, then .

In this paper we use , a map of which is shown in Fig. 1 (top), to distinguish between sightlines whose dust emission is dominated by a single component and those where multiple components might be contributing to the signal. Panopoulou & Lenz (2020) showed that is anticorrelated with the degree of linear polarization at 353 GHz, suggesting that LOSs where multiple components contribute to the polarization signal exhibit larger LOS depolarization than the rest of the sky. However, a simple selection on alone does not imply a high ratio of column densities between clouds; a value of can be achieved by two clouds or by an arbitrary number of clouds, the former case being in general more likely to induce measurable LOS frequency decorrelation. Thus in one variation of our pixel selection we consider a different metric (see Sect. 4) involving the ratio of dominant cloud column densities, , defined as follows: for pixels with at least two clouds (), (2)

where is the column density of the cloud with the highest N_HI, and is that of the cloud with second-highest N_HI. We use the cloud column densities provided by Panopoulou & Lenz (2020).

Fig. 1 Orthographic projections in Galactic coordinates. Longitude zero is marked by the vertical thick lines. The Galactic poles are at the center of each disk. Galactic longitude increases counter-clockwise in the northern hemisphere (left) and clockwise in the southern one (right). We show the maps of the effective number of clouds (top) and the map of Δ(θIVC, θLVC) used in ‘Implementation 1’ (second row). We also show the map of sky positions of pixel samples from ‘Implementation 1’ and ‘Implementation 2’ (bottom). White pixels are both in target1 and target2 samples. Green pixels are target2 pixels not found in target1 and purple pixels are target1 pixels not found in target2. Red pixels belong to the control sample. Black pixels are those that belong to all, but neither to control nor to target1 or target2.

3.2 Orientation of HI structures

The morphology of HI emission encodes properties of the ambient magnetic field in two measurable ways. First, high-resolution HI channel maps reveal thin, linear structures that are well aligned with the magnetic field as traced by starlight polarization (Clark et al. 2014) and polarized dust emission (Clark et al. 2015; Martin et al. 2015). These magnetically aligned HI structures are associated with anisotropic cold HI gas (McClure-Griffiths et al. 2006; Clark et al. 2019; Peek & Clark 2019; Kalberla & Haud 2020; Murray et al. 2020). Second, the degree of alignment of linear HI structures as a function of LOS velocity traces LOS magnetic field tangling, and therefore the observed dust polarization fraction (Clark 2018).

These insights were synthesized into a formalism by Clark & Hensley (2019) that defines 3D maps of the Stokes parameters of linear polarization. These maps are based purely on the morphology of HI emission. The distribution of linear HI emission as a function of orientation on the sky is quantified by the Rolling Hough Transform (RHT; Clark et al. 2014). The RHT is applied to discrete HI velocity channels in an HI data cube to calculate maps of R(v, θ), the linear intensity as a function of LOS velocity v, and orientation θ. Here, R(v, θ) is normalized such that it can be treated analogously to a PDF for the orientation of HI in each pixel. The HI-based Stokes parameters are then defined as: (3) (4)

where I_HI(v) is the HI intensity as a function of LOS velocity. Integrating Q_HI(v) and U_HI(v) over the velocity dimension yields HI-based Stokes Q_HI and U_HI maps that reproduce the Planck 353 GHz Q and U maps with remarkable fidelity. Clark & Hensley (2019) also demonstrate consistency with a tomographic determination of the magnetic field orientation along one LOS based on measurements of optical starlight polarization and Gaia stellar distances (Panopoulou et al. 2019).

We therefore use the Clark & Hensley (2019) maps as a probe of the local magnetic field orientation as a function of LOS velocity. We use their HI4PI-based maps, which use a non-uniform LOS velocity bin size and cover the full sky at the HI4PI angular resolution of 16.2′ (see Clark & Hensley 2019 for map details). To match the resolution and pixelization of the N_c map, we apply a Gaussian filter to degrade the Clark & Hensley maps to a uniform 30′ resolution, and use the healpy function ud_grade to bin the smoothed maps to N_side = 128. We can use these 3D maps to measure the HI-based polarization angle in a specified velocity range by summing Q_HI (v) and U_HI (v) over the desired velocity bins and computing θ_HI = 1∕2 arctan(−U_HI, Q_HI), where arctan is the four-quadrant inverse tangent function here and throughout this paper. In this paper we use θ to denote the position angle of HI structures and ψ for polarization position angles.

3.3 Polarization data from the Planck satellite

In this work we employ two full-sky sets of submillimeter polarization data, both obtained by the Planck satellite. First we utilize the third data release of the Planck collaboration (PR3). We use the 217 GHz single-frequency maps and the 353 GHz single-frequency maps from the polarization-sensitive bolometers only, as recommended in Planck Collaboration III (2020) and Planck Collaboration XII (2020), which we downloaded from the Planck Legacy Archive³ (PLA).

Second, we use a more recent set of high-frequency polarization maps obtained from Planck data but processed through the upgraded map-making algorithm SRoll2, which corrects data for known residual systematic errors in Legacy maps down to the detector noise level (Delouis et al. 2019). We use the full-dataset Polarization Sensitive Bolometers SRoll2 polarization maps at frequency 353 and 217 GHz available at their website⁴ . We note that most of the analysis presented in this paper was completed before the Npipe maps became available (Planck Collaboration Int. LVII 2020). Analyzing this new set of maps would require the implementation of a different analysis pipeline than the one developed and used in this work because per-pixel block-diagonal covariance matrices are not available. However, we note that preliminary studies using Npipe maps yield results that are consistent with those obtained in this paper. As for the case with SRoll2 maps, the detection of LOS frequency decorrelation is more significant with Npipe maps than that obtained using PR3 maps.

We apply the same post-processing to both sets of polarization maps. We smooth the I, Q, and U maps to a resolution of 30′ in order to increase the signal-to-noise ratio. We smooth the per-pixel block-diagonal polarization covariance matrices following the analytical prescription in Appendix A of Planck Collaboration XIX (2015). This formalism neglects correlations between neighboring pixels, but takes into account the off-diagonal covariance between the Q and U Stokes parameters. These terms can be substantial at high Galactic latitudes.

When necessary, we propagate the observational uncertainties in our analysis by making use of Monte Carlo (MC) realizations of correlated noise using a Cholesky decomposition of the smoothed per-pixel block-diagonal covariance matrix (see e.g. Appendix A of Planck Collaboration XIX 2015 or Appendix B of Skalidis & Pelgrims 2019). To assess the observational uncertainty on a measurement, we repeat our analysis on those simulated Stokes parameters and study the resulting per-pixel distribution. We validated this approach by comparing the uncertainties obtained for the polarized intensity and the polarization position angle to analytical estimates.

3.4 Maps of CMB polarization

We make use of the CMB polarization maps obtained from the application to the PR3 data set of the four component-separation algorithms (commander, nilc, sevem, and smica Planck Collaboration XII 2014; Planck Collaboration IX. 2016; Planck Collaboration IV 2020, and references therein).

We downloaded the CMB maps from the PLA and smoothed them so that they all have an effective resolution corresponding to a Gaussian beam with a full width at half maximum (FWHM) of 30′ , just as we do with the single-frequency maps used in this work.

4 Analysis framework

4.1 Sample selection

In order to determine statistically if LOS frequency decorrelation is present and measurable in the Planck high-frequency polarization data, we construct astrophysically selected samples of pixels on the sky based only on HI data. We distinguish between our samples using the labels all (all the pixels in the high Galactic latitude LOS cloud decomposition of Panopoulou & Lenz 2020); control (pixels that should not exhibit LOS frequency decorrelation); and target (pixels that are likely to exhibit large LOS frequency decorrelation). According to Tassis & Pavlidou (2015), the degree of LOS decorrelation between two frequencies depends on (a) how the ratio of polarized intensities contributed by distinct components along a LOS changes between frequencies; and (b) the degree of magnetic field misalignment between these contributing components. The first factor above depends nontrivially on both the temperaturedifference between components, and on the amount of emitting dust (column density) in each. Our physical understanding of these dependencies motivates our definition of control and target samples from HI data:

• control: if the dust emission is strongly dominated by a single component (cloud), no LOS frequency decorrelation is expected, regardless of the other criteria above. For this reason, we construct our control sample using HI data to select pixels where a single component dominates the HI emission (proxy for the emitting dust).

• target: for LOS frequency decorrelation to be significant, there must be (a) more than one contributing component, and (b) significant misalignment (≳60°) between the orientations of plane-of-sky magnetic field that permeate the components. Both criteria are required for a pixel to be included in the target sample. We do not attempt to use the HI data to make predictions about the shape of the dust SED.

The nature of the control sample allows a simple selection criterion: requiring that pixels contain a single cloud along the LOS. We therefore select those pixels that have a column-density-weighted number of clouds (see Sect. 3) equal to unity (). However, for the target sample, there are different ways in which these selection criteria can be implemented in practice. For this reason, we performed the analysis using two distinct implementations of the sample selection, so as to ensure that our particular choices do not qualitatively affect our results. Our selection criteria are described below and are summarized in Table 1 and Fig. 2.

Fig. 2 Cartoon illustration of the pixel selection described in Sect. 4.1. Left panel: HI intensity spectrum of a representative pixel from our control group. The control sample targets sightlines defined by a single HI cloud, parameterized by . Right panel: HI intensity spectrum of a representative pixel that is included in both target1 and target2. Pixels in the target samples are selected to have multiple HI clouds along the LOS, as parameterized by either (target1) or (target2). HI orientations are determined for two clouds along each target LOS by summing the Clark & Hensley (2019) HI-based Stokes parameters over the indicated velocity ranges, and we require that the angles in these clouds differ by at least 60°. Cloud orientations in the target1 sample are determined from predefined IVC and LVC velocity ranges. Cloud orientations in the target2 sample are determined from the 1σ velocity range around the two most prominent HI clouds identified in Panopoulou & Lenz (2020).

Implementation 1

The first criterion for constructing the target sample in this implementation (hereafter target1) selects pixels for which .

This ensures that there is a significant contribution to the dust emission signal in intensity that is not from the dominant component. The same will hold for polarized intensity, with the exception of special cases where the magnetic field in one of the clouds lies mainly along the LOS (which would result in very little, if any, polarized emission from the specific cloud). While we cannot control for the unknown 3D geometry of the magnetic field in each cloud, this unknown simply adds noise to the LOS frequency decorrelation signal we are seeking; our selection of a statistically large sample of pixels likely contains all possible relative orientations between the 3D magnetic field of clouds along the same LOS.

In addition to the requirement that , target1 pixels must also satisfy a misalignment condition. To impose such a condition, we first post-process the Clark & Hensley (2019) HI-based Stokes parameter data (provided in pre-defined discrete velocity bins) to obtain orientation information on a per-cloud basis. We make use of the commonly used distinction of high-latitude HI clouds with respect to their velocity: low-velocity clouds (LVCs) are found in the range − 12 km s⁻¹ ≤ v₀ ≤ 10 km s⁻¹ while intermediate-velocity clouds (IVCs) are found in the range − 70 km s⁻¹ ≤ v₀ ≤−12 km s⁻¹ or 10 km s⁻¹ ≤ v₀ ≤ 70 km s⁻¹ (where v₀ is the cloud centroid velocity and the velocity ranges are defined as in Panopoulou & Lenz 2020). These two classes of clouds are found to show systematic differences in their dust properties, with IVCs, for example, having higher dust temperatures than LVCs on average (e.g., Planck Collaboration XXIV 2011; Planck Collaboration XI 2014; Panopoulou & Lenz 2020). Pixels in which the HI orientation changes significantly between the LVC and IVC ranges likely satisfy all necessary conditions for the LOS frequency decorrelation effect: varying dust SED and magnetic field orientation along the LOS (in addition to the requirement of ).

For each pixel we therefore compute the orientation of two ‘effective’ clouds: an LVC and an IVC. For this, we sum the HI Stokes parameters within the LVC and IVC velocity ranges separately, and then calculate a single HI orientation within the LVC range, θ_LVC, and within the IVC range, θ_IVC. For a pixel to be included in the target1 sample, the misalignment criterion requires that the angles θ_LVC and θ_IVC differ by at least 60°. The (unsigned) angle difference between two angles expressed in radians is computed as (5)

where ξ_1,2 are position angles (either θ’s or ψ’s) defined in the range and where the consecutive absolute values take into account the π degeneracy of orientations.

Implementation 2

We modify the criteria for constructing the target sample in order to test for the robustness of our results against sample selection. First, we identify pixels with at least two significant HI components by requiring (a) that and (b) the ratio of column densities of the two main HI components, , is high (see Eq. (2)). Specifically, candidate pixels for the target sample in this implementation (hereafter target2) are selected so that the column density of the second most prominent component is at least one-third of the dominant component, i.e., . By using instead of a higher threshold in the value of (as was done in Implementation 1) we ensure that the dust emission signal (in intensity at least) arises mainly from two clouds of comparable N_HI , rather than a larger number of low-N_HI clouds (as discussed in Sect. 3).

We also modify the construction of the per-cloud HI orientation compared to Implementation 1. For each cloud, we consider the velocity range within v₀ ± σ₀, where v₀ is the cloud centroid velocity and σ₀ is the second moment of its spectrum. We sum the HI Stokes parameters of the Clark & Hensley maps within this velocity range creating maps of per-cloud Stokes parameters, and . For each pixel, we use these per-cloud Stokes parameters to calculate the HI orientation of the highest-N_HI cloud, θ₁ , and that of the second-highest-N_HI cloud, θ₂ . The target2 sample is constructed by requiring pixels to have Δ(θ₁, θ₂) ≥ 60°, in addition to the aforementioned column-density-based criteria. This cloud-based definition of the misalignment condition avoids relying on the predefined velocity ranges for the LVC and IVC components.

Statistical properties of the samples

The samples contain N_all = 83374, N_control = 7328, N_target1 = 5059, and N_target2 = 5755 high-latitude pixels on a HEALPix map (Górski et al. 2005) of N_side = 128. The pixels in target1 (target2) represent about 6.1% (6.9%) of the high-latitude sky defined by the data and about 2.6% (2.9%) of the full sky. target1 and target2 have 2383 pixels in common. This overlap is to be expected, because despite the different specific criteria, both Implementations 1 and 2 are motivated by the same astrophysical requirements.

In Fig. 1 we show polar projections of the map (top), the difference of position angle between the IVC and LVC effective clouds (second row), and the sky position of the pixels of our control, target1, and target2 samples (bottom). We note that there is a significant difference between the locations of target and control pixels: the former are preferentially found in the northern hemisphere (in both implementations), while the latter are mostly found in the southern hemisphere. This uneven distribution is inherited from the spatial distribution of . As noted in Panopoulou & Lenz (2020), is spatially correlated with the column density of IVCs. The presence of these clouds primarily in the northern hemisphere was previously noted in earlier studies of Galactic HI surveys (e.g., Danly 1989; Kuntz & Danly 1996), and is tied to their astrophysical origin (e.g., Shapiro & Field 1976; Bregman 1980; Wesselius & Fejes 1973; Heiles 1984; Verschuur 1993).

4.2 Statistical methodology

We select pixels from the Planck 353 and 217 GHz polarization maps for each of our three samples, and compute the signed-difference between the EVPAs according to (6)

where the EVPA at both frequencies is determined from the Stokes Q_ν and U_ν according to ψ_ν = 1∕2 arctan(−U_ν, Q_ν) and has a value in the range [0°, 180°). Δ_s (ψ₃₅₃, ψ₂₁₇) is defined in the range [−90°, 90°]. The subscript s in Δ_s is used to denote the signed difference of EVPA from Eq. (5), the unsigned position angle difference (the two are related through Δ(ξ₁, ξ₂) = |Δ_s(ξ₁, ξ₂)|).

We choose to use the signed angle difference rather than the unsigned version because an ensemble of signed angle differences is centered on and symmetric about zero in the absence of systematic offsets. For an ensemble of N 2-circular quantities {ξ_1,2,...,N}, the circular mean and the circular standard deviation are defined as (7)

and (8)

For a sample of pixels, the distribution of Δ_s(ψ₃₅₃, ψ₂₁₇) is expected to have a circular mean close to zero and a finite circular standard deviation. The latter encodes a decorrelation of EVPAs between frequencies due to (i) uncorrelated noise at different frequencies; (ii) the relative contribution of dust and CMB at the two frequencies; and (iii) LOS frequency decorrelation due to the polarized intensity contribution from distinct misaligned dust clouds with SEDs varying between frequencies (the effect we are seeking to detect).

Because the target samples are selected to have a higher likelihood of large LOS frequency decorrelation, we predict a larger circular standard deviation for the target sample than for the control sample. Therefore, we adopt the spread of the distribution of polarization angle differences as our test statistic: (9)

where a detection of LOS frequency decorrelation would correspond to a larger for the target sample than for control. Any inference of the presence of LOS frequency decorrelation has to account for the other sources of increased scatter in the distribution of Δ_s(ψ₃₅₃, ψ₂₁₇), i.e., residual systematic errors, CMB polarization, sampling uncertainties, and data noise, and must consider the possibility that these properties differ between target and control.

We address the first two effects (residual systematic errors and CMB polarization) by repeating our analysis on maps that are derived from the same raw Planck data but processed differently.

One plausible concern is that spatially correlated systematic errors in PR3 maps affect target and control differently, resulting in a false-positive detection of LOS frequency decorrelation. To exclude this possibility, we repeat our analysis using the improved version of Planck HFI polarization maps obtained from the SRoll2 map-making algorithm, which better corrects for known residual systematic errors down to the detector noise level (Delouis et al. 2019). The difference between the Δ_s (ψ₃₅₃, ψ₂₁₇) distributions computed from the PR3 and SRoll2 maps particularized to our samples is shown in Fig. 3. We find that this difference distribution is offset from zero for the target samples, indicating that the region of sky containing the target pixels differed systematically between the PR3 and SRoll2 maps; a conclusion also reached from inspection of Fig. 7 of Delouis et al. (2019).

A second plausible concern is that the contribution of the CMB to the polarized intensity changes between 353 GHz (where it is largely negligible) and 217 GHz (where it might be considerable, especially for pixels with low 217 GHz polarized intensity). This would result in measurable decorrelation of the total emission between 353 and 217 GHz in pixels of low 217 GHz polarized intensity. This is an especially worrisome possibility because target pixels are selected for their misaligned LOS magnetic field structure, and are therefore expected to have systematically lower dust polarized intensity in both frequencies. This is indeed the case, as demonstrated by histograms of the polarized intensities at 353 and 217 GHz in Fig. 4. To exclude the possibility of detecting CMB-induced frequency decorrelation and incorrectly attributing it to frequency decorrelation induced by misaligned magnetic fields in distinct dust components, we perform our analysis on maps from which the CMB contribution has been subtracted. To control against differences between component-separation algorithms, we repeat the analysis on maps obtained using four different algorithms: commander, nilc, sevem, and smica (Planck Collaboration XII 2014; Planck Collaboration IX. 2016; Planck Collaboration IV 2020).

The two remaining effects (sampling uncertainties and data noise) are statistical, and we deal with them through the formulation and statistical testing of two null hypotheses, discussed below. Both null hypotheses express the same physical conclusion: no LOS frequency decorrelation is detectable in Planck data. Rejection of these null hypotheses, consistent across different maps and implementations of the target sample, will constitute evidence for the presence of frequency decorrelation induced by multiple dust components permeated by misaligned magnetic fields along selected LOSs.

We quantify the per-pixel multi-frequency data noise by propagating the observational uncertainties on the individual Stokes parameters at the two frequencies to the measurement of the EVPA difference (). For pixel i we thereforedefine the multi-frequency data noise as (10)

where the ensemble is obtained through the computation of EVPA difference on 10 000 MC simulations of noise-correlated Stokes parameters at each frequency. Therefore, in computing Eq. (10), the sum in Eq. (8) is over realizations, rather than over sample pixels as Eq. (9).

Fig. 3 Normalized histogram of the difference between PR3 and SRoll2 maps of the EVPA difference between 353 and 217 GHz ( for sky pixels of all, control, and target1. For most pixels, the results agree within ~ ± 5°; however, pixels of target1 exhibit larger differences between map versions, centered at 2.3°. This suggests that the sky area covered by our target1 sample received more correction from the systematic cleaning. A similar picture is obtained considering target2 instead of target1.

Fig. 4 Histograms of debiased polarized intensity at 217 GHz (left) and 353 GHz (right) (Plaszczynski et al. 2014) of all (black), control (blue) and target1 (orange). (Dark) Gray shaded areas mark (68) 95 percent of the CMB contribution to the polarized intensity as inferred by smica for a FWHM beam of 30′ and for the footprint. The CMB contribution is negligible at 353 GHz but not at 217 GHz, especially for pixels of low . Histograms correspond to PR3 polarization maps with no CMB subtraction.

4.3 Null hypotheses

Null Hypothesis I

“.” The selection of the target and control samples is astrophysical and “gnagnostic” to other sources of frequency decorrelation. Once the CMB is subtracted, residual systematic errors are corrected, and sample size is accounted for, any significant difference in between the two samples should therefore have an astrophysical explanation. There are two astrophysical reasons why would differ in these samples. First, LOS frequency decorrelation (the effect we are looking for) induces an EVPA change between 217 and 353 GHz in target. This directly increases in target compared to control. Second, LOS frequency decorrelation results in depolarization in target pixels. This increases indirectly in target compared to control, because a lower polarization fraction leads to a lower polarized intensity and thus a higher level of noise (e.g., see Fig. 5). This difference is also attributable to the effect we are looking for. The fact that target pixels are more depolarized than control pixels reflects the anti-correlation between and p₃₅₃ already found in Panopoulou & Lenz (2020). The dissimilarity of p₃₅₃ in the two samples is shown in the left panel of Fig. 6. The misalignment criterion used to select target pixels means that these LOSs experience more LOS depolarization. The polarized intensity and multi-frequency polarization angle uncertainty are anti-correlated (Fig. 5). Thus, the preferentially depolarized target pixels have systematically higher (Fig. 6). We confirmed that there is no systematic difference in the distribution of total intensity between the target and control samples at either frequency.

Therefore, we conclude that, once we have accounted for sample size, any deviation of from zero that persists across all PR3/SRoll2 CMB-subtracted maps should be astrophysical in origin; if the direction of such a deviation is , this would constitute evidence for LOS frequency decorrelation. In practice, we calculate and report: the best-guess value ; its uncertainty, calculated from the individual uncertainties in and ; the p-value of the null hypothesis, . If we were to find the p-value to be improbably low, this would reject the null Hypothesis I and constitute evidence for LOS frequency decorrelation (from a combination of depolarization and direct EVPA change) caused by misaligned magnetic fields in distinct dust components.

Null Hypothesis II

“The observed target sample is a coincidental high-noise draw from the same parent sample as control.” The physical consequence of this hypothesis is that any excess of over is entirely due to target being smaller and noisier⁵ than control (see Figs. 5 and 6); any direct EVPA change between 217 and 353 GHz because of LOS magnetic-field misalignment is below the noise level of Planck data. To test this hypothesis, we generate draws from control that are as small and as noisy as target, and we compare them with the observed target using the test statistic. Clearly, these “target-like” Monte-Carlo-generated draws do not include any EVPA change between 353 and 217GHz due to LOS-frequency decorrelation, because all control pixels feature only a single cloud along that LOS. To match the noise properties of target, we weight the probability of choosing a specific pixel j by its value of , according to the distribution of in target (see Fig. 7). We then construct the distribution of in these simulated target-like LOS-decorrelation–free draws, hereafter referred to as target-like MC, and calculate and report the one-sided p-value of drawing the observed from that distribution (i.e., the probability that in that distribution). If the observed is improbably high compared to typical values in target-like MC (i.e., if the p-value is improbably low), this means that we can reject null Hypothesis II and constitutes evidence for EVPA change due to LOS-induced frequency decorrelation in excess of any increased noise in highly depolarized pixels.

Fig. 5 Two-dimensional normalized histograms of the uncertainties in EVPA differences () and debiased polarized intensity at 353 GHz () for the control sample (top) and the target1 sample (bottom) using PR3 maps with no CMB subtraction. Both histograms are normalized and bounded to the same color scale. The two quantities are correlated: target1 has noisier EVPA differences than control, because of the lower polarized intensities in its pixels.

Fig. 6 Histograms of polarization fraction at 353 GHz (left) and per-pixel inter-frequency uncertainty (, Eq. (10)) (right), for all (black), control (blue), and target1 (orange). Histograms correspond to PR3 polarization maps with no CMB subtraction. The target1 sample is distinctly less polarized (left) and noisier (right) than all and control.

Fig. 7 Effectiveness of weighted resampling in producing target-like MC draws from control with noise properties matched to target: means and standard deviations per bin of normalized histograms of for target1-like MC samples (blue), overplotted on the distribution of those uncertainties for the observed target1 sample (orange). The shaded blue area marks the plus and minus one standard deviation around the mean calculated in each bin from 10000 target1-like MC draws obtained through weighted bootstrap resampling of control. This corresponds to sampling uncertainties. The continuous blue line marks the mean in each bin. Very similar results are obtained forthe target2 sample and for all combinations of sets of polarization maps and CMB estimates.

5 Detection of LOS frequency decorrelation

Figure 8 shows the normalized distributions of Δ_s (ψ₃₅₃, ψ₂₁₇), the signed difference of EVPAs between 353 and 217 GHz frequency bands, for the all, control, and target1 samples in CMB-subtracted PR3 polarization maps. The distributions for control and all are similar, while the distribution for target1 differs noticeably from the other two by being much less peaked around zero and much more spread out.

To test Null Hypothesis I, we calculate , the spread of the distribution of EVPA differences, for target and control, for both implementations of target, both sets of Planck polarization maps, and all CMB estimates from the four component-separation algorithms. We also calculate uncertainties of through unweighted bootstrapping for each of these cases. The left panel of Fig. 9 shows and , with their respective uncertainties, for the PR3 map, from which the smica CMB estimate has been subtracted. It is obvious that is very significantly larger than , and so we expect to be able to reject Null Hypothesis I at very high significance, providing clear evidence forthe presence of LOS frequency decorrelation in Planck data. As the distributions of the and obtained from the bootstrapped samples are very nearly Gaussian, the mean of their difference will be the difference of their means, and the uncertainty of their difference can be obtained from their individual uncertainties added in quadrature. These values are given in Table 2, together with the one-sided p-value of Null Hypothesis I (“”). Indeed, NullHypothesis I is very strongly rejected for both sets of maps (PR3 vs SRoll2), all CMB subtraction algorithms, and both target implementations.

We have therefore established that target has statistically greater polarization angle differences between frequencies than control, and that this is not an effect of increased CMB contribution in target pixels, nor an artifact of the CMB estimate produced by any specific component-separation algorithm, nor an artifact of spatially correlated residual systematic errors.

We now proceed to test whether this excess decorrelation is also significant beyond what would be justified by the increased noise level of target compared to control (i.e., test whether Null Hypothesis II is also rejected). For each case considered, we generated 10 000 target-like MC draws through noise-weighted subsampling from control, as described in the previous section; we calculated the test-statistic for each; we constructed the distribution of ; and we computed the one-sided p-value that describes the probability that the measured for target could be measured for a random pixel sample that is as small as target, as noisy as target, but completely free of LOS decorrelation according to the best current knowledge of the 3D magnetized ISM.

One example of this process is visually represented in the right panel of Fig. 9 for the case of Implementation I of target and PR3 maps from which the smica CMB estimate has been subtracted. It is clear that the observed target is highly decorrelated, even compared to comparably high-noise draws from control. Summary statistics for the distributions and p-values obtained from all samples in both our implementations are reported in Table 3, while a visual representation of these results for all combinations of maps and CMB-subtraction algorithms and for Implementation I of target is shown in Fig. 10. For comparison, both in Table 3 and in Fig. 10, we also providethe results of our analysis on maps without any subtraction of the CMB contribution. Indeed, the p-value of Null Hypothesis II is low for both sets of maps (PR3 vs SRoll2), both implementations of target, and all CMB estimate subtractions, with p-values ranging from 4 × 10⁻³ to 5 × 10⁻² for PR3 maps, and from 4 × 10⁻⁷ to 4 × 10⁻³ for SRoll2 maps.

Null Hypothesis II is systematically rejected at a higher significance for SRoll2 maps than for PR3 maps, if all other features of the analysis remain the same. The most straightforward way to interpret this trend is that residual systematic errors in PR3 maps act as an additional source of noise; the SRoll2 corrections for these systematic errors reduces the noise, and the LOS frequency decorrelation stands out more. Additionally, the significance of the effect is always higher when the CMB has not been subtracted, confirming that the CMB indeed makes a distinct contribution to the difference between 353 and 217 GHz EVPAs, and that difference is more pronounced in the pixels with lower polarization intensity of target.

The robustness of the low p-value of Null Hypothesis II across maps, CMB-subtraction algorithms, and target selections gives us confidence that the effect is real, and that LOS frequency decorrelation due to multiple dust components is present in Planck data and detectable above the noise level – as long as one knows where in the sky to look for it.

Fig. 8 Normalized histograms of Δs(ψ353, ψ217) for the all, control, and target1 samples in black, blue, and orange, respectively. CMB has been subtracted from the PR3 maps using smica. The shaded area results from the propagation of observational uncertainties in Qν and Uν down to the computation of Δs(ψ353, ψ217). The shaded areas mark the plus and minus one standard deviation around the means obtained in each bin of width 2° through the MC simulations. Continuous lines show the means of the three samples.

Fig. 9 Left: rejecting Null Hypothesis I. Summary statistics of values obtained through 10 000 bootstrap resampling of control (cyan) and target1 (orange) samples. The means and one standard deviation are reprensented by the thick vertical lines and shaded area, respectively. Right: rejecting Null Hypothesis II. The blue histogram shows the distribution of values obtained through 10 000 resampling of control with weights that guarantee the same level of EVPA difference uncertainties in the resampled samples as in target1. The shaded blue distribution is a Gaussian fit to the histogram. The vertical orange arrow indicates the value computed for the observed full target1 sample. The examples shown in both panels make use of the PR3 polarizationmaps from which we have subtracted the smica CMB estimate. Results are consistent with all other implementations as shown in Tables 2 and 3.

Fig. 10 Summary statistics of distributions for target1-like simulations obtained from control while subtracting different CMB estimates from the PR3 polarization maps (left) and SRoll2 polarization maps (right) compared to the value of target1. This illustrates part of the information given in Table 3.

6 Validation

In this section we discuss additional validation tests, both statistical and physical, to increase our confidence that we have in fact detected LOS-induced frequency decorrelation in Planck data.

6.1 Sky distribution of target and control pixels

Pixels of target and control samples cover largely disjoint parts of the sky (see Sect. 4.1). It is therefore conceivable that their difference in observed might stem from different local properties, and most notably different instrumental noise or systematic properties of the data. In principle, our test of Hypothesis II, where the observed is compared tothat of noise-matched subsamples of control, should take into account the difference in noise properties; and the comparison between PR3 and SRoll2 maps is performed exactly to evaluate the impact of the residual systematic errors. Nevertheless, we performed two additional tests to verify that some additional, hidden, spatial correlation bias is not generating a false-positive detection of LOS frequency decorrelation.

First, we repeated our analysis inside two sky patches that contain intermixed target1 and control pixels, and are sufficiently small so that instrumental systematic errors should not vary considerably within each patch. These sky patches were defined as regions with an angular radius of 15°, centered on (l, b) = (70°, 50°) in the north, and (l, b) = (−110°, − 50°) in the south. These regions were visually identified and are indicated by green and magenta outlines, respectively, in Fig. 11. These patches as a whole contain 3352 (north) and 3353 (south) pixels. Of those, in the northern (southern) patch, 202 (162) are target1 pixels, and 214 (461) are control pixels. The noise properties within each patch are overall consistent between samples (Fig. 12), unlike the full target and control samples (Fig. 6). We find that for all combinations of maps, CMB subtraction algorithms, and target sample implementations, is larger than . The sample sizes are now too small for Null Hypothesis II to be rejected through the weighted-resampling analysis discussed in Sect. 4.3; we have however verified through subsampling of the full target1 and control samples, that the behavior of both the distribution of and the observed in these sky patches is consistent with what we would expect given the local noise properties and the decrease in sample size.

Second, having observed that noise properties differ systematically between northern and southern hemispheres, we repeated our analysis in the northern hemisphere alone. We chose the northern hemisphere because it contains more target pixels: LOSs intersecting multiple, misaligned clouds are evidently more common in the northern Galactic sky. We find that, despite the modest decrease in sample size for target, in this case the significance with which Null Hypothesis II is rejected increases (p-value decreases), because in general pixels in the north are less noisy.

Fig. 11 Map showing the location of sky pixels belonging to the target1 (white) and control (orange) samples. Black pixels are those in all but neither in target1 nor control. The gray area are pixels where has not be determined (see Panopoulou & Lenz 2020). The locations of the northern and southern sky patches studied in order to investigate the effect of target and control sampling different sky regions are shown with the green and magenta circles, respectively.

Fig. 12 Normalized histograms of as measured on PR3 maps in the northern (left) and southern (right) sky patches for the different subsamples of target1 (in orange), control (in blue), and all (in black).

6.2 Projected Rayleigh statistic

In order to strengthen our analysis and confirm that our results do not depend critically on our choice of as our test statistic, we repeated our analysis using the projected Rayleigh statistic (PRS) to quantify the degree of alignment of EVPA between frequencies. The PRS (Z_x) is computedas (e.g., Jow et al. 2018): (11)

where ξ_i is defined in the range , so that Z_x takes values between −1 and 1. Computing the PRS for a sample of signed difference angles Δ_s(ψ₃₅₃, ψ₂₁₇) defined in Eq. (6), we can quantify the level of alignment of EVPA between 353 and 217 GHz. We expect Z_x to be smaller for samples with statistically larger EVPA differences. We reproduced the analysis presented in Sect. 4 using the PRS in place of the circular standard statistic (S) in order to quantify the degree of alignment/misalignment in our samples and quantitatively compare them. We found that the significance with which our hypotheses are rejected in each case are generally consistent, with no strong dependence on the choice of test statistic.

6.3 versus Δ(θ_LVC, θ_IVC)

According to the simplest two-cloud model (Tassis & Pavlidou 2015), if the EVPA differences between frequencies are due to SED differences and magnetic field misalignment between the dust clouds, then for an ensemble of sky pixels we expect to see (i) a decrease of degree of polarization and (ii) an increase of LOS frequency decorrelation (which we quantify using ) as Δ(θ_LVC, θ_IVC) increases.

To test this simple scenario, we consider all LOSs showing a sufficient degree of complexity in terms of number of clouds, namely . We bin the sky pixels according to their Δ(θ_LVC, θ_IVC) values as measured from HI orientation data in the scheme of Implementation I. Then, for each bin, we examine the distribution of p₃₅₃ and compute the statistic. As expected, for increasing Δ(θ_LVC, θ_IVC), we observe a small but systematic decrease in the degree of polarization and a clear rise of values. The latter is shown in Fig. 13 for the PR3 polarization maps from which the smica CMB was subtracted. We obtain similar conclusions when we use other combinations of polarization maps and removed CMB estimates, as well as when we consider the HI orientation as in Implementation II of the selection of target pixels (i.e., at the peak of the two dominant clouds) rather than the scheme used in Implementation I.

In the simple two-cloud model of Tassis & Pavlidou (2015), LOS decorrelation is expected to be more pronounced towards LOSs where the magnetic fields of the clouds form an angle of 60° or more. As noted by these authors, smaller differences in angle can also result in LOS decorrelation, but at a lower level. LOS decorrelation is therefore not expected to abruptly appear at some large misalignment angle, but should qualitatively match the observed smooth trend in Fig. 13. A more quantitative comparison of this observation with analytic models should take into account a number of factors. First, in general, LOSs might be composed of more than two dust clouds that contribute to the polarized signal. Second, changes in the spectral index of the dust SED (and not simply the dust temperature, as assumed in the Tassis & Pavlidou 2015 model) can alter the frequency dependence of the dust emission EVPA for a given misalignment angle. Finally, the difference between HI filament orientation and the plane-of-the-sky (POS) magnetic field orientation shows an intrinsic astrophysical scatter, which should also be taken into account as an extra source of uncertainty. Such detailed comparisons with models will require further work beyond that presented in this paper.

We note that in our analysis we have not optimized our cutoff in Δ(θ_LVC, θ_IVC) for the selection of our target pixels; rather, we adopted 60° based on our a priori physical expectations. Had we decreased the cutoff to Δ(θ_LVC, θ_IVC) ≥ 45°, the size of the target sample, and hence the significance with which we detected LOS frequency decorrelation, would have increased, as subsequent analysis confirms.

Fig. 13 Increase of the spread of EVPA differences between 353 and 217 GHz as a function of offset angle between HI structures from integration in LVC and IVC ranges (‘Implementation 1’). All sky pixels with Nc > 1.5 are binned according to their Δ(θIVC, θLVC) values and the statistic is computed for each subsample with observational uncertainties propagated. The error bars in each bin represent the 1σ value of from a bootstrap resampling of the data 103 times per bin.

6.4 A case study using starlight polarization

In this study we use HI morphology as an indirect probe of the direction of magnetic fields in individual clouds. Starlight polarization, induced by the same dust grains that produce polarized emission, is a more direct probe of the dust polarization position angle. Currently available starlight polarization measurements are sparse, but large-scale starlight polarization surveys like PASIPHAE (Tassis et al. 2018) are planned for the near future. Nevertheless, data do exist in a small sky patch that we can use for a proof-of-principle analysis using starlight polarization instead of HI data.

Panopoulou et al. (2019) used starlight polarization data from the RoboPol polarimeter (Ramaprakash et al. 2019) to study a sky region where several Galactic dust components are present along the LOS. Based on these stellar polarization data the authors inferred the number of dust clouds and the POS orientation of the magnetic field permeating those for two nearby observing beams of 16′ radius. The latter two were pre-selected based on HI data to likely harbor two dust clouds (two-cloud LOS) and one dust cloud (one-cloud LOS).

The authors demonstrated that the two clouds exhibit significant differences in terms of column density and polarization properties, and that their mean POS magnetic field orientations differ by about 60°. In principle, the different SEDs in those significantly misaligned clouds could lead to a measurable effect of LOS frequency decorrelation in Planck data towards the two-cloud LOS. However, if the effect is weak it could be hidden in the noise, as suggested in Sect. 6.3 of Panopoulou et al. (2019) based on a set of polarization maps from the second Planck data release.

Here, we investigate further the polarization data for those particular LOSs. We retrieve the polarized emission at 217 and 353 GHz measured by Planck towards the sky region of interest (see Sect. 3) smoothed to a 16′ FWHM beam, and we compute the signed difference of EVPA in each pixel (see Eq. (6) in Sect. 4). We thus obtain the EVPA-difference map presented in Fig. 14 (top) where we highlight the two LOSs studied in Panopoulou et al. (2019). Interestingly, the two-cloud region at the center, which is known to feature a complex magnetized ISM structure with at least two dust components (Panopoulou et al. 2019, Clark & Hensley 2019), displays a higher EVPA difference than the nearby one-cloud region.

We quantify the level of uncertainty of the EVPA difference induced by the observational uncertainties on the Stokes Q_ν and U_ν () through MC simulations (see Sect. 3 for details). For each MC draw, we compute the EVPA difference and build the histograms shown in the right panel of Fig. 14. Even when accounting for Planck noise, the two-cloud LOS deviates significantly from zero EVPA difference in the two frequencies, suggesting a LOS frequency decorrelation of the polarization data. In contrast, the distribution corresponding to the one-cloud LOS is compatible with zero EVPA difference in the two frequencies (no LOS frequency decorrelation). Although less significant, the offset from zero of the EVPA difference for the two-cloud LOS survives the subtraction of the CMB estimates. This is reported in Table 4.

This tentative result demonstrates how starlight polarization data can be used to identify sky pixels that experience LOS-induced frequency decorrelation.

Fig. 14 Tomography region of Panopoulou et al. (2019). Top panel: map of EVPA differences computed from 353 and 217 GHz polarization maps from Planck. The two-cloud and one-cloud sight lines are marked respectively by orange and blue crosses northeast and southwest of the map center. The circles have 16′ radius and mark the beams within which starlight polarization data have been taken and studied by Panopoulou et al. (2019). Bottom panel: histograms of EVPA differences computed through 10 000 MC simulations to propagate observational uncertainties on (Qν , Uν ). The two-cloud LOS histogram is shown in orange, the one-cloud LOS in blue. The vertical lines with corresponding colors show the EVPA differences from the data.

7 Estimation of required SED variation

Frequency decorrelation of dust emission is ultimately the result of spatial variations of the dust SED. The detection of rotation of the dust polarization angle between frequencies is evidence for variation of the dust SED along the LOS. We can therefore use the observed magnitude of this effect to estimate the intrinsic variability of the dust SED.

Let us divide the LOS into N clouds such that the ith cloud has column density . Then the observed Stokes parameters of the polarized dust emission at a frequency ν are given by (e.g., Hensley et al. 2019):

where fⁱ, , , , γ_i , and ψ_i are the alignment fraction, dust-to-gas mass ratio, polarized opacity at frequency ν, dust temperature, angle between the magnetic field and the plane of the sky, and polarization angle of the ith cloud, respectively, and m_p is the proton mass. When there are multiple clouds along the LOS, Eqs. (12) and (13) make clear that the ratio U_ν ∕Q_ν, and thus the polarization angle ψ_ν, is generally not constant with frequency.

For a single cloud, the ratios of Q_ν and U_ν at 217 and 353 GHz are given by (14)

If dust everywhere had the same temperature and same opacity law, then this ratio would be constant across the sky and ψ_ν would be constant with frequency. As this is inconsistent with what is observed, let us assume that this quantity has a mean value α and that cloud-to-cloud variations are described by a parameter ρ with a mean of zero, that is, (15)

A modified blackbody with T_d = 19.6 K and β = 1.55, typical parameters for high-latitude dust (Planck Collaboration XI 2020), has α = 0.21, though our analysis is not sensitive to the value of α. Here, σ_ρ quantifies the intrinsic variation in the dust SED between 217 and 353 GHz, regardless of whether those variations arise from temperature, composition, or other effects. Modeling ρ as Gaussian distributed with a mean of zero and variance , we seek the value of σ_ρ that can account for the enhanced dispersion of polarization angles on multi-cloud sightlines (Fig. 8).

To estimate the effect of σ_ρ on the dispersion in polarization angles, we use the HI maps to constrain both the LOS distribution of clouds and their relative orientations. To simplify the analysis we consider the data from our Implementation 2 (Sect. 4). For each sightline we thus consider only the two dominant clouds (in HI column density) as identified by Panopoulou & Lenz (2020) and create maps of per-cloud Stokes parameter by integrating the HI-based Q and U maps of Clark & Hensley (2019) in the velocity range within v₀ ± σ₀, where v₀ is the cloud centroid velocity and σ₀ is the second moment of its spectrum (see Sect. 3). We then estimate ψ₃₅₃ on each sightline as: (16)

where Q_HI and U_HI are given by Eqs. (3) and (4), respectively, and the superscripts denote integration over clouds 1 and 2. The angles between the magnetic field and the plane of the sky γ₁ and γ₂ are unknown, and so we draw sinγ uniformly from the interval [−1, 1] for each; γ = 0 when the magnetic field is in the plane of the sky. This equation does not explicitly model variations in the 353 GHz dust emissivity per H atom, although marginalizing over different values of γ₁ and γ₂ achieves a similar effect numerically. Rather, as we are interested only in the variability of the polarized dust SED between 353 and 217 GHz, we model such effects through the ρ₁ and ρ₂ parameters when computing the 217 GHz polarization angle only.

Using Eq. (15), ψ₂₁₇ on each sightline can be modeled as (17)

On each sightline, ρ₁ and ρ₂ are drawn from a Gaussian distribution of mean zero and variance . For each sightline we can then compute (Eq. (6)) and finally the dispersion (Eq. (9)) over all target2 sightlines as a function of σ_ρ.

If the difference in dispersion between the target and control samples is attributed entirely to varying dust SEDs, then we can estimate (18)

from the resampling analysis presented in Sect. 4 and Table 3. This range is indicated by the horizontal band in Fig. 15.

The ± 1σ range of over 1000 simulations for each value of σ_ρ is presented in Fig. 15. We see that σ_ρ = 0.15 matches the observed enhancement in dispersion between the target and control samples. This is consistent with the dust SEDs varying in the ratio of 217 to 353 GHz polarized intensity at the level of 15% from cloud to cloud over the region analyzed. As we model contributions from only the two most dominant clouds on each sightline, we may be slightly overestimating the true dispersion.

Fig. 15 Dispersion (Eq. (18)) resulting solely from variations in the dust SEDs between two clouds along theLOS in the target2 sample. We quantify the level of SED variation by the parameter σρ (Eq. (15)), finding that σρ = 0.15 can account for the excess dispersion in the target sample. Thus, we estimate that the ratio of 353 to 217 GHz polarized intensity is varying at roughly the 15% level from cloud to cloud. The blue shaded regions indicate the observed range of estimated in Sect. 4 and the ±1σ confidence interval from 1000 realizations of γ1, γ2 , ρ1 , and ρ2 in each pixel. The red shaded region is the resultant constraint on σρ.

8 Discussion

In this paper we report the detection of the effect of LOS-induced frequency decorrelation – the combined effect of varying dust SEDs and magnetic field orientations along the LOS – in Planck polarization data. This detection was made possible by the use of HI datasets, which allowed us to construct our target and control samples a priori in an astrophysically motivated way. The consistency of the results between our two implementations and between the different sets of polarization maps and CMB estimates reinforces our confidence in the robustness of our finding. Our analysis additionally shows that the significance of the effect becomes higher when we use maps cleaned of residual systematic errors that were present in Planck PR3 polarization maps.

We emphasize that we have not in any way optimized our analysis choices to maximize the significance with which the effect is detected. Rather, whenever a choice had to be made, we made it based on astrophysical arguments. There are several examples where different choices in our analysis would have increased the significance of the detection of the effect (decreased the p-value of Null Hypotheses I and II). These include:

(a) Definition of target: defining target as the union of target1 and target2 increases the significance.

(b) Cutoff in HI orientation misalignment: changing the misalignment requirement for inclusion in target from ≥ 60° to ≥ 45° increases the significance.

(c) Localization: restricting our analysis to the northern hemisphere increases the significance.

This work finds evidence for LOS frequency decorrelation, and does not directly address the question of decorrelation in the dust power spectra. Our findings show that frequency decorrelation of the dust polarization signal is not an effect that is uniform throughout the sky, because the change in polarization pattern is more severe for sightlines that pass through more convoluted magnetized ISM, and where those particular sightlines are distributed unevenly on the sky (see the bottom map in Fig. 1). This may have implications in power-spectrum-based estimates; which future work will investigate.

From a CMB perspective, it would be interesting to estimate the level of LOS-induced frequency decorrelation using cross-power spectra (as in, e.g., Planck Collaboration XXX 2016; Planck Collaboration XI 2020) on maps that have been corrected for residual systematic errors (Delouis et al. 2019; Planck Collaboration Int. LVII 2020) and in sky regions thatare dominated by pixels comprising our target samples. Such an analysis is beyond the scope of this paper.

Our Implementation I of the target pixel selection focused on the distinction between LVCs and IVCs based on the physical expectation that IVCs might feature different dust SEDs than LVCs due to differences in temperature and/or dust grain properties (Planck Collaboration XXIV 2011; Planck Collaboration XI 2014). Our Implementation II imposed no such constraint on the velocities of the identified distinct peaks in HI emission. This enables us to use target2 to a posteriori test whether IVC–LVC cloud pairs exhibit a stronger LOS frequency decorrelation effect than LVC–LVC pairs.

Concentrating on the two dominant clouds (the ones corresponding to the two highest-HI-column-density components), we split pixels in target2 in two groups: pixels where both dominant clouds have a velocity centroid in the LVC range (513 pixels), and pixels dominated by LVC-IVC pairs (5242 pixels), with LVC and IVC ranges defined as in Sect. 4.1. We infer the relative strength of LOS frequency decorrelation in the two groups through a uniform, unweighted resampling analysis of each subset of pixels, with N_Boot = 500. The values obtained for PR3 polarization maps and smica CMB subtraction are = 1.05 ± 0.04 and = 1.13 ± 0.05. The statistic is thus found to be higher for LVC–LVC pairs than for IVC–LVC pairs, although the two values are consistent within sampling uncertainties. The same trend is observed for all combinations of polarization maps and subtracted CMB estimate. As a result, it appears that the LOS frequency decorrelation induced by dust clouds does not only involve LOSs passing through IVCs, and, on the contrary, significantly misaligned LVCs may be a substantial source of LOS frequency decorrelation.

Observables that can provide insight into the 3D structure of the magnetized ISM will play a critical role in any effort to improve CMB dust polarization foreground modeling and subtraction accounting for LOS frequency decorrelation. Such observables include HI data (as used in this paper) and starlight polarization. Starlight polarization originates in dichroic absorption by the same dust grains that produce polarized emission, and thus traces the same physical processes of grain alignment with the magnetic field, but for the LOS between observer and star. Large-scale starlight polarization surveys like PASIPHAE (Tassis et al. 2018) will therefore soon provide an independent, direct probe of dust grain orientations in individual clouds.

9 Conclusions

It is certain that the SEDs of the dust clouds vary to some extent between different parts of the Galaxy. It is equally certain that there are in general multiple dust clouds along a large fraction of LOSs, and that the magnetic field of the Galaxy is not uniform andmay vary along the LOS. Consequently, decorrelation between polarized dust emission at different frequencies, both in the plane of the sky and along the LOS, must be present to some extent. The relevant question is whether the magnitude of this frequency decorrelation effect is high enough to be detected by an instrument of given specifications.

In this work, we pursue a new approach that specifically targets LOS frequency decorrelation. Physically, we expect that LOS frequency decorrelation does not occur at a uniform level throughout the sky, but rather should be more severe where the orientations of the magnetic field permeating different dust clouds superposed along the LOS are strongly misaligned. Therefore, we used HI velocity and orientation data to select pixels that are most likely to exhibit significant LOS frequency decorrelation induced by multiple dust SED components. Each of these target sightlines has an HI emission structure consistent with multiple LOS clouds with misaligned magnetic fields. We compare these to a control sample of sightlines that contain a single HI cloud. The use of HI allows us to distinguish these two sets of pixels using data that are entirely independent of polarization measurements. The pixels that maximize the likelihood of showing a LOS frequency decorrelation signal are highly unevenly distributed on the sky.

We quantify LOS decorrelation using the dispersion of inter-frequency EVPA differences. We find that this dispersion is larger for our target sample than for the control sample in Planck data. We detect the LOS frequency decorrelation effect at a level above the Planck noise (see Fig. 9). We confirm that our finding is robust to inhomogeneous data noise level, residual systematic errors, CMB contamination, and the specifics of sky pixel selection. We find that trends in polarization data closely follow the phenomenology expected from the simplest modeling of the effect (Fig. 13). Additionally, relying on a model-independent approach, we estimate that an intrinsic variability of the dust SED of ~ 15% can lead to the observed magnitude of the effect that we measured from polarization maps at 353 and 217 GHz (see Fig. 15). Finally, we demonstrate that LOS superposition of both LVC–LVC and LVC–IVC pairs of clouds contributes to the signal detection.

In this paper we present the first detection of LOS frequency decorrelation in the Planck data. This detection was made possible thanks to the use of ancillary datasets, HI emission data and starlight polarization data, which allow us to identify sky regions that are potentially most susceptible to this effect.

Acknowledgements

We thank Vincent Guillet and Aris Tritsis for insightful discussions. We thank our anonymous referee for her/his report. This work has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement Nos. 771282, 772253, and 819478. G V.P. acknowledges support by NASA through the NASA Hubble Fellowship grant HST-HF2-51444.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. S E.C. acknowledges support by the Friends of the Institute for Advanced Study Membership. V.P. 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 and 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). We acknowledge the use of data from the Planck/ESA mission, downloaded from the Planck Legacy Archive, and of the Legacy Archive for Microwave Background Data Analysis (LAMBDA). Support for LAMBDA is provided by the NASA Office of Space Science. Some of the results in this paper have been derived using the HEALPix (Górski et al. 2005) package. This work is partially based on publicly released data from the HI4PI survey which combines the Effelsberg-Bonn HI Survey (EBHIS) in the northern hemisphere with the Galactic All-Sky Survey (GASS) in the southern hemisphere.

References

All Tables

All Figures

Fig. 1 Orthographic projections in Galactic coordinates. Longitude zero is marked by the vertical thick lines. The Galactic poles are at the center of each disk. Galactic longitude increases counter-clockwise in the northern hemisphere (left) and clockwise in the southern one (right). We show the maps of the effective number of clouds (top) and the map of Δ(θIVC, θLVC) used in ‘Implementation 1’ (second row). We also show the map of sky positions of pixel samples from ‘Implementation 1’ and ‘Implementation 2’ (bottom). White pixels are both in target1 and target2 samples. Green pixels are target2 pixels not found in target1 and purple pixels are target1 pixels not found in target2. Red pixels belong to the control sample. Black pixels are those that belong to all, but neither to control nor to target1 or target2.

In the text

Fig. 2 Cartoon illustration of the pixel selection described in Sect. 4.1. Left panel: HI intensity spectrum of a representative pixel from our control group. The control sample targets sightlines defined by a single HI cloud, parameterized by . Right panel: HI intensity spectrum of a representative pixel that is included in both target1 and target2. Pixels in the target samples are selected to have multiple HI clouds along the LOS, as parameterized by either (target1) or (target2). HI orientations are determined for two clouds along each target LOS by summing the Clark & Hensley (2019) HI-based Stokes parameters over the indicated velocity ranges, and we require that the angles in these clouds differ by at least 60°. Cloud orientations in the target1 sample are determined from predefined IVC and LVC velocity ranges. Cloud orientations in the target2 sample are determined from the 1σ velocity range around the two most prominent HI clouds identified in Panopoulou & Lenz (2020).

In the text

Fig. 3 Normalized histogram of the difference between PR3 and SRoll2 maps of the EVPA difference between 353 and 217 GHz ( for sky pixels of all, control, and target1. For most pixels, the results agree within ~ ± 5°; however, pixels of target1 exhibit larger differences between map versions, centered at 2.3°. This suggests that the sky area covered by our target1 sample received more correction from the systematic cleaning. A similar picture is obtained considering target2 instead of target1.

In the text

Fig. 4 Histograms of debiased polarized intensity at 217 GHz (left) and 353 GHz (right) (Plaszczynski et al. 2014) of all (black), control (blue) and target1 (orange). (Dark) Gray shaded areas mark (68) 95 percent of the CMB contribution to the polarized intensity as inferred by smica for a FWHM beam of 30′ and for the footprint. The CMB contribution is negligible at 353 GHz but not at 217 GHz, especially for pixels of low . Histograms correspond to PR3 polarization maps with no CMB subtraction.

In the text

Fig. 5 Two-dimensional normalized histograms of the uncertainties in EVPA differences () and debiased polarized intensity at 353 GHz () for the control sample (top) and the target1 sample (bottom) using PR3 maps with no CMB subtraction. Both histograms are normalized and bounded to the same color scale. The two quantities are correlated: target1 has noisier EVPA differences than control, because of the lower polarized intensities in its pixels.

In the text

	Fig. 6 Histograms of polarization fraction at 353 GHz (left) and per-pixel inter-frequency uncertainty (, Eq. (10)) (right), for all (black), control (blue), and target1 (orange). Histograms correspond to PR3 polarization maps with no CMB subtraction. The target1 sample is distinctly less polarized (left) and noisier (right) than all and control.
In the text

Fig. 7 Effectiveness of weighted resampling in producing target-like MC draws from control with noise properties matched to target: means and standard deviations per bin of normalized histograms of for target1-like MC samples (blue), overplotted on the distribution of those uncertainties for the observed target1 sample (orange). The shaded blue area marks the plus and minus one standard deviation around the mean calculated in each bin from 10000 target1-like MC draws obtained through weighted bootstrap resampling of control. This corresponds to sampling uncertainties. The continuous blue line marks the mean in each bin. Very similar results are obtained forthe target2 sample and for all combinations of sets of polarization maps and CMB estimates.

In the text

Fig. 8 Normalized histograms of Δs(ψ353, ψ217) for the all, control, and target1 samples in black, blue, and orange, respectively. CMB has been subtracted from the PR3 maps using smica. The shaded area results from the propagation of observational uncertainties in Qν and Uν down to the computation of Δs(ψ353, ψ217). The shaded areas mark the plus and minus one standard deviation around the means obtained in each bin of width 2° through the MC simulations. Continuous lines show the means of the three samples.

In the text

Fig. 9 Left: rejecting Null Hypothesis I. Summary statistics of values obtained through 10 000 bootstrap resampling of control (cyan) and target1 (orange) samples. The means and one standard deviation are reprensented by the thick vertical lines and shaded area, respectively. Right: rejecting Null Hypothesis II. The blue histogram shows the distribution of values obtained through 10 000 resampling of control with weights that guarantee the same level of EVPA difference uncertainties in the resampled samples as in target1. The shaded blue distribution is a Gaussian fit to the histogram. The vertical orange arrow indicates the value computed for the observed full target1 sample. The examples shown in both panels make use of the PR3 polarizationmaps from which we have subtracted the smica CMB estimate. Results are consistent with all other implementations as shown in Tables 2 and 3.

In the text

	Fig. 10 Summary statistics of distributions for target1-like simulations obtained from control while subtracting different CMB estimates from the PR3 polarization maps (left) and SRoll2 polarization maps (right) compared to the value of target1. This illustrates part of the information given in Table 3.
In the text

Fig. 11 Map showing the location of sky pixels belonging to the target1 (white) and control (orange) samples. Black pixels are those in all but neither in target1 nor control. The gray area are pixels where has not be determined (see Panopoulou & Lenz 2020). The locations of the northern and southern sky patches studied in order to investigate the effect of target and control sampling different sky regions are shown with the green and magenta circles, respectively.

In the text

	Fig. 12 Normalized histograms of as measured on PR3 maps in the northern (left) and southern (right) sky patches for the different subsamples of target1 (in orange), control (in blue), and all (in black).
In the text

Fig. 13 Increase of the spread of EVPA differences between 353 and 217 GHz as a function of offset angle between HI structures from integration in LVC and IVC ranges (‘Implementation 1’). All sky pixels with Nc > 1.5 are binned according to their Δ(θIVC, θLVC) values and the statistic is computed for each subsample with observational uncertainties propagated. The error bars in each bin represent the 1σ value of from a bootstrap resampling of the data 103 times per bin.

In the text

Fig. 14 Tomography region of Panopoulou et al. (2019). Top panel: map of EVPA differences computed from 353 and 217 GHz polarization maps from Planck. The two-cloud and one-cloud sight lines are marked respectively by orange and blue crosses northeast and southwest of the map center. The circles have 16′ radius and mark the beams within which starlight polarization data have been taken and studied by Panopoulou et al. (2019). Bottom panel: histograms of EVPA differences computed through 10 000 MC simulations to propagate observational uncertainties on (Qν , Uν ). The two-cloud LOS histogram is shown in orange, the one-cloud LOS in blue. The vertical lines with corresponding colors show the EVPA differences from the data.

In the text

Fig. 15 Dispersion (Eq. (18)) resulting solely from variations in the dust SEDs between two clouds along theLOS in the target2 sample. We quantify the level of SED variation by the parameter σρ (Eq. (15)), finding that σρ = 0.15 can account for the excess dispersion in the target sample. Thus, we estimate that the ratio of 353 to 217 GHz polarized intensity is varying at roughly the 15% level from cloud to cloud. The blue shaded regions indicate the observed range of estimated in Sect. 4 and the ±1σ confidence interval from 1000 realizations of γ1, γ2 , ρ1 , and ρ2 in each pixel. The red shaded region is the resultant constraint on σρ.

In the text