Magnetic Misalignment of Interstellar Dust Filaments

We present evidence for scale-independent misalignment of interstellar dust filaments and magnetic fields. We estimate the misalignment by comparing millimeter-wave dust-polarization measurements from Planck with filamentary structures identified in neutral-hydrogen (HI) measurements from HI4PI. We find that the misalignment angle displays a scale independence (harmonic coherence) for features larger than the HI4PI beam width ($16.2'$). We additionally find a spatial coherence on angular scales of $\mathcal{O}(1^\circ)$. We present several misalignment estimators formed from the auto- and cross-spectra of dust-polarization and HI-based maps, and we also introduce a map-space estimator. Applied to large regions of the high-Galactic-latitude sky, we find a global misalignment angle of $\sim 2^\circ$, which is robust to a variety of masking choices. By dividing the sky into small regions, we show that the misalignment angle correlates with the parity-violating $TB$ cross-spectrum measured in the Planck dust maps. The misalignment paradigm also predicts a dust $EB$ signal, which is of relevance in the search for cosmic birefringence but as yet undetected; the measurements of $EB$ are noisier than of $TB$, and our correlations of $EB$ with misalignment angle are found to be weaker and less robust to masking choices. We also introduce an HI-based dust-polarization template constructed from the Hessian matrix of the HI intensity, which is found to correlate more strongly than previous templates with Planck dust $B$ modes.


MOTIVATION
We continue the investigation of magnetic misalignment from Clark et al. (2021), which sought an explanation for the parity-violating T B correlation measured in Galactic dust polarization by the Planck satellite at millimeter wavelengths (Planck Collaboration et al. 2020a). A polarization field can be generically decomposed into parity-even E modes and parity-odd B modes (Seljak & Zaldarriaga 1997;Kamionkowski et al. 1997). The T B cross-spectrum is a measure of the correlation between the total intensity T and the B-mode polarization and indicates a net chirality in the polarization field. The T E cross-spectrum is a correlation with the E-mode polarization and is non-chiral. Using the dust-dominated frequency channel centered at 353 GHz, Planck Collaboration et al. (2020a) reported T B/T E ∼ 0.1 in the multipole range 40 < < 600, which roughly corresponds to angular scales of 1-10 • . Parity-probing cross-spectra such as T B and T E are of interest both in studies of the interstellar medium (ISM), for which the observed cross spectra may constrain magnetohydrodynamic (MHD) models (Caldwell et al. 2017;Kandel et al. 2017;Kritsuk et al. 2018;Kim et al. 2019), and in measurements of the cosmic microwave background (CMB), for which asymmetries in the Galactic foregrounds can bias polarization calibration (Abitbol et al. 2016) and confound searches for cosmic birefringence (Minami & Komatsu 2020). A magnetic helicity in the local ISM (Brandenburg & Subramanian 2005;Blackman 2015) could produce a nonzero T B correlation, and Bracco et al. (2019) produced toy models with positive T B and T E on large scales (multipoles 20). The polarization of interstellar dust emission is a probe of Galactic magnetic fields, the observed polarization orientation being perpendicular to the plane-ofsky (POS) magnetic field (Stein 1966;Hildebrand 1988;Martin 2007). At the same time, the dust in the diffuse ISM is organized partially in filamentary structures that are preferentially aligned to the magnetic field (Planck Collaboration et al. 2016a,b). Filamentary structures can also be identified in neutral hydrogen (Hi), which is well mixed with dust (Lenz et al. 2017) and has the advantage of three-dimensional information from spectroscopic separation into velocity bins. The alignment between Hi filaments and magnetic-field orientations has been confirmed by comparison to millimeter-wave polarization and optical starlight polarization (McClure-Griffiths et al. 2006;Clark et al. 2014Clark et al. , 2015Martin et al. 2015;Kalberla et al. 2016).
Whereas previous work, e.g., Clark & Hensley (2019), has assumed a perfect alignment between interstellar dust filaments and magnetic-field lines, it was suggested in Huffenberger et al. (2020) that a small misalignment could act as a mechanism for parity violation, i.e., a tendency toward features of one chirality (or handedness) over the other. In Clark et al. (2021), this idea was extended to allow dust filaments and magnetic-field orientations to display a scale-dependent misalignment, which could potentially account for the observed T B.
In this work, we directly compute the misalignment angle in many regions of the sky and in many multipole bins for > 100. We find evidence for scale independence of the misalignment angle and also for a correlation with the observed dust T B.

Observed T B
The dust T B was reported in Planck Collaboration et al. (2016c), where it was noted that a positive signal in the multipole range 60 < < 130 became more significant as the sky area was increased. The investigation was continued in Planck Collaboration et al. (2020a) with the observation that T B/T E ∼ 0.1 for 40 < < 600. The EB signal was reported to be consistent with null.
In Weiland et al. (2020), the T B signal was further confirmed by using WMAP K-band polarization (Page et al. 2007) in place of the Planck B modes and also by using the magnetic-field template from Page et al. (2007) that is based on optical starlight-polarization catalogs (Heiles 2000;Berdyugin et al. 2001;Berdyugin, A. & Teerikorpi, P. 2002;Berdyugin, A. et al. 2004). The K-band measurement is dominated by synchrotron rather than dust emission but is also a probe of Galactic magnetic fields. The starlight measurements largely probe the same magnetic dust-grain alignment that produces polarized millimeter-wave emission. Both choices are independent of the Planck polarization calibration, and both show positive T B.

Magnetic misalignment
Magnetic misalignment is a discrepancy between the orientation of filamentary density structures and the polarization-inferred magnetic-field lines. In the case of perfect magnetic alignment, we expect T E > 0 and T B = 0 (Zaldarriaga 2001). A misalignment of 45 • would produce T E = 0 and T B = 0, where the sign depends on the chirality of the misalignment. The robustly positive T E measured by Planck can be interpreted as supportive evidence for magnetic alignment of dust filaments (Clark et al. 2015;Planck Collaboration et al. 2016b;Kalberla et al. 2016). In Planck Collaboration et al. (2020a), the T E correlation over sky regions and multipoles is reported as r T E = D T E / D T T D EE = 0.357 ± 0.003, where D XY denotes the cross-spectrum of X and Y ; the T B correlation is reported as r T B = D T B / D T T D BB ≈ 0.05. Since the T B correlation is much smaller than the T E correlation, the magnetically aligned model need only be perturbed a small but coherent amount in order to produce the observed T B, and this perturbation would also produce a positive EB (Huffenberger et al. 2020), though this EB would be obscured by Planck noise (Clark et al. 2021).
In Clark et al. (2021), it was suggested that an Hi-based filamentary polarization template could be used as a comparison point in the search for magnetic misalignment. The template of Clark & Hensley (2019) is constructed by 1) quantifying the orientation of linear Hi structures with the Rolling Hough Transform (RHT, Clark et al. 2014) in velocity-channel maps from Hi4PI (HI4PI Collaboration et al. 2016), 2) assuming perfect alignment between the RHT-measured Hi orientation and the POS magnetic-field orientation and thereby obtaining a prediction for the dust polarization angle, 3) applying weights based on the Hi intensity, and 4) integrating the channel maps to form a template that can be compared to the measured millimeter-wave dust polarization. A strong correlation with the Planck 353-GHz maps is detected in both E and B modes up to the Hi4PI beam scale of 16.2 ( 1000). The Hi4PI-based template was used in Clark et al. (2021), but Clark & Hensley (2019) also constructed polarization templates with observations from the Galactic Arecibo L-Band Feed Array Hi Survey (GALFA-Hi, Peek et al. 2018), which has higher angular resolution (FWHM = 4.1 ) but smaller sky coverage (32% of the celestial sphere).
Given the alignment between Hi and dust filaments, a difference between the Planck -measured dust polarization angles and the Hi-inferred angles is a potential indication of magnetic misalignment and could be used as a tracer of dust T B and EB (Clark et al. 2021). In extending the work of Clark et al. (2021), we measure the aggregate misalignment angle in different sky regions by using the Hi template as a reference. We study the observed properties of this misalignment angle and its correlation with the measured dust T B and EB. Clark et al. (2021) also introduced a scale-dependent effective misalignment angle ψ , which is a function of multipole . This effective misalignment angle is given explicitly by (cf. Eq. 11 of Clark et al. 2021) where we see that the ratio T B/T E is the controlling quantity. 1 As noted in Planck Collaboration et al. (2020a), the ratio T B/T E is approximately constant across a broad range of multipoles at high Galactic latitudes, and this is related to the observation of Clark et al. (2021) that ψ ∼ 5 • in the range 100 500 on a similar sky area.
Equation 1 provides an estimator for the effective misalignment angle. Because it is formed from the T B and T E cross-spectra, we will call this type of estimator "spectrum-based". In this work, we present several additional spectrum-based estimators by considering the auto-and cross-spectra of the Planck dust maps and the Hi templates. We also present a map-based estimator that is similar to the projected Rayleigh statistic (PRS) of Jow et al. (2017). We test for consistency among these estimators.
Although Clark et al. (2021) allowed for scale dependence in the misalignment angle, we find in this work that ψ tends to display a scale independence even when measured on small regions of sky. Equivalently, we find that ψ is roughly constant with , which we will occasionally refer to as harmonic coherence.
It is important to note that the dust is likely organized only partially in filaments, which are in turn only partially captured by the Hi template. We expect, therefore, that there are contributions to the dust polarization that are unrelated to the Hi template and, more generally, unrelated even to the true underlying filamentary structure. An estimator like that of Eq. 1 may be influenced by these non-filamentary contributions, since it depends only on the T B and T E cross-spectra of the full dust maps. Some of the estimators we will introduce in later sections will be defined by reference to the Hi template, which will partially but imperfectly restrict the analysis to modes which are related to filaments.
In contrast to the previous paragraph, the DUSTFILAMENTS code of Hervías-Caimapo & Huffenberger (2022) constructs a phenomenological dust model, which is composed entirely of filaments and which reproduces the main features of the angular power spectra measured by Planck. Using this model, it was recently shown in Huang (2022) that the measured T B is unlikely to be a statistical fluctuation of an underlying parity-even distribution, if the assumptions of the DUSTFILAMENTS code represent the true sky.

Cosmic birefringence
Cosmic birefringence is an observable consequence of certain types of parity-violating physics beyond the Standard Model and manifests as a rotation of the plane of linear polarization of photons (Carroll et al. 1990;Harari & Sikivie 1992;Carroll 1998). A popular source of cosmic birefringence is an electromagnetically-coupled axion-like field, which can behave as both dark matter and dark energy (Marsh 2016). In the CMB, the polarization rotation can be detected as an EB correlation (Lue et al. 1999;Feng et al. 2005Feng et al. , 2006Liu et al. 2006). A T B correlation should also be produced, but it is typically a less sensitive observable on account of the large cosmic variance in T .
There are several species of cosmic birefringence that have been investigated in the literature. An isotropic, static cosmic birefringence manifests as an overall polarization rotation by the same angle along every line of sight. This observable is, unfortunately, degenerate with a miscalibration of the instrumental polarization orientation (Yadav et al. 2010). The degeneracy is sometimes exploited as a means of "self-calibration" by assuming a standard cosmology in which the true EB vanishes (Keating et al. 2013). Although this type of calibration removes sensitivity to an isotropic, static cosmic birefringence, it is still possible to search for cosmic birefringence which is anisotropic (Ade et al. 2015;BI-CEP2 Collaboration et al. 2017;Namikawa et al. 2020;Bianchini et al. 2020) or time-variable (BICEP/Keck et al. 2021, 2022Ferguson et al. 2022). Through a campaign of modeling and calibration, it is possible to account for instrumental systematics and measure the isotropic, static cosmic-birefringence angle. Recent measurements of this kind are consistent with a standard cosmology (Kaufman et al. 2014;Gruppuso et al. 2016;Planck Collaboration et al. 2016d;Choi et al. 2020).
A new technique was proposed in Minami et al. (2019), which exploits the fact that the Galactic foregrounds are subject only to polarization miscalibration and not to cosmic birefringence. The observed CMB is rotated by both miscalibration and a possible cosmic birefringence. With measurements at multiple frequencies, the calibration angles and the cosmic-birefringence angle can be extracted simultaneously. Applied to Planck 2018 polarization data (Planck Collaboration et al. 2020b), a cosmic-birefringence angle β = 0.35 ± 0.14 • , a discrepancy with the null hypothesis with a significance of 2.4σ, was reported in Minami & Komatsu (2020) under the assumption of a vanishing dust EB. With the newer Planck maps produced by the NPIPE pipeline (Planck Collaboration et al. 2020c), the same prescription produced β = 0.30 ± 0.11 • as reported in Diego-Palazuelos et al. (2022). Recently, a similar analysis that includes WMAP polarization data (Bennett et al. 2013) produced the consistent but stronger result Eskilt & Komatsu 2022). In these two recent cosmic-birefringence analyses, the impact of a possible foreground EB correlation was incorporated by two different approaches, one of which was based on the filamentary misalignment paradigm of Huffenberger et al. (2020) and Clark et al. (2021). When accounting for a possible foreground EB, the birefringence angle varies as a function of sky fraction but remains positive. Diego-Palazuelos et al. (2022) refrain from an estimate of statistical significance due to the currently limited understanding of foreground polarization, while Eskilt & Komatsu (2022) quote a significance of 3.6σ but acknowledge that the foreground polarization must be better understood to be confident that the measured EB is cosmological rather than Galactic. The study of magnetic misalignment is, therefore, of central importance in the search for cosmic birefringence.

Outline
In Sec. 2, we describe the data products used throughout the analysis. In Sec. 3, we introduce a new filamentary polarization template that relies on the Hessian matrix of Hi intensity maps. In Sec. 4, we present our misalignment ansatz, i.e., our assumptions of how misalignment perturbs the dust polarization in both map space and harmonic space. We derive misalignment estimators in terms of the auto-and cross-spectra of the Planck dust maps and the Hi-based polarization templates, and we test some immediate consequences of these relations. In Sec. 5, we describe a set of mock skies that we have used to check our estimators. These mock skies are constructed to match the 2-point statistics of the Planck dust maps including cross-spectra with the Hi template. In Sec. 6, we introduce a map-based misalignment estimator and present tentative evidence for a global misalignment angle of ∼ 2 • . In Sec. 7, we divide the sky into small patches and present evidence for scale independence (harmonic coherence) of magnetic misalignment as well as evidence of spatial coherence. In Sec. 8, we present evidence for a scale-independent relation between magnetic misalignment and parity-violating cross-spectra such as T B and EB. We close in Sec. 9 with suggestions for improvements in our analysis and new directions to further the investigation of parity violation in Galactic dust polarization.

DATA
We use the Planck Commander dust maps (Planck Collaboration et al. 2020d) as our fiducial measurements of the on-sky thermal dust emission in Stokes T , Q and U . The maps are constructed by component separation applied to the nine Planck frequency maps, whose passband centers span 30-857 GHz, though polarization is available only for the seven bands spanning 30-353 GHz.
Half-mission maps are available and are constructed using data exclusively from either the first or the second half of the Planck observation period. When forming cross-spectra, we will often use these half-mission splits in order to avoid positive-definite noise biases.
In this work, our Hi template is derived from the Hi4PI survey (HI4PI Collaboration et al. 2016), a set of full-sky maps of the 21-cm hyperfine transition with an angular resolution of 16.2 , a sensitivity σ rms = 43 mK and a velocity (spectral) resolution of δv = 1.49 km/s. The Hi4PI survey is a combination of the Parkes Galactic All-Sky Survey (GASS, McClure-Griffiths et al. 2009) and the Effelsberg-Bonn Hi Survey (EBHIS, Winkel et al. 2016). The GASS observations cover the southern sky in the velocity range −470 ≤ v lsr ≤ 470 km/s, and the EBHIS observations cover the northern sky in the velocity range −600 ≤ v lsr ≤ 600 km/s. At high Galactic latitudes, nearby dust is generally expected to be associated with lower-velocity Hi emission, i.e., with small |v lsr |, so our dust-polarization template is drawn from the range −15 ≤ v lsr ≤ 4 km/s, a choice that is motivated in more detail in Sec. 3.1.
We compute purified power spectra with NaMaster (Alonso et al. 2019) using a C 2 apodization window (Grain et al. 2009) with a scale of 1 • . Before computing power spectra, we smooth the Commander maps to 16.2 , the Hi4PI beam width. We use a HEALPix pixelization scheme (Górski et al. 2005) and downgrade all maps to N side = 256 for faster power-spectrum estimation. We spot-checked some of our results at higher N side and find that they are consistent.

Galaxy masks
We use the Galaxy masks provided by the Planck Legacy Archive. 2 These masks are constructed to limit Galactic emission to varying levels. The masks with smaller sky fraction f sky restrict the analysis to relatively high latitudes. The masks with larger f sky allow more contributions from nearer the Galactic plane. The set of Galaxy masks is shown in Fig. 1. Our fiducial mask in much of the analysis is defined by f sky = 70%, and we will refer to it as the "70% Galaxy mask".

Notation
We use the subscript "Hi" to denote quantities derived from the Hi-based polarization template. For example, the Hi-based prediction for dust E modes is denoted by E Hi . It is important to note that these quantities are describing Hi-based predictions for the polarization of dust rather than polarization properties of the Hi itself. The Hi is measured in total (unpolarized) intensity, and prescriptions like the Hessian method of Sec. 3 convert those intensity maps into dust-polarization templates.
We use the subscript "d" to denote quantities related to Galactic dust. Usually, this will refer specifically to the Planck Commander maps described above.

Bandpass filtering
Much of our analysis is restricted to > 100, and we often form maps which are bandpass filtered. We filter by applying an -dependent Tukey window to the spherical-harmonic representation of the maps. We use a taper of length t = 50, which produces a flat-topped passband when the window width is larger than 2t . We tested these filters on full-sky Planck dust maps and find the out-of-band response to be suppressed by a factor of more than 10 4 . In particular, the out-of-band leakage is below the level of the high-latitude dust power (computed on the Planck 70% Galaxy mask), even when the filtered power spectra are computed on the full sky, i.e., including the Galactic Plane.

HESSIAN METHOD
We introduce a new Hessian-based filament-finding algorithm (similar to those of, e.g., Planck Collaboration et al. 2016b; Kalberla et al. 2021). Whereas previous work on misalignment (Clark et al. 2021) used a filamentary model based on the Rolling Hough Transform (RHT, Clark et al. 2014;Clark & Hensley 2019), we find that our new Hessian-based polarization template correlates more strongly with Planck measurements of B-mode dust polarization for 100 (Fig. 17). Furthermore, whereas the RHT loses its correlation with Planck B modes for 400, the Hessian maintains a correlation of ∼ 10% up to our highest multipoles ( max = 767). In E modes, the two methods correlate with Planck at roughly equivalent levels.
We additionally prefer the Hessian method for its relative computational efficiency. The Hessian method requires only two operations in spherical-harmonic space, while the RHT requires a suite of convolutions to sample polarization angles. Direct comparisons will be presented in future work (Halal et al. in prep.).
The Hessian matrix contains information about the local second derivatives. By searching for regions of negative curvature in an intensity map, we find candidate filaments. Negative curvature implies that at least one of the Hessian eigenvalues is negative. The orientation of the filament is determined by the local eigenbasis. As in, e.g., Clark et al. (2015), we assume that the planeof-sky (POS) filament is aligned with the POS magnetic field. The dust polarization is taken to be orthogonal to the filament. With these assumptions, we can convert an intensity map into a polarization template.
Hessian-based filament identifications have been performed, e.g., on 353-GHz maps in Planck Collaboration et al. (2016a,b), on Hi4PI Hi and Planck 857-GHz maps in Kalberla et al. (2021), on Herschel images of molecular clouds (Polychroni et al. 2013) and on simulations of the cosmic web (Colombi et al. 2000;Forero-Romero et al. 2009). In addition to the RHT, some non-Hessian filament-finding algorithms that have been applied to studies of the ISM include DisPerSE (Sousbie 2011;Arzoumanian et al. 2011) andgetfilaments (Men'shchikov 2013). See Sec. 3.10 of Hacar et al. (2022) for a more comprehensive review.

Prescription
We use the Hessian matrix to identify filament orientations. To construct Hi-based templates for dust polarization, we form weights from the Hessian eigenvalues.
We analyze the Hi maps in individual velocity bins. Our final polarization template is produced by summing over velocities. The Hi intensity in velocity channel v is denoted I v . We work in spherical coordinates with polar angle θ and azimuthal angle φ. The local Hessian matrix is given by where The eigenvalues are where The candidate polarization angle is then but we will enforce conditions below to ensure this identification is sensible. First, for the local curvature to be negative along at least one axis, we need λ − < 0. Second, we want this negative curvature to be the dominant local morphology, so we require λ − to be the larger of the two eigenvalues in magnitude. Define Then we define the velocity-dependent weight where the Iverson brackets on the right-hand side enforce conditions on ∆λ and λ − . 3 Our velocity-dependent Stokes templates are given by The Hessian method is susceptible to small-scale noise and scan artifacts, so we restrict our analysis to the Hi4PI velocity bins with greatest sensitivity. We start with the binning of Clark & Hensley (2019). As a proxy for noisiness, we search for pixels with intensities that are reported to be negative. We remove any velocity slice that contains negative-intensity pixels on the Planck 70% Galaxy mask. This leaves a continuous range between −15 and 4 km/s. The velocity selection is intended to avoid numerical pathologies and should be revisited in a future iteration of the Hessian algorithm.
Most of the Hi emission is at low velocities (see, e.g., Fig. 1 of Clark & Hensley 2019), so we are retaining the dominant contributions even with the current velocity cuts. The velocity cut affects our analysis mainly in terms of sensitivity, since there is potentially useful information about dust filaments in the velocity slices that are discarded. In general, sensitivity is greater when the Hi template correlates more strongly with Planck dust polarization. We defer to future work an investigation of the potential improvements from differently chosen velocity cuts (Halal et al. in prep.). We repeated a majority of the following calculations with the RHT-based template (Sec. A) that uses the much broader velocity range of −90 to 90 km/s as in Clark & Hensley (2019), and we find consistent results. The full templates are given by for X ∈ {T, Q, U }. An illustration of our Hessian method is provided in Fig. 2, where we analyze a region with area 10 • ×10 • centered on (l, b) = (12 • , 45 • ). The velocity bin is centered on −4.4 km/s with a width of 1.3 km/s. The panels of Fig. 2 show how the raw Hi4PI intensity map is transformed into a filamentary intensity w v and how the filament orientations determine the inferred magnetic-field orientations.
Additional material related to our Hessian method is provided in Sec. A.

MISALIGNMENT ANSATZ
As an ansatz for the observable signature of magnetic misalignment, we assume a multipole-dependent rotation angle ψ as in Clark et al. (2021). We denote the observed E and B modes by E( ) and B( ) and the unphysical modes that would be observed in the absence of misalignment byẼ( ) andB( ), where identifies a particular spherical harmonic with multipole moment . Our ansatz takes the form (15) For the purposes of the ansatz, we are imagining that all of the dust polarization participates in the misalignment. As mentioned in Sec. 1.2, this assumption is likely inaccurate, since some of the dust morphology is nonfilamentary. Later, we will form estimators by comparing the observed dust polarization with the predictions of the Hi template, and this will restrict the analysis to the filamentary modes that we do expect to be described by the ansatz of Eq. 15 (in the misalignment paradigm).
To make magnetic misalignment less abstract, we provide an illustration in Fig. 3. We consider perfect alignment (black) and scale-independent misalignment (green). Perfect alignment is assumed by the Hibased filamentary model of dust polarization (Sec. 1.2). For scale-independent misalignment, we take ψ = 20 • for all . This is a much larger misalignment than we expect to measure on the true sky, but the exaggeration is useful for visualization. In this case, the magnetic field shows a consistent rotation by the same amount and with the same sense relative to the Hi template.

Assumptions
The Hi-based filamentary model of dust polarization assumes perfect alignment. We observe in Sec. A.3 that the Hi model displays no intrinsic parity violation. We, therefore, assume E Hi correlates withẼ d but not Illustration of magnetic misalignment. (Black ) Magnetic-field orientations derived from our Hi template for the same sky region and velocity bin as in Fig. 2 but downgraded to N side = 64. The pseudovector lengths are proportional to the template-implied polarization intensity. (Green) The same after applying a global misalignment angle ψ = 20 • , which is an unrealistically large amplitude for better visualization. (Greyscale colormap) The raw Hi intensity map (identical to the leftmost panel of Fig. 2). withB d , and we make a symmetric assumption for B Hi . Our assumptions are summarized by

Implications for cross-spectra
With Eqs. 15 and 16, we can derive the following relations between observable cross-spectra in terms of the misalignment angle ψ : and From the known positive T d B d and T d E d measured by Planck (Sec. 1.1), we expect ψ to be mostly positive in the range 100 500 on large sky areas away from the Galactic plane, e.g., with a 70% Galaxy mask (Clark et al. 2021). We also know that Planck Collaboration et al. 2016c) across the same multipole range and on the same sky area. Our qualitative expectations, then, are to find E Hi B d , T Hi B d and E d B d to be positive but to find B Hi E d to be negative.
We can make simple estimates of ψ with Eqs. 17, 18 and 19, though each is potentially biased by noise in the denominator: We could form a similar estimate from Eq. 20, but the E d B d measurement from Planck is especially noisy, so we ignore it for the remainder of this section. The four cross-spectrum ratios in Eq. 21 allow for tests of the misalignment ansatz without explicit calculation of ψ .
While positive T Hi B d and E Hi B d might be anticipated on account of the known positive nal to be entirely decoupled from the Hi-correlated components of the dust maps. In Sec. 5, we describe how to construct mock skies with exactly this property. These mock skies show positive T d B d but zero T Hi B d and zero E Hi B d . While we consider T Hi B d , E HI B d > 0 to be the most plausible expectation, it is formally nontrivial.
The negativity of B Hi E d is a new prediction of the misalignment ansatz. We can make quantitative predictions for this signal (Eq. 21), and comparisons are shown in Fig. 4 for our fiducial 70% Galaxy mask (Sec. 2.1). We mentioned above that we expect ψ to be smooth over the multipole range 100 500, so we expect B Hi E d to be smooth over similar multipoles. We can, therefore, gain in per-bandpower sensitivity by using the relatively large bin width of ∆ = 100.
In Fig. 4, we find that B Hi E d tends negative and is broadly consistent with the expectations of Eq. 21 over the full mask and in the northern and southern hemispheres independently. Due to the unavailability of suitable dust and Hi simulations, we do not attempt a statistical evaluation of the consistency. The plotted error bars are derived from Gaussian variances. As the dust field displays both non-Gaussianity and statistical anisotropy, these variances are meant only as a rough indication of the fidelity of the measurements.
We expect B Hi E d to be noisier than E Hi B d and T Hi B d , because r BHiB d is smaller (by roughly a factor of 2-3 for > 100) than r EHiE d and r THiT d , i.e., B Hi is a less accurate representation of B d than E Hi or T Hi is of E d or T d , respectively. The Hi-based polarization template is, therefore, more sensitive toẼ d modes mixed into the observed B d than toB d modes mixed into E d (Eq. 15).
We consider the results of Fig. 4 to be a first step in confirming that the misalignment ansatz of Eq. 15 is at least a partial description of the true sky. These re-sults avoided an explicit calculation of the misalignment angle ψ . In later sections, we will compute ψ directly.

MOCK SKIES
We construct a set of mock sky realizations in order to check for biases and spurious signals in the estimators that we will introduce in subsequent sections. We maintain the 2-point statistics of the true sky including correlations with the Hi-based polarization templates. These mock skies are phenomenological in the sense that they produce realistic observables without explicit appeal to the underlying ISM physics; in particular, these are not numerical ISM simulations.
Our mock skies include Gaussian noise, an Hi-based filamentary component and Gaussian dust. We arrange for all of the 2-point statistics to be the same as for the true sky, i.e., the mock skies replicate the measured X a Y b for X, Y ∈ {T, E, B} and a, b ∈ {d, Hi}. The Hi-based component is the same for all realizations and is derived from the true-sky Hessian template (Sec. 3).
In harmonic space, we express the mock-sky (S) map as a linear combination of an -filtered Hi template, a Gaussian dust component (G) and a Gaussian noise component (N ): for X ∈ {T, E, B}, where X Hi ( ) is the harmonicspace representation of the Hessian template (Sec. 3).
The -dependent coefficient in the Hi term is necessary, because D XHiXHi ∝ D XHiX d . To maintain D XHiX S = D XHiX d , we modify with the transfer function k (X) (Sec. A.2), which ensures consistency with the true Hi-Planck cross-spectra. While the Hi term is constant across realizations and based on the true sky, the Gaussian dust and noise are stochastic.
The power-spectra of the Gaussian dust and noise components are estimated from the measured dust and Hi power-spectra. We calculate these spectra after applying the 70% Galaxy mask. We compute X Hi ( ) from a masked map as well. As a result, the mock skies are well-defined only on the unmasked 70% of the celestial sphere.
Unlike the CMB, Galactic dust emission is statistically anisotropic, i.e., the statistics of the dust are different in different regions of sky. We approximate the non-stationarity by beginning with Gaussian noise and Gaussian dust that are isotropic and then modulating based on the statistical anisotropies in the Commander maps. The modulation is performed on scales much larger than those used in our analysis. The modulation field is smoothed to 14.7 • , twice the side length of a pixel with N side = 8, while most of our analysis is con-cerned with multipoles > 100, i.e., degree scales and smaller. We, therefore, expect negligible mode mixing.
A realization of the modulated dust mock sky is shown in Fig. 5 after highpass filtering to > 100, the multipole range targeted by most of our analysis. Before filtering, the mock skies are dominated by large-scale modes, which are bright and relatively poorly estimated, but those large-scale modes are irrelevant for most of our analysis. Visually, we find greater non-Gaussianity in the real T d map than in the mock sky. The polarization maps bear a greater resemblance to each other. A higher level of realism is unnecessary, since we use these mock skies only to check for biases in our estimators.
A breakdown of the mock-sky components is shown in Fig. 6, where we see that most of the dust power is in the Gaussian rather than the Hi component. We note that the unbiased power-spectrum estimator is crucial for avoiding a large noise bias in polarization, especially for 200. The Hi component shows negligible T B and EB but strong T E, which is arguably a defining characteristic of the filamentary polarization model. Although the Gaussian component dominates in T T , EE and BB, the Hi component accounts for roughly half of the T E signal.
We verified that the 2-point statistics of the mock skies are approximately equivalent to those of the true dust maps. In particular, we check that X S Y S (where "S" denotes a mock sky) matches X d Y d and that X Hi Y S matches X Hi Y d . The agreement is sufficient to test the estimators that we will introduce below.
We emphasize that our mock-sky framework is not intended to represent a null hypothesis for the purposes of statistical inference. In particular, the mock skies are missing much of the non-Gaussian structure in the true sky, even beyond the Hi-correlated component. Instead, because no aggregate misalignment has been input, these mock skies are useful for testing our estimators for spurious signals.

MISALIGNMENT ESTIMATOR
We present an estimator for the misalignment angle of a region of sky containing multiple pixels. In Clark et al. (2021), the angle difference between the dust and the Hi template was computed by The vertical-axis units are µK 2 RJ . The annotations indicate the spectrum type, e.g., T T for the upper-left plot. Blue shows the half-mission cross-spectrumD XY (Eq. 23). Orange shows the half-mission noise cross-spectrum. The Gaussian-dust spectrum is shown in green, and the Hi spectrum is shown in red. We note that the Hi-correlated component represents a minority contribution to all of the spectra except T E, which is a characteristic signature of a filamentary magnetically-aligned polarization model.
where c x ≡ cos (2θ x ) and s x ≡ sin (2θ x ). 4 While Eq. 24 measures the misalignment angle of a single pixel, care must be taken in computing the mean over multiple pixels, because ∆θ is a circular statistic. The values of ∆θ are restricted to [−90 • , 90 • ], but the endpoints of this range are physically identical. Naively averaging random values from this range will produce mean values that cluster near 0 instead of being uniformly distributed. As a result, noise fluctuations produce a multiplicative bias that suppresses the magnitude of ∆θ .
To account for the circularity of ∆θ, we use a modified version of the projected Rayleigh statistic (Jow et al. 2017), which is itself a form of alignment estimator (cf. Sec. 6.2 of Clark & Hensley 2019). The essence of the method is to consider terms of the form cos [2 (θ d (n) − θ Hi (n)) − ψ], where θ d (n) is the polarization angle measured by Planck, θ Hi (n) is the angle predicted by the Hi template and ψ is a free parameter independent ofn and representing the misalignment angle. We sum such terms over the selected map pixelsn and maximize with respect to ψ. Denote the maximizing value byψ. When θ d (n) − θ Hi (n) is random,ψ is also random. This is a plausibility argument thatψ is unbiased, but we will describe an explicit test below.
Rather than simply summing the cosine terms described above, we upweight pixels with higher signal-to-noise ratio in polarization. The weights are proportional to the product of the signal-to-noise ratios for Planck and Hi4PI. Denote the per-pixel weight by w(n).
We form the alignment metric where ψ is a free parameter and W = n w(n). The non-uniform weighting of the contributing pixels distinguishes our alignment metric from that of Sec. 6.2 of Clark & Hensley (2019), but it is an estimator for the same quantity. We maximize ξ(ψ) with respect to ψ and denote the maximizing value byψ. We can calculateψ analytically with the following prescription. Form where c x ≡ Q x /P x and s x ≡ U x /P x . Then we can express the alignment metric as from which the maximizing value can be found to bê where, because A(θ d = θ Hi ) > 0 and B(θ d = θ Hi ) = 0, this choice of arctangent ensuresψ = 0 in the case of perfect alignment. In the limit of a single pixeln, the estimatorψ is equivalent to Eq. 24, the ∆θ used in Clark et al. (2021). The added benefit ofψ is in aggregating pixels into patches without biasing the estimates low (as described at the beginning of this section).

Misalignment maps
We present maps ofψ in Fig. 7. We computeψ on masks defined by HEALPix pixels with various values of N side . For the lowpass-filtered ( < 702) results in the left column of Fig. 7, the misalignment angles are partially correlated between patches due to the presence of large-scale polarization modes. Part of the motivation for the bandpass filtering (101 < < 702) implemented for the right column of Fig. 7 is to remove these correlations and acquire approximately independent estimates in each patch. For N side = 32, the smallest patch size we consider, the side length of each mask is 1.8 • , which means that the above-mentioned bandpass filtering suppresses modes with wavelengths larger than a single patch. Most of the patch-to-patch correlations are removed by the bandpass filtering. (We will show in Sec. 7.2, however, that there is evidence for nontrivial spatial coherence ofψ that cannot be simply attributed to large-scale modes.) As the patch size decreases, regions of higher and lower variance emerge at all latitudes in a pattern that is similar to that of ∆θ in Fig. 3 of Clark et al. (2021). The above observations are broadly consistent between the lowpass-( < 702) and bandpass-filtered (101 < < 702) maps, but the former are visually smoother.

Test for estimator bias
To check for biases in our misalignment estimator, we measureψ on masks defined by HEALPix pixels as described above, artificially rotate the angles of the Planck polarization map by a known amount ψ 0 ∈ [−90 • , 90 • ] and then recomputeψ. We track the median of the distribution and find that it follows ψ 0 . We conclude that ψ is an unbiased estimator of misalignment angle.

Positive misalignment tendency
We observe a tendency toward positive misalignment angles in Fig. 7. To estimate the statistical significance, it is tempting to appeal to the central limit theorem. Unfortunately, the values ofψ(n), wheren here represents a particular patch, are neither completely independent nor identically distributed. By bandpass filtering to 101 < < 702 as in the right column of Fig. 7, we can achieve approximate independence of the estimates for differentn. We cannot, however, guarantee that the estimates are identically distributed.
Nevertheless, because the calculation is simple, we estimate a mean and standard error by appealing to the central limit theorem. For N side = 32, we findψ <702 = 1.9 ± 0.3 • , but we caution that the patches are nontrivially correlated with each other by the bright, large-scale polarization modes and, therefore, refrain from claiming any statistical significance. Restricting to the sky area allowed by our fiducial 70% Galaxy mask (cf. Fig. 11), we findψ <702 = 1.7 ± 0.3 • . After bandpass filtering, the patches are more (but not completely) independent, and we findψ 101< <702 = 0.9 ± 0.3 • , which implies a statistical significance of 3σ. Restricting to our fiducial 70% Galaxy mask, we findψ 101< <702 = 0.8 ± 0.4 • , which implies a significance of 2σ. We have deliberately limited ourselves to a single significant figure, because we consider these calculations to be crude estimates.

Relationship to dust properties
We note that the large-scale features ofψ are similar to those of the dust polarization fraction p d ≡ In the left column, the Planck and Hi maps have been lowpass filtered to < 702, which is approximately the Hi4PI beam scale. In the right column, the maps have been bandpass filtered to 101 < < 702, which is the multipole range used for much of the following analysis. The color scales are in degrees. We see greater variance in the northeast (top left) and southwest (bottom right), which are the regions with the lowest dust intensity, and we see a tendency for the misalignment angles to be positive.
For visual comparison, we provide in Fig. 8 a panel of maps displaying misalignment angleψ, dust intensity T d and dust polarization fraction p d . We omit an estimate of the correlation strengths, because the bright, large-scale modes induce covariances that are difficult to model. The low-column (small-T d ) sky regions in the northeast and southwest are also regions of increased variance inψ, as evidenced by the large fluctuations between neighboring patches, but a much clearer visual correspondence appears betweenψ and p d . The regions of larger p d show smaller variance inψ, and those regions also tend toψ > 0. For both T d and p d , the correspondence with regions of lower variance may be related to the signal-to-noise ratio of the misalignment measurement. In regions with higher polarized intensity, there is less variance inψ.
There may also be a connection to S d , the angle dispersion of the dust polarization, and to S Hi , that of the Hi polarization template. The former anticorrelates with p d , i.e., regions of greater polarization-angle coherence have larger polarization fractions (Planck Collaboration et al. 2020e). The variation in polarization-angle coherence may be related to the magnetic-field orientation relative to the line of sight (e.g., Hensley et al. 2019). The Hi-based dispersion S Hi and polarization fraction p Hi also anticorrelate, and the alignment of the dust polarization angle to the Hi template anticorrelates with S Hi (Clark & Hensley 2019). We would, therefore, expect that regions of larger polarization fraction correlate with regions of coherent magnetic misalignment, which is indeed what we observe.

Large sky areas
A major motivation for this study is to understand the origin of the parity-violating T B correlation measured by Planck on large fractions of the high-latitude sky (Clark et al. 2021). In addition to measuring the variation in misalignment angle across relatively small patches of sky (Fig. 7), we can apply our estimator (Eq. 29) to large sky areas and compare to expectations based on measured cross-spectra (Eq. 21). On large sky areas, the variation inψ is suppressed, and we can safely make a small-angle approximation. Then we expect Noise in the denominators may bias these expressions, but we are here looking only for a broad consistency and for the approximate level of aggregate misalignment on large sky areas. Sinceψ is estimated by reference to the Hi template, we expect greatest consistency with the dust-Hi cross-spectra, e.g., E Hi B d as opposed to the Planck -only T d B d and E d B d .
In Fig. 9, we compare the misalignment estimates from Eq. 30 for a 70% Galaxy mask (Sec. 2.1). We find a misalignment angle of ∼ 2 • that is coherent in the range 101 < < 702. As expected,ψ tends to be more consistent with the dust-Hi cross-spectra, especially E Hi B d and T Hi B d , which are more sensitive than B Hi E d (Sec. 4.2). The Planck -only T d B d is more discrepant (though not dramatically so) but reproduces the coherently positive behavior for < 500.
Theψ estimates are broadly consistent between the northern and southern Galactic hemispheres. In particular, the magnitude of the misalignment and its scale (multipole) independence are consistent. The Fig. 9) also shows a similar consistency between the hemispheres. The positiveψ, the approximate scale independence ofψ and the consistency ofψ between hemispheres are robust to the choice of Galaxy mask; we checked this for f sky ∈ {40%, 60%, 70%, 80%, 90%} (Sec. 2.1) and present some related results below in Sec. 6.7. The consistency between hemispheres begs for an explanation, which should be a target for future investigations.
The uncertainties onψ in Fig. 9 are derived from the scatter of our mock skies (Sec. 5) and are used only for visualization purposes, namely, to give a rough indication of the expected variance. We do not rely on these uncertainties for statistical inference.
Because the large-scale (low-) modes are difficult to reproduce in our mock-sky framekwork (Sec. 5), we have restricted Fig. 9 to 101 < < 702. We can calculateψ for < 101, but we cannot form a reliable uncertainty based on mock skies. Nevertheless, it is interesting to report the values forψ <101 . We find 1.4 • when using both hemispheres, 2.9 • in the northern hemisphere and 0.7 • in the southern hemisphere. The statistical weight cannot be evaluated in the present analysis, but it is noteworthy that the large scales show the same tendency for misalignment to be positive.

Aggregate global misalignment?
An intriguing possibility is that there is an aggregate global misalignment of ∼ 2 • . An aggregate misalignment would appear as an isotropic, multipoleindependent rotation of the dust polarization relative to the filamentary structures, i.e., ψ = const. The implied magnetic-field structure relative to the dust intensity field would be qualitatively similar to that depicted by the green pseudovectors in Fig. 3. A global polarization rotation can also be produced by a miscalibration We note the similarities in the large-scale features ofψ and p d .  Comparison of misalignment-angle estimates fromψ (Eq. 29) and from parity-violating crossspectra (Eq. 30) for a 70% Galaxy mask (Sec. 2.1) and for each hemisphere independently (bottom two rows). Thê ψ estimates are formed after filtering the Planck and Hi maps to each multipole bin (∆ = 100). The parity-violating spectra are shown with Gaussian uncertainties. The green and orange bands are nearly identical and difficult to separate visually. The error bars on theψ points are from the standard deviations of our mock skies (Sec. 5). We find an approximately scale-independent positiveψ that persists at the same level in both hemispheres and is broadly consistent with the spectrum-based estimators, especially those incorporating the Hi template.
of the absolute polarization angle or in the CMB by the phenomenon of cosmic birefringence (Sec. 1.3).
We consider miscalibration to be unlikely, since Planck estimates a systematic uncertainty of 0.28 • (Planck Collaboration et al. 2016d), nearly an order of magnitude smaller than our measured misalignment. In the following sections, we measureψ in small sky regions and search for correlated variations with other interesting quantities. The relative variation from region to region is insensitive to an overall miscalibration.
A global misalignment signal in the dust, which acts as a foreground for CMB measurements, would need to be accounted for in searches for cosmic birefringence, especially with methods that rely on the symmetry properties of the dust polarization, e.g., Minami et al. (2019).
As a consistency check, we modified our Hi template (Sec. 3) by imposing a global polarization rotation of 2 • . This rotation mixes E and B modes. Because the Hi template is dominated by E modes (E Hi E Hi /B Hi B Hi ∼ 5 for > 100), the effect is fractionally stronger in the B modes. We can estimate the expected impact of this modification by considering that 2 • ≈ 0.03 rad, so this should produce a percentlevel change in the correlations. We correlate with the Planck dust maps and find that the B-mode correlation for 100 < < 700 increases by 0.1-0.5% in addition to the original correlation of 10-25%, which is indeed a fractional increase of O(1%). We performed the same exercise with the opposite rotation, i.e., by −2 • , and we find an approximately symmetric decrease in the Hi-Planck correlation. These results are consistent with the estimates of Fig. 9 and increase our confidence in a true on-sky aggregate misalignment of approximately 2 • .
The T d B d -based estimate of ψ (red in Fig. 9) is coherent only over the range 100 500. Where it is nonzero, it also tends to be larger thanψ. The discrepancy may be an indication that T d B d and T d E d are affected by additional polarization sources that are missed by the filamentary misalignment model, and it is conceivable that the simultaneous positivity of T d B d andψ is merely coincidental. In Sec. 8, we seek further evidence of a relationship by analyzing small sky patches.

Varying sky fraction
We track the dependence of these misalignment estimates on the sky fraction f sky .
In Fig. 9, we considered only f sky = 70%, whereas we now allow f sky ∈ {20%, 40%, 60%, 70%, 80%, 90%, 97%} (Fig. 1). In Fig. 9, we considered multiple estimators for ψ . For simplicity, we now downselect to only two. One is Fig. 9), which uses only the dust field (cf. Fig. 9 of Clark et al. 2021). We contrast the dust-only estimator with one that includes some Hi filamentary information: ψ Hi×d ≡ atan2 (T Hi B d , T Hi E d ) /2 (green in Fig. 9). Whereas ψ d×d includes information from the entire dust field, ψ Hi×d collapses the misalignment estimate onto the filamentary modes. When the two are in agreement, the filamentary magnetic misalignment is representative of the parity violation in the full dust field. When they deviate from each other, the Hi template may be incomplete or inaccurate, or the full dust field may contain parityviolating contributions which are non-filamentary.
In Fig. 10, we show ψ d×d and ψ Hi×d for a variety of sky fractions f sky . We show the estimates for individual multipole bins (cf. Fig. 9), and we also highpass filter to form the broadband ψ >100 , which may potentially average away signal but is less noisy. We find that ψ Hi×d is consistently positive over all and remains in the range of 0-5 • , while ψ d×d is much more variable, especially at large f sky . We note that the two estimates display closer agreement at small f sky , i.e., when restricting to high Galactic latitudes. At the same time, we find that ψ Hi×d >100 steadily decreases from ∼ 3 • to ∼ 1 • as f sky increases from 20% to 97% (right plots of Fig. 10), a phenomenon which is observed in both hemispheres independently. The decline may be related to the fact that the Hi becomes a less robust tracer of dust at low Galactic latitudes where the column densities are relatively large (e.g., Lenz et al. 2017), so it may be that the Hi template is simply less representative of the dust field for large f sky .
Interestingly, all of the variations considered in Fig While random deviations constitute a form of "misalignment" relative to the Hi-defined filaments, it is unsurprising that such deviations are detected. The Hi-based polarization templates (presented here and in, e.g., Clark & Hensley 2019) correlate strongly with the Planck dust maps, but they are not identical, even within the limits of the Planck noise. If the term "magnetic misalignment" is to refer to any kind of deviation from the Hi template, then a detection of misalignment teaches us only that the Hi template is an incomplete description of the dust polarization field.
As a result of these considerations, we focus much of the rest of our analysis on a search for magnetic misalignment that displays certain types of coherence, which is less likely to be mimicked by random deviations from an Hi template. We search for coherence both in harmonic and map space. Harmonic coherence indicates that ψ is approximately constant with and manifests as a uniform rotation of the dust polarization pseudovectors relative to the Hi predictions (Fig. 3). We also refer to harmonic coherence as scale independence.
We restrict the analysis to the high-latitude sky by masking the Galactic plane to varying levels. Our fiducial choice is the Planck 70% Galaxy mask. For example, when dividing the sky into patches defined by pixels with N side = 8, we consider only those shown in Fig. 11. In subsequent sections, we will consider other, similarly-parameterized Galaxy masks and patch sizes.
In Sec. B, we introduce a number of cross-power and correlation metrics, which are used throughout the following sections.

Harmonic coherence
The relative constancy of the Planck T B/T E at high latitudes (see, e.g., Fig. 9 of Clark et al. 2021) is perhaps a hint that magnetic misalignment, if it is to address the mystery of the positive dust T B, ought to display a harmonic coherence in the multipole range 100 500. Indeed, we find that a direct calculation of the misalignment angleψ yields an apparent coherence over an even larger multipole range, tentatively across all < 702 ( Fig. 9 and Sec. 6.5). Is the apparent harmonic coherence an emergent phenomenon that appears only when aggregating large sky areas? In this section, we divide the sky into smaller patches and test whether harmonic coherence is a generic feature of magnetic misalignment.
To look for coherence in harmonic space, we bandpass filter the maps into two disjoint multipole ranges. For all of our results, 101 < < 702, so min = 101 and max = 702 can be taken as lower and upper limits, respectively, on the multipole ranges. We form a set of maps with < c − ∆ /2 and a set with > c + ∆ /2, where c is a transition multipole and ∆ is a multipole buffer between the two ranges. Misalignment angle ψ as a function of multipole in bins of width 100 for three example sky fractions f sky . We note that the "Dust only" estimates are more variable over and that the agreement with the "Dust × Hi" estimates tends to be better with smaller f sky , i.e., at higher latitudes. The 70% curves in the upper panel also appeared in Fig. 9 Highpass-filtered misalignment angle ψ >100 as a function of sky fraction f sky with each hemisphere also considered individually. We note that the "Dust × Hi" estimates are all fairly consistent with a steady decline with increasing f sky , while the "Dust only" estimates fluctuate strongly for f sky 70%.  The red area is allowed by the Galaxy mask, and the yellow patches are those that lie within the red area.
Letψ lo (n) be the misalignment angle estimated in the patch centered on sky coordinaten after filtering to < c −∆ /2, and letψ hi (n) be similarly defined after filtering to > c + ∆ /2. We will refer toψ lo (n) andψ hi (n), respectively, as the "lowpass-filtered" and "highpassfiltered" misalignment estimates. Recall, however, the multipole limits min and max mentioned above, so these estimates are, in fact, products of bandpass filtering.
Our correlation calculations must consider the circularity of the misalignment angle. We expectψ to cluster around 0 • mod 180 • . Even for the smallest masks of Fig. 7, which are defined by N side = 32, the majority of ψ values lie within [−45 • , 45 • ] mod 180 • . In our correlations, therefore, we ignore the circularity ofψ and instead force the values to their physical equivalents in the range [−90 • , 90 • ]. In a small minority of cases, we will miss correlations between angles that lie at opposite extremes of this range. In testing for a correlation, our choice is conservative. Since we will be using Spearman correlations, which operate on rank variables, a convenient ordering strategy is to form tanψ.
We form the Spearman cross-power (Spearman version of Eq. B3) where the sum is taken over patches labeled byn. Note that these Spearman cross-powers are not correlation coefficients, so the numerical values range outside of [−1, 1]. Correlation coefficients are less numerically stable in the presence of noise, so we prefer crosspowers for the purposes of establishing a relationship. In Fig. 12, we showψ lo ×ψ hi for several choices of c and ∆ for the patches of (2) hi Figure 12. Cross-power ss between misalignment angles estimated from lowpass-and highpass-filtered maps (Eq. 31): ψ lo ×ψ hi for the patches of Fig. 11. We consider three values for ∆ , and the data points have been offset from c for visual purposes. We also show the cross-power between halfmission splits in the shaded bands for lowpass-filtered (ψ (2) hi ). Thê ψ hi estimates become noise dominated for c 400. This is a test for harmonic coherence (scale independence) of magnetic misalignment, and we find a positive signal for c 450.
these measurements is mainly in the highpass-filtered misalignment estimatesψ hi .
We do not expect our mock skies (Sec. 5) to show a coherence over , because the Gaussian modes are resampled independently of each other and also independently of the Hi template. One concern might be that the masking creates mode correlations, and this was part of the motivation for introducing the multipole buffer ∆ . As ∆ increases, the two multipole passbands are further separated, and spurious correlations between the two are less likely.
As a null test, we calculateψ lo ×ψ hi for an ensemble of mock skies (Sec. 5), and we find the mean values to be consistent with zero. Recall that these mock skies are simplified in the sense that they are designed to reproduce only the 2-point statistics of the dust field, both in correlation with itself and with the Hi template. Nonetheless, they are helpful in checking that our estimators produce sensible results. The Hi template appears in these mock skies with the observed amplitude, and the only magnetic misalignment that has been input is due to random scatter. We see no positive bias in the mock-sky cross-powers. The positive signal seen in the real map (Fig. 12) must be due to a feature that it is absent in the mock skies.
Due to noise in theψ estimates, it is nontrivial to determine the fraction of the misalignment signal that is harmonically coherent. In computing a correlation coefficient, noise tends to dilute the true signal. Using half-mission splits as in Eq. B4 may lead to numerical pathologies when the denominators are small. From the decay ofψ (1) hi ×ψ (2) hi in Fig. 12, we see that the highpass-filtered estimateψ hi is especially noisy. The Spearman cross-powersψ (1) lo ×ψ (2) lo (red in Fig. 12) andψ (1) hi ×ψ (2) hi (purple) are limited only by noise and set rough upper limits on the cross-powerψ lo ×ψ hi . Even if the misalignment angles were perfectly coherent across multipoles, noise would suppress the cross-power. That ψ lo ×ψ hi is of the same order asψ (2) hi is an indication that, within the limits of the noise, the harmonically coherent component is contributing a nonnegligible fraction of the misalignment signal.
We estimate the statistical significance of the apparently positive signal in Fig. 12. For each choice of c and ∆ , we construct permutation tests (Sec. B.1), where covariances are preserved by using the same permutations for all choices. We combine the results with weights based on the half-mission cross-powers (bands in Fig. 12) and produce a single overall estimate of the statistical significance. For the case of Fig. 12, we estimate the statistical significance to be 2.3σ, where most of the sensitivity comes from the cross-powers with smaller c and smaller ∆ . This is becauseψ hi quickly becomes noise dominated as c increases, and increasing ∆ pushes the filter cutoff even higher. The data points in Fig. 12 are computed from the same maps but with different filtering parameters, so we expect them to be highly correlated. As such, combining the data points increases the overall significance only modestly. To a rough approximation, the overall significance can be estimated from the lowestc data point.
The results of Fig. 12 are based on the patches shown in Fig. 11, which are defined by N side = 8 and f sky = 70%. We can compute similar quantities for other values of N side and f sky , and the results are compiled in Tab. 1, where we see that the significances are generally between 2σ and 4σ for N side ∈ {2, 4, 8, 16, 32} and f sky ∈ [60%, 90%]. With smaller f sky , the significances tend to be smaller, but this may be simply a consequence of a decreased signal-to-noise ratio. On the full sky, the significances also tend to decrease, and this may be due to the inclusion of longer, denser sightlines at low Galactic latitudes. We do not attempt to combine the results from different choices of N side and f sky , because the covariances are difficult to capture.
We consider the results of Fig. 12 Fig. 12, for different values of N side (rows with side lengths provided parenthetically) and f sky (columns). All of the results are positive with little dependence on N side and f sky .
of misalignment angles measured in small regions of sky away from the Galactic plane.

Spatial coherence
We additionally search for spatial coherence of misalignment angles by considering neighboring pairs of sky masks. Although we bandpass filtered to 101 < < 702 in order to include only modes with wavelengths smaller than each patch, there is still a residual correlation between neighboring patches, which we detect with the mock skies of Sec. 5. To avoid the coherence due to common modes between neighboring patches, we again construct highpass-and lowpass-filtered maps as in Sec. 7.1. We correlate the lowpass-filtered estimate from each patch with the highpass-filtered estimate from each of its neighbors, and we simultaneously correlate with the opposite application of filters. The misalignment estimates that enter the correlation calculations are separated in both harmonic and map space.
We utilize the Spearman version of the 4-variable cross-power (Eq. B6) ψ lo (n) ×ψ hi (n ) ≡ S s tanψ lo (n), tanψ hi (n ), wheren andn are the central sky coordinates of neighboring patches. As in Sec. 7.1, we consider several choices for N side and f sky , where Fig. 11 shows one example. The sum is taken over all pairs of neighboring patches. Each patch appears multiple times in this sum, but each pair appears only once. Equation 32 measures a simultaneous correlation betweenψ lo (n) andψ hi (n ) and betweenψ hi (n) andψ lo (n ). The entire multipole range is being used in both patches but in two splits. Without multipole separation, we cannot pass a null test based on our mock skies (Sec. 5). With multipole separation, however, the mock skies show no significant cross-power. . Cross-power Ss (Eq. 32) between misalignment estimates from neighboring patches defined by N side = 16 and f sky = 70%. This cross-power splits the misalignment estimates in both harmonic and map space and is a test for spatial coherence. The plotting conventions are the same as in Fig. 12. We find a positive signal for c 500.
The results forψ lo (n) ×ψ hi (n ) (Eq. 32) are shown in Fig. 13 for patches defined by N side = 16. This patch area is 4 times smaller than that used for the measurement of harmonic coherence (Fig. 12). Measuring neighbor correlations at N side = 16 probes the spatial coherence within patches defined by N side = 8, so we are approximately measuring the spatial coherence within the patches of Fig. 12. For the particular example of Fig. 13, we find positive spatial coherence for c 500.
We estimate the statistical significance of the positive signal shown in Fig. 13 by following a prescription similar to that of Sec. 7.1. Combining all of the measurements in a manner that accounts for covariances, we estimate the statistical significance to be 3.6σ, where most of the sensitivity, as for harmonic coherence, comes from the lowc , low-∆ cross-powers.
We can compute similar quantities with other values of N side and f sky , and the results are compiled in Tab. 2, where we see that the spatial coherence tends to be stronger as the resolution is made finer. For N side ∈ {16, 32}, the significances are mostly between 2σ and 5σ. As the side length associated with N side = 32 is 1.8 • , these results may imply that magnetic misalignment displays a coherence length of O(1 • ).
We consider the results of Fig. 13 and Tab. 2 to represent tentative evidence for spatial coherence of misalignment angles. One perspective, however, is to consider the spatial coherence to be a necessary implication of the harmonic coherence that was established in  Table 2. Statistical significance (in units of σ) of measurements of spatial coherence ofψ (Eq. 32), e.g., those presented in Fig. 13, for different values of N side (rows with side lengths provided parenthetically) and f sky (columns). The significances are mostly positive and tend to increase with N side .
Sec. 7.1. We explained at the beginning of this section that the multipole split in Eq. 32 helps to evade correlations between neighboring patches, which appear even in our statistically aligned mock skies (Sec. 5).
The claim is thatψ lo (n) is correlated withψ lo (n ) even in the mock skies. So we chose to correlateψ lo (n) withψ hi (n ), and the mock skies show no correlation in this case. But the mock skies are also lacking harmonic coherence. In the real maps, harmonic coherence appears to correlateψ lo (n) withψ hi (n) andψ lo (n ) withψ hi (n ) (Sec. 7.1), but then we should expect, on the basis of the residual neighbor-to-neighbor correlations in the mock skies, a correlation betweenψ lo (n) andψ hi (n ) and betweenψ hi (n) andψ lo (n ). So there may indeed be a spatial coherence, but the crucial ingredient might be the harmonic coherence.

PARITY-VIOLATING CROSS-SPECTRA
We now investigate connections between the misalignment angleψ and the parity-violating cross- Secs. 4 and 6.5, we described the expected relationships. In Figs. 4 and 9, we tested the implications of these relationships on relatively large sky areas. In those particular cases, we used a 70% Galaxy mask and also checked the robustness of the results by restricting to the northern and southern hemisphere separately.
We now consider finer masks to search for coordinated variation in misalignment angle and parity-violating cross-spectra.

Random vs. harmonically coherent misalignment
As described in Sec. 7, we are interested in distinguishing between random and harmonically coherent misalignment. In both cases, we expect to find that misalignment angle is correlated with T B and EB, but we can impose additional constraints to isolate harmonically coherent correlations. Random misalignment is exemplified by the mock skies of Sec. 5. The mock skies show deviations from the Hi template, but the deviations are incoherent across multipoles, and there is no aggregate misalignment on large sky areas (cf. Sec. 6.5). We find the mock skies show significant correlations betweenψ and D T d B d and betweenψ and D E d B d . If we instead correlate between disjoint multipole bins, e.g.,ψ < c with D T d B d > c for some cutoff multipole c , we find that the correlations vanish. The real data, as we will show below in Secs. 8.3 and 8.4, display both types of correlations.
The mock skies display correlations betweenψ and D T d B d and betweenψ and D E d B d as a direct consequence of the known strong correlation between the Planck dust maps and the Hi templates (e.g., Clark & Hensley 2019). The Hi-dust correlation is maintained in the mock skies. The Hi component contributes nonnegligibly to the dust polarization, and the Planck maps can be viewed as perturbed versions of the Hi templates. From Fig. 6, we see that the perturbations need not be especially small; in fact, the Gaussian-dust component dominates over the Hi component, though only modestly. If the perturbations are random, the dust polarization angles are symmetrically distributed relative to the Hi template, andψ = 0. If, however, there is a region of sky in which the dust polarization angles are distributed asymmetrically relative to the Hi template, thenψ = 0; in this case, there will be a net chirality, which will in turn produce non-zero contributions to T d B d and E d B d . Even in the mock skies, there are regions of sky that fluctuate to non-zeroψ, and these regions tend to contribute non-zero T d B d and E d B d with a corresponding sign. Our estimators avoid noise biases, so the relevant fluctuations are likely due to on-sky dust components that deviate from the Hi template. This is the expected contribution of magnetic misalignment to the parity-violating dust polarization quantities, but we aim to investigate whether the observed T d B d and E d B d are consistent with random fluctuations away from the filament orientations -as exemplified by the mock skies -or show evidence for harmonic or spatial coherence, which might be expected from a physical misalignment between the magnetic field and dusty filaments.
It is important to note that, whileψ > 0 implies a tendency toward T d B d , E d B d > 0, the converse is not guaranteed. It is possible to have T d B d , E d B d > 0 but ψ = 0. Our mock skies (Sec. 5) illustrate this point. They are constructed to retain the T d B d spectrum of the true dust maps, but this property is placed entirely in the Gaussian component, which is statistically independent of the Hi component. As such, the mock skies display no aggregate misalignment (beyond realization-dependent scatter). For example, on a 70% Galaxy mask, the ensemble mean ofψ is zero, although the mean T S B S spectrum is positive for 100 500 as for the true T d B d (Sec. 1.1).
To the extent that the observed T d B d is related to magnetic misalignment, an outstanding question is whether the real dust T B is a consequence of physical misalignment or random scatter. Thus, we search for harmonically coherent relationships between misalignment angle and parity-violating cross-spectra. The mock skies will help us to make the distinction, since harmonic coherence is not included in them.

Misalignment controls T B
We begin the investigation with sky areas that are only modestly smaller than those of Secs. 4 and 6.5. In this limit, the aggregate misalignment angles are small, and the expected relationship to parity-violating crossspectra can be approximated as (cf. Eq. 30) We will focus more on the dust-only spectra T d B d and E d B d as opposed to the dust-Hi spectra T Hi B d , E Hi B d and B Hi E d , but similar operations can be performed for either set. We divide the sky into patches defined by N side = 2 within an overall 70% Galaxy mask. In each patch, we measureψ as in Fig. 7 (with an additionally imposed Galaxy mask). We form a combined mask from the patches withψ larger than the median value, and we form an analogous mask for the patches withψ smaller than the median; the masks are shown in the lower right of Fig. 14, where the inputψ values are from maps that have been filtered to 101 < < 702 (as in the right column of Fig. 7). We repeat for maps that have been lowpass filtered to 101 < < c (which will be labeled by the subscript "lo") and for maps that have been highpass filtered to c < < 702 (subscript "hi") for c ∈ {202, 302, 402}.
We calculate auto-and cross-spectra for the full combination of patches, for the large-ψ samples and for the small-ψ samples. To avoid sharp mask features at shared vertices of the HEALPix-defined patches, we use a conservative apodization scale of 5 • for the results in this subsection. Note that, due to the HEALPix pixelization and the increased apodization scale, the full combination of patches represents a smaller overall sky area than produced by the fiducial 70% Galaxy mask of Secs. 4 and 6.5. The reduction is somewhat severe and leaves a sky fraction of only 37%.
We convert the auto-and cross-spectra to misalignment estimates according to Eq. 33. The results are shown in the panel of spectra in Fig. 14.
We find that the T d B d -and E d B d -based misalignment estimates increase and decrease in a manner consistent with theψ-based mask definition. The large-ψ masks tend to produce larger , though the latter is much noisier. The smallψ masks tend to produce smaller spectrum-based estimates; interestingly, the resulting T d B d /T d E d (blue in the top row of Fig. 14) is broadly consistent with zero rather than negative. This may suggest that the positive T d B d measured at high Galactic latitudes is due to a few regions of sky with positive misalignment and that the rest of the sky respects parity.
Whenψ is estimated over the same multipole range as the spectra, as in the leftmost column of Fig. 14, we cannot distinguish between the case of random fluctuations and that of harmonically coherent misalignment (Sec. 8.1). To isolate the harmonically coherent signal, we compare estimates from disjoint multipole ranges as in Secs. 7.1 and 7.2. The three rightmost columns of the panel in Fig. 14 show the results when estimatingψ from restricted multipole ranges. Because the dust is brighter at low multipoles, the selections based onψ lo tend to be similar to those based on the unfilteredψ. For T d B d and E d B d to respond to theψ selection in the disjoint multipole range is an indication of the harmonic coherence of magnetic misalignment, which was demonstrated in Sec. 7.1 but has now been explicitly connected to T d B d . With multipole splits, the data are too noisy to make a confident claim about E d B d . The T d B d results, which are less vulnerable to noise fluctuations, are similar for all choices of multipole filtering; furthermore, the small-ψ lo and small-ψ hi results tend to track each other, as do the large-ψ lo and large-ψ hi results. This may be yet another indication of the harmonic coherence ofψ, i.e., theψ-based patch selections are broadly similar in all multipole ranges.
We estimate the statistical significance of the multipole-split results by making random patch selections to define the masks. We preserve covariances by using the same randomization for all c . In combining the results from all c , we weight by a measure of signal-to-noise ratio (cf. Sec. 7.1), and we estimate the overall significance of the harmonic coherence to be 2.2σ for T d B d and 0.8σ for E d B d .

Correlations with misalignment angle
We search for correlations between misalignment angleψ and the parity-violating cross-spectra In general, we see an increase in the spectrum-based misalignment estimates for large-ψ masks and a decrease for small-ψ. Black points are for the total patch collection (combined red and blue in the lower-right map) and are the same across each row. Red points are from large-ψ selections; blue from small-ψ. The leftmost column shows the results from filtering to 101 < < 702. In the remaining columns, the patch selections are based onψ lo , computed after lowpass filtering to 101 < < c (darker colors and downward markers), orψ hi , computed after highpass filtering to c < < 702 (lighter, upward), and the background shaded regions indicate the multipole splits. The brown band shows the binned difference between the large-and small-ψ lo estimates; yellow shows the same forψ hi . (Lower right) Patch selection based onψ after filtering to 101 < < 702. Large-ψ patches are in red; small-ψ in blue. The underlying white region is our fiducial 70% Galaxy mask. tions in Eqs. 17, 18, 19 and 20. For small angles,ψ is expected to track the parity-violating cross-spectra, but there are additional scaling factors. Rather than correlating the spectra withψ directly, we transformψ according to Eqs. 17,18,19 and 20. We compute Spearman half-mission correlation coefficients (cf. Eq. B5), because both variables entering each of the following calculations are derived from the same Planck dust modes and, therefore, subject to covariant noise fluctuations. Due to the aforementioned transformations, the variables are different for each correlation calculation. We compute (cf. Eqs. 17,18,19 and 20) r (HM) for x ∈ {d, Hi} and where Eq. 35 is expected to be negative and all others positive. We show these correlations and maps of the parity-violating quantities in Fig. 15 for patches defined by N side = 8 (as in Fig. 11). The correlations have the expected sign in all cases. Half-mission cross correlations (Eq. B5) avoid noise covariance but not sample variance in the dust measurements. The map features that produce positiveψ also produce positive T d B d , E d B d , T Hi B d and E Hi B d and negative B Hi E d (Sec. 8.1). While the effective mode weighting in calculatingψ is different than in the crossspectra, we nevertheless find a correlation between the two in our mock skies (Sec. 5), for which the non-Hi component is statistically independent of the Hi component. In particular, the mock skies approximately reproduce the results of Fig. 15.
That our mock skies show correlations between misalignment angle and parity-violating cross-spectra is an indication that the correlations of Fig. 15 could be attributed to random fluctuations away from the Hi template (Sec. 8.1). This is yet another motivation to restrict the search to signals that are coherent in either harmonic or map space rather than correlating identical patches with identical multipole bins.
Independent of the distinction between random and harmonically coherent misalignment, the correlations of Fig. 15 disfavor the presence of significant confounding contributions to parity violation in the polarization field. A priori, we might have expected nonfilamentary contributions to dilute the relationship be-  Figure 15. Maps of misalignment angleψ and the five parity-violating cross-spectra binned to 101 < < 702 for patches defined by N side = 8 and f sky = 70%. Forψ, the units are degrees; for T d B d and E d B d , µK 2 RJ ; and, for THiB d , EHiB d and BHiE d , µKRJ K km/s. For the cross-spectra, the correlation with the transformedψ (Eqs. 36, 37, 34 and 35) is given in the subtitle with a 1σ uncertainty (Sec. B.1). We find all of the correlations predicted by the misalignment ansatz (Eqs. 17,18,19 and 20).
tween Hi-based misalignment angle and parity-violating cross-spectra, especially when we consider that the Hicorrelated component is a minority contributor to the dust field (Fig. 6). Results like those of Fig. 15 and the leftmost panel of Fig. 14 suggest that the Hi template is sufficiently significant and representative to provide a reference for searches for parity violation.

Harmonic coherence of parity violation
Instead of directly correlatingψ with, e.g., D T d B d , we define disjoint multipole ranges and correlate the lowpass-filtered quantities with the highpass-filtered. This is similar to the multipole splits described in Sec. 7.1 and better extracts a signal that is coherent across multipoles. For the misalignment angle, we lowpass or highpass filter the map to formψ lo orψ hi , respectively. For the spectra, we simply bin the lower or higher multipoles to form D XY lo or D XY hi , respectively. As in Sec. 7.1, the "lowpass-filtered" multipole range is 101 < < c − ∆ /2, and the "highpass-filtered" is c + ∆ /2 < < 702, where c is a transition multipole and ∆ is a multipole buffer between the two ranges.
We now modify Eqs. 34, 35 36 and 37 to correlate across multipole splits. We filter the spectra in the same way and the misalignment angle in the opposite way, e.g., we correlate and our hypothesis is that the connection is provided byψ hi , even though the latter is estimated in a disjoint multipole range. We look for a simultaneous correlation when the multipole ranges are switched.
For this purpose, we use the Spearman version of the 4-variable cross-power (Eq. B6) and similar combinations for the other four parityviolating cross-spectra. The cross-powers are presented in Fig. 16 for patches defined by N side = 8 (as in Fig. 11). We also form these cross-powers for an ensemble of mock skies. As noted earlier, our mock skies cannot be used for null-hypothesis testing, but they are useful for testing basic properties of our correlation metrics and misalignment estimators. We find that the mock skies produce null results within the realization-dependent scatter.
To assess the noise level in the parity-violating cross-spectra, we consider the half-mission crosspowers (Spearman version of Eq. B3) where X (i) and Y (i) are the X and Y fields, respectively, from the ith half mission. We form a similar quantity for the highpass-filtered observables: Both XY (1) lo ×XY (2) lo and XY (1) hi ×XY (2) hi are limited only by noise. When noise is subdominant to the sky components, they will show strong positive signals; when noise is significant, they will decay to zero. We plot XY hi as red and purple bands, respectively, in Fig. 16, where we find that, in general, the high-quantities are substantially noisier than the low-quantities.
As discussed in Sec. 7.1, half-mission cross-powers likê ψ (2) hi set rough upper limits on the observable strength of the signals we are seeking. In the case of Fig. 16, we must consider the fidelity of botĥ ψ (red/purple bands in Fig. 12) and XY (red/purple in Fig. 16). Whenψ × XY is of the same order as the halfmission cross-powers, a non-negligible fraction of the variation in XY is associated with harmonically coherent misalignment. In Fig. 16, this is the case for E Hi B d , hi , are too noisy to make a reliable comparison.
We estimate the statistical significance of measurements like those of Fig. 16 by following a prescription similar to those of Secs. 7.1 and 7.2. The weights used in combining the measurements now account for noise in bothψ (bands in Fig. 12) and the cross-spectra (bands in Fig. 16). For the example of Fig. 16 hi cross-powers (Eqs. 39 and 40 shown as red and purple bands in Fig. 16) are 2-3 times smaller than the corresponding T d B d quantities at low c and have fractionally larger uncertainties. Furthermore, 2) 8 (7.3 • ) 1.9 (−0.9) 1.5 (−1.5) 2.9 (−0.8) 3.7 (−0.9) 3.9 (0.8) 4.4 (0.2) Table 3. Statistical significance (in units of σ) of measurements of harmonically coherent correlations betweenψ and T d B d (E d B d ), e.g., those shown in Fig. 16, for different values of N side (rows with side lengths provided parenthetically) and f sky (columns). Theψ × T d B d significances are all positive and tend to increase with N side and f sky . Theψ × E d B d significances are mostly positive but are generally smaller.
The real data, within the limits of the noise, are broadly consistent with our expectations, namely, positive correlations betweenψ and T d B d , E d B d , T Hi B d and E Hi B d and negative betweenψ and B Hi E d , though these signals disappear for some choices of N side , c and ∆ . In particular, the signal tends to decay as c increases, which we expect due to increased noise in the highpass-filtered quantities. For N side = 4, the expected negative correlation with B Hi E d appears only for c 350 and is fairly weak. For N side = 8, the expected correlation with E d B d disappears.
These results build confidence in our picture of harmonically coherent magnetic misalignment. Alternatively, these correlations can be considered necessary implications of the harmonic coherence ofψ (Sec. 7.1) coupled with the expected relationship betweenψ and the parity-violating cross-spectra (Secs. 8.1 and 8.3). With this perspective, the cross-powers in Fig. 16 are merely tests for consistency.
We draw special attention to the correlations betweenψ and T d B d , as the positive T d B d measured by Planck has been recently discussed in the literature (Huffenberger et al. 2020;Weiland et al. 2020;Clark et al. 2021;Huang 2022). Although correlations betweenψ and E d B d yielded only weak results, we note the positivity of (E d B d ) (2) lo (Eq. 39 shown as the red band in Fig. 16), which indicates on-sky variation in E d B d that rises above the noise level. Variation in E d B d may be attributable to sample-variance fluctuations of underlying parity-even statistical processes, but the particular dust realization that we observe is a foreground that must be mitigated for, e.g., measurements of the CMB. Spatial variation in E d B d may be of relevance for measurements of cosmic birefringence (Minami et al. 2019;Minami & Komatsu 2020;Diego-Palazuelos et al. 2022;Eskilt & Komatsu 2022). If future measurements can more confidently establish a relationship betweenψ and E d B d , foreground removal could be performed more robustly by, e.g., relating E d B d to T d B d and other observables (e.g., Eqs. 17,18,19 and 20).

CONCLUSION AND OUTLOOK
We have extended the work of Clark et al. (2021) in establishing a connection between dust T B correlations and the magnetic misalignment of interstellar dust filaments. We have introduced a new version of a Hessianbased Hi polarization template, which correlates more strongly with dust B modes than the RHT-based template used previously (Sec. 3). We introduced several spectrum-based misalignment estimators formed from the auto-and cross-spectra of Planck dust maps and Hi polarization templates (Sec. 4.2), and we also introduced a map-based estimator for misalignment angle (Sec. 6). We have presented maps of the misalignment angle (Sec. 6.1), which show a tendency to positive values and a visual correlation with dust polarization fraction. We have provided evidence for the scale independence (harmonic coherence) of the misalignment angle for multipoles 700 (Secs. 6.5 and 7.1) and for spatial coherence on angular scales of ∼ 1 • (Sec. 7.2). On large sky areas at high Galactic latitudes, we find a scale-independent misalignment angle of ∼ 2 • , which is robust to a variety of masking choices (Sec. 6.7). We have described a set of mock skies (Sec. 5) containing Hibased filamentary structure as well as Gaussian-random components, and we have used these mock skies to refine our notion of magnetic misalignment. In particular, we have explored the question of whether the measured misalignment between Hi filaments and the magnetic-field orientation is consistent with random fluctuations in the polarization field (Sec. 8.1). This question motivated a search for scale independence (harmonic coherence) as a salient physical property of magnetic misalignment. We find evidence for a scale-independent correlation between misalignment angle and dust T B (Sec. 8.4). With the noisier EB, we find a correlation for some but not all of our masking choices. We also find that the observed positive dust T B may be due to a few regions with strong positive misalignment while the rest of the sky largely respects parity (Sec. 8.2).
In general, the picture that is beginning to emerge contains the following features: • On large scales at high Galactic latitudes, there is a global tendency toward an aggregate misalignment of ∼ 2 • (Secs. 6.5, 6.6 and 6.7).
• Magnetic misalignment is a reliable predictor of parity violation in the dust polarization (Secs. 8.2 and 8.3).
We now provide suggestions for potential improvements to our analysis.
1. The ansatz (Sec. 4) could be modified to allow for only a fraction of the dust to participate in misalignment. In this work, it is assumed that all of the dust is misaligned, but this may dilute the sensitivity of our estimators. In Clark et al. (2021), this type of concern was addressed in estimating the misalignment-induced EB in Eq. 12.
2. The Hi template (Sec. 3) could be improved to correlate more strongly with the Planck dust maps.
In this work, we introduced a new Hessian-based template, which correlates more strongly with the dust B modes than the RHT-based template used in Clark et al. (2021), but the correlation is still less than 20% for 200. While the Hi-based filamentary model may be fundamentally limited due to diffuse non-filamentary dust or other dust morphologies, we consider it more likely that a dedicated exploration will yield stronger correlations with the measured dust polarization (Halal et al. in prep.). Magnetic misalignment is a perturbation to the filamentary model, so an increased correlation would improve the sensitivity of all of the Hi-related estimators presented in this work.
3. More realistic mock skies or simulations (Sec. 5) will aid in the physical interpretation of our estimators. For example, the magnetohydrodynamic (MHD) simulations of Kim & Ostriker (2017), which can be converted to dust polarization maps (Kim et al. 2019) that were considered in Clark et al. (2021), model the Solar neighborhood and are publicly available at a resolution of N side = 128. Similar simulations with higher resolution and synthetic Hi observations could be analyzed with our estimators. Misalignment could also be investigated in synthetic dust polarization observations directly by searching for, e.g., scale independence in T B and EB. In Sec. 4 of Clark et al. (2021), this type of analysis was performed on a limited multipole range (60 < < 120).
Higher-resolution simulations will enable an extension to higher multipoles and further investigation of the link to underlying physics.
4. The pixels weights w(n) (Sec. 6) that enter the calculation ofψ are likely suboptimal. We checked that our choice reduces variance relative to a uniform weighting, but we have not explored the full space of possibilities. A better choice may be a Wiener filter that prevents a few bright pixels from dominating. Similarly, the correlation metrics used in Secs. 7 and 8 could be defined with weights to suppress noisy regions of sky.
5. The large-scale (low-) misalignment should be considered more rigorously, because dust polarization is dominated by these modes. We have mostly limited our investigation to 100 to avoid large-scale covariances. But there appears to be a strong positive misalignment on large scales (Fig. 7), and we speculate that this may be related to the magnetic-field structure in the vicinity of the Local Bubble (e.g., Lallement et al. 2003;Alves et al. 2018;Leike et al. 2020;Pelgrims et al. 2020;Vergely et al. 2022).
6. Other sources of parity violation should be considered, since magnetic misalignment alone may be insufficient to account for, e.g., the observed T B. We mention in Sec. A.3 that the distribution of dust filaments may itself display a chiral asymmetry even in the limit of perfect magnetic alignment. Both the Hessian-based and RHT-based Hi templates, which assume perfect alignment, predict a rise in EB for 100 ( Fig. 18), though the expected signal is below the Planck noise levels. We defer the investigation of this morphological parity violation to future work. 7. Other magnetic-field tracers such as starlight polarization and Faraday rotation could be incorporated to better understand the threedimensional manifestation of magnetic misalignment. With stellar distance measurements from, e.g., Gaia (Gaia Collaboration et al. 2016), starlight polarization measurements can enable a tomographic reconstruction of Galactic magnetic fields, though this technique is sensitive only to the plane-of-sky component (Panopoulou et al. 2019). Faraday rotation measures probe the line-of-sight magnetic-field component (Hutschenreuter & Enßlin 2020) and can be combined with model expectations or plane-of-sky observations to constrain the three-dimensional magnetic-field structure (e.g., Tahani et al. 2019Tahani et al. , 2022. Our misalignment analysis can be applied to a variety of ISM environments. As a method of studying the relative orientations of magnetic fields (not necessarily with dust polarization) and density structures (not necessarily with Hi), our approach is complementary to those of, e.g., Planck Collaboration et al. (2016a), Soler et al. (2017) and Fissel et al. (2019), which consider both the diffuse ISM and molecular clouds.
The study of parity violation in Galactic dust polarization is of central importance both for cosmology and for ISM physics. Our investigation has been limited by noise in the Planck polarization maps, and we, therefore, recommend follow-up surveys at millimeter and submillimeter wavelengths covering large sky fractions with resolution similar to or finer than the Hi4PI beam width (16.2 ). More sensitive measurements will become available from upcoming projects including the spacebased LiteBIRD (Hazumi et al. 2020)   Here we provide additional information related to the Hessian method (Sec. 3) that supports some of the analysis choices made in this work. In the comparisons below, we will occasionally use the subscript "H" to refer to the Hessian method and "RHT" for the Rolling Hough Transform. With the subscript "Hi", we implicitly refer to the Hessian method as in the main text.

A.1. Comparison with the Rolling Hough Transform
The Rolling Hough Transform (RHT) is another filament-finding algorithm (Clark et al. 2014) from which we can produce polarization templates that correlate strongly with the dust polarization measured by Planck (Clark et al. 2015;Clark & Hensley 2019).
We find that the Hessian correlates more strongly with Planck dust B modes than the RHT for 100, and we show a comparison in Fig. 17. In E modes, the two perform similarly. In these comparisons, the RHT is constructed from velocities spanning −90 to 90 km/s (as in Clark & Hensley 2019), while the Hessian template is constructed from the restricted range of −15 to 4 km/s, (Sec. 3.1).  Figure 17. Correlation r between Planck dust maps (d) and the RHT, the Hessian method (H) and the raw Hi4PI intensity map (IHi4PI) with multipole bin width ∆ = 25 on a 70% Galaxy mask.
For the T template, we consider two choices for the RHT. One might use I Hi4PI , i.e., the Hi4PI intensity map (HI4PI Collaboration et al. 2016) without any processing. Since we are especially interested in the fil-amentary component of the dust intensity, it may be preferable to highpass filter the Hi4PI intensity as in the first step in the RHT algorithm of Clark & Hensley (2019). The filter is implemented as an unsharp mask with FWHM = 30 . Denote the highpass-filtered intensity by T RHT . In Fig. 17, we find that I Hi4PI correlates more strongly than T RHT with Planck T modes. This is not necessarily the relevant metric, however, since we are specifically targeting filaments. We also present in Fig. 17 the correlation with Planck E modes, since the T E correlation is a signature of filamentary polarization. We find that r TRHTE d is generally larger than r IHi4PIE d . For this reason, we prefer T RHT as a template for filamentary dust intensity.
The Hessian intensity T H (Eqs. 11 and 14) is defined mainly by the Hessian eigenvalues. In Fig. 17, we find that r T H T d is smaller than r TRHTT d but that r T H E d is similar to and, at high , slightly larger than r TRHTE d .
On account of the greater correlation in B modes, we have selected the Hessian-based template as our baseline. The T templates considered above all correlate strongly with both T d and E d and at roughly the same level. The E templates for both algorithms correlate with E d at roughly the same level.
We defer to future work a more detailed investigation of these and related filament-finding alogirthms for the construction of polarization templates (Halal et al. in prep.). Each can be modified and tuned by making different choices for, e.g., velocity binning, weighting, spatial filtering, etc.

A.2. Transfer function
The Hi-based polarization templates have different mode structures than the Planck dust maps. For example, the Hessian method upweights small-scale features; the E H and B H power spectra increase with . The RHT also upweights small-scale features but especially emphasizes the multipole range 300 500 for E RHT and 150 350 for B RHT (with RHT parameters set to those of Clark & Hensley 2019).
Correlations, e.g., those presented in Clark & Hensley (2019), are insensitive to differences in mode structure, because they are evaluated in individual multipole bins. The upweighting of one multipole bin relative to another is normalized out of the calculation. The difference in mode structure can be viewed as a multipole-dependent transfer function.
For the purposes of converting our Hi-based polarization templates into quantities that are directly comparable to the Planck dust maps, we assume a transfer function that depends only on multipole : This transfer function converts an Hi-based quantity into dust-intensity units with a mode structure that is directly comparable to the observed dust field. We find that k (X) rises strongly at low multipoles, which is an indication that the Hi templates tend to under predict the amplitude of large-scale dust polarization relative to small-scale. In spite of this underprediction, the correlations, which normalize out the dependence, are actually stronger at low . There is no guarantee that the transfer function k (X) provides a representative estimate of the amplitude of Hi-based modes in the real dust maps. The amplitudes may depend on both and m, the spherical-harmonic eigenvalues. We use k (X) as a rough conversion factor to make direct comparisons between the Hi templates and the real dust maps.
Because the Hessian-based template correlates nonvanishingly with Planck up to at least = 750, the transfer function remains usable across the entire multipole range considered in our analysis.

A.3. Parity in the templates
The Hi templates are produced under the assumption of perfect magnetic alignment. Even so, chirality in the filament morphology could produce parity-violating signatures such as non-zero T Hi B Hi and E Hi B Hi .
We computed these parity-violating cross-spectra for both the Hessian and the RHT templates. 5 To determine if the results are significant, we compare to the T d B d and E d B d spectra from the Planck dust maps. To make this comparison, we applied the transfer function k (X) introduced in Sec. A.2, which converts the Hi templates into dust-intensity units. The results are shown in Fig. 18. The Planck T d B d displays a positive signal, and T Hi B Hi is negligible in comparison. The Planck E d B d appears to be consistent with noise, and we find that E Hi B Hi is negligible in comparison with the fluctuations. The above comparisons are restricted to > 100, which is the target multipole range of the analysis presented in this work. Based on these observations, the intrinsic (or morphological ) Hi-based T Hi B Hi and E Hi B Hi are assumed in our analysis to vanish.
Intriguingly, however, the Hi-based E Hi B Hi shows a rise for < 150. With finer multipole binning ∆ , we find that this rise persists down to = 17 with ∆ = 10, the lowest bin center with the finest binning that we checked. There is a corresponding rise in T Hi B Hi that persists down to = 27. In all cases, the Hi-based predictions are subdominant to the expected noise in the Planck measurements but only by a factor of ∼ 10. We defer to future work an investigation of these low-Hi-based predictions, which could represent a source of parity violation independent of magnetic misalignment.
We make use of a variety of correlation metrics. We wish for these metrics to be numerically stable, unbiased by noise or other covariances and, in some cases, sensitive to two different effects simultaneously.
Let X be the sample mean for a set of n measurements {X i }. The mean-subtracted observable is The index i will be labeling the central sky coordinate of small patches. Unless the sky mask has been chosen to retain relatively isotropic dust power, the dust intensity can vary dramatically across the observing region. It is, therefore, likely that the set of observables {x i } is dominated by the brightest patches. Because many of the quantities of interest range over several orders of magnitude, we prefer metrics related to the Spearman rank correlation coefficient, for which the observables are converted to rank variables. This avoids overweighting bright sightlines. By collapsing the observables onto ranks, the absolute magnitudes are less important, and both large and small values of X contribute equally. We consider correlation metrics for two data vectors X and Y . When the data vectors are noisy, it is more useful to consider a cross-power, which is essentially the numerator of a Pearson correlation coefficient. We define the Pearson cross-power of X and Y as half-mission cross-correlation coefficient r (HM) p (X, Y ) ≡ s p X (1) , Y (2) + s p X (2) , Y (1) 2 s p (X (1) , X (2) ) s p (Y (1) , Y (2) ) , (B5) which measures a simultaneous correlation between opposite data splits of X and Y . The convenient splits for Planck are half-mission (HM) maps. As above, we form the Spearman version simply by using s s in place of s p .
In some cases, we will wish to measure simultaneous correlations between two pairs of observables, e.g., a correlation between W and X and a correlation between Y and Z. We define the 4-variable cross-power As above, we can construct a Spearman version of this cross-power, which we denote S s (W, X, Y, Z), by converting to rank variables.

B.1. Statistical inference
For statistical inference regarding our correlation metrics, we use permutation tests. For 2-variable metrics, we randomly permute Y to obtain π(Y ), where the function π defines a random permutation that here acts on the data vector Y . We then compute, e.g., s s (X, π(Y )). With a large number of permutations, we can build a null-hypothesis distribution for s s (X, Y ). Similar permutation tests can be formed for the other 2-variable correlation metrics. In this work, each ensemble of permutations contains 200 realizations. Our results change only negligibly with larger permutation ensembles.
For the 4-variable correlation metrics, we also appeal to permutation tests, but we coordinate the permutations of X and Z. We form many realizations of, e.g., S s (W, π(X), Y, π(Z)).
We will often estimate uncertainties on our correlation coefficients by converting to z scores and taking half the difference between the value at 1σ and the value at −1σ. For s s (X, Y ), call this uncertainty estimate σ[s s (X, Y )], and we maintain similar notational conventions for the other correlation metrics.