Mapping the Magnetic Interstellar Medium in Three Dimensions over the Full Sky with Neutral Hydrogen

Neutral hydrogen is observed via the λ21 cm spin-flip transition, and is a ubiquitous tracer of interstellar gas: no known sightline lacks Galactic H i emission. Recent technological advances have enabled high dynamic range observations of the H i sky. The Leiden/Argentine/Bonn Survey (LAB; Kalberla et al. 2005), long the gold standard full-sky H i survey with angular resolution ${\vartheta }_{\mathrm{fwhm}}=36^{\prime}$, was recently superseded in both sensitivity and resolution by the H i 4π Survey (H i4PI; HI4PI Collaboration et al. 2016). H i4PI, with angular resolution ${\vartheta }_{\mathrm{fwhm}}=16\buildrel{\,\prime}\over{.} 2$ and spectral resolution δv = 1.49 $\mathrm{km}\,{{\rm{s}}}^{-1}$, is the highest-resolution full-sky H i survey to date, achieved by combining the Effelsberg–Bonn H i Survey of the northern sky (EBHIS; Winkel et al. 2016) with the Parkes Galactic All-Sky Survey in the south (GASS; McClure-Griffiths et al. 2009). The Galactic Arecibo L-Band Feed Array H i Survey (GALFA-H i; Peek et al. 2018) is a higher-resolution (${\vartheta }_{\mathrm{fwhm}}=4\buildrel{\,\prime}\over{.} 1$, δv = 0.184 $\mathrm{km}\,{{\rm{s}}}^{-1}$) survey of the entire sky visible to the 305 m Arecibo telescope, ∼32% of the celestial sphere. For a comprehensive comparison of modern large-area H i surveys, see Winkel et al. (2016).

In this work, we make use of the alignment between slender H i features in narrow spectral channels and the plane-of-sky magnetic field (Clark et al. 2015). Analysis of both GASS and GALFA-H i data found that the alignment between H i features and the projected magnetic field is sharper for higher spatial resolution H i data (Clark et al. 2014). For our purposes, this amounts to a trade-off between the superior angular resolution of GALFA-H i and the optimal sky coverage of H i4PI. This work focuses on H i4PI data, in order to study the 3D distribution of the magnetic ISM over the full sky. We also build maps from GALFA-H i, and compare the two in Section 5.

We bin the H i4PI data into velocity channels such that there is approximately equal integrated intensity in each pair of channels moving symmetrically outward from the center line. This binning scheme, illustrated in Figure 1, uses the native Δv = 1.3 $\mathrm{km}\,{{\rm{s}}}^{-1}$ channel width for velocity channels near v_lsr = 0, and channel widths of up to 20.6 $\mathrm{km}\,{{\rm{s}}}^{-1}$ in the outermost bins. This velocity channel binning is similar to the one used by Lenz et al. (2019). Throughout this work, we denote the integrated H i intensity in a velocity bin centered on velocity v with width dv as $I\left(v\right)$, i.e.,

$$\begin{eqnarray}&&I\left(v\right)={\int }_{v-{dv}/2}^{v+{dv}/2}{T}_{b}\left(v^{\prime} \right){dv}^{\prime} ,\end{eqnarray} \tag{ 1 }$$

where $I\left(v\right)$ is measured in K km s⁻¹. We also define the H i column density

$$\begin{eqnarray}&&{N}_{{\rm{H}}{\rm{I}}}=1.82\times {10}^{18}{\int }_{-90}^{+90}{T}_{b}\left(v\right){dv},\end{eqnarray} \tag{ 2 }$$

where N_{H i} is in cm⁻² for T_b in K and v in $\mathrm{km}\,{{\rm{s}}}^{-1}$, and the integral over −90 $\mathrm{km}\,{{\rm{s}}}^{-1}$ < v < +90 $\mathrm{km}\,{{\rm{s}}}^{-1}$ is chosen so that N_{H i} mostly represents Galactic emission. Equation (2) is correct under the assumption that the H i emission is optically thin—which is not valid for all sightlines in these data, but is generally a good approximation at high Galactic latitudes (Murray et al. 2018).

Zoom In Zoom Out Reset image size

We also use the publicly available H i intensity and RHT output of Δv = 3.7 $\mathrm{km}\,{{\rm{s}}}^{-1}$ GALFA-H i channel maps from −36.4 to +37.2 $\mathrm{km}\,{{\rm{s}}}^{-1}$ described in Peek et al. (2018), and convert to Galactic coordinates. Maps constructed from GALFA-H i data not only use independent 21 cm observations that are higher-resolution than H i4PI, but also a different velocity range and binning scheme. The GALFA-H i observations span −1° ≲ decl. ≲ +38°, but we make a conservative cut on the area considered in this work to 15 < decl. < 355 to avoid telescope scan artifacts at the edges of the Arecibo decl. range.

We make use of several data products released by the Planck collaboration for comparison with our H i-based maps. We use the 80' R3.00 353 GHz Stokes I, Q, and U maps post-processed with the Generalized Needlet Internal Linear Combination algorithm (GNILC; Remazeilles et al. 2011) to remove the anisotropic Cosmic Infrared Background (CIB) from the Galactic dust emission. Following the fiducial offset correction adopted by the Planck collaboration, we subtract 452 μK_CMB from the GNILC Stokes I map to correct for the CIB monopole at 353 GHz. We also add a Galactic offset correction of 63 μK_CMB to the I map (see discussion in Planck Collaboration XII 2018). We compute the modified asymptotic estimator introduced by Plaszczynski et al. (2014) to create maps of the noise-debiased polarization fraction. The noise debiasing is a minor adjustment at 80' resolution. Hereafter, we refer to this 80', debiased, offset-corrected, GNILC-derived polarization fraction map as p₃₅₃.

For the power spectra comparisons in Section 6.5, we employ the 353 GHz R3.01 half-mission I, Q, and U maps. We subtract the CMB contribution from these maps using the R3.00 SMICA CMB maps, and then smooth to 162 to match the native resolution of the H i4PI data.

In this work, we present a framework for mapping the three-dimensional structure of the magnetic neutral ISM using spectroscopic observations of the 21 cm line. We rely on three observational facts that relate the structure of H i to the properties of dust emission. These are as follows:

• 1.  
  The H i column density traces dust column density in the diffuse ISM (e.g., Lenz et al. 2017).
• 2.  
  H i features in narrow spectral channel maps are aligned with the plane-of-sky magnetic field (Clark et al. 2015).
• 3.  
  The coherence of H i orientation as a function of velocity traces the degree of line-of-sight magnetic field tangling (Clark 2018).

Dust and gas are well-mixed in the diffuse ISM as evidenced by the tight empirical correlation between the H i column density (${N}_{{\rm{H}}{\rm{I}}}$) and the dust extinction and emission (Burstein & Heiles 1978; Boulanger et al. 1996). Although this relationship breaks down at higher column densities, where molecular gas can constitute a significant fraction of the hydrogen column, ${N}_{{\rm{H}}{\rm{I}}}$ is a reliable tracer of the total dust column for low-column lines of sight. Recent analysis of the H i4PI data found that ${N}_{{\rm{H}}{\rm{I}}}$ and the dust reddening $E(B-V)$ are well-fit by a linear relationship for ${N}_{{\rm{H}}{\rm{I}}}$ < 4 × 10²⁰ cm⁻² with a scatter of ∼10% (Lenz et al. 2017).

High-resolution images of H i emission display prominent linear structure that is aligned well with the plane-of-sky magnetic field, as traced by both optical starlight polarization (Clark et al. 2014) and FIR polarized dust emission (Clark et al. 2015; Martin et al. 2015). Magnetically aligned H i intensity structures are a ubiquitous feature of the diffuse ISM. The aligned structures are thought to be components of the cold neutral medium (CNM) on the basis of temperature estimates from linewidth measurements (Clark et al. 2014; Kalberla et al. 2016), the existence of similar structures observed in H i absorption in the Galactic plane (McClure-Griffiths et al. 2006), an enhanced FIR/${N}_{{\rm{H}}{\rm{I}}}$ ratio at the locations of these structures in the diffuse ISM (Clark et al. 2019), and the dependence of the equivalent width of interstellar Na i absorption on the column density in small-scale channel map structures (Peek & Clark 2019).

Similar alignment between diffuse FIR dust intensity structures and the projected magnetic field was measured by Planck (Planck Collaboration Int. XXXII 2016). This alignment is thought to give rise to some of the statistical properties of the Planck polarized dust emission, notably the positive TE correlation and the nonunity EE/BB ratio (Planck Collaboration Int. XXXVIII 2016). Indeed, template E- and B-mode maps constructed from H i orientation alone (with no information on the polarized intensity) measure EE/BB ∼ 2 in the diffuse ISM, in rough agreement with the Planck measurement (Clark et al. 2015). If the E/B asymmetry and positive TE correlation are a consequence of the alignment of density structures and the magnetic field, we expect that our maps will naturally reproduce these statistical properties (see also Ghosh et al. 2017; Adak et al. 2019).

The final observational point allows us to model variable magnetic field orientations along the line of sight in a data-driven manner. Clark (2018) computed a metric, termed "H i coherence," quantifying the degree of order or disorder in H i orientation in different velocity channels along the same line of sight. Lines of sight where the orientations of H i features are relatively aligned across velocity channels correspond to an ordered field. Conversely, sightlines with relatively misaligned H i features in different velocity channels correspond to sightlines probing multiple field orientations. Indeed, the H i coherence was found to be a strong predictor of the Planck 353 GHz fractional polarization (p₃₅₃) in the diffuse ISM.

We emphasize that H i LOS velocity is not a one-to-one probe of distance, particularly at high Galactic latitude where the H i is not strongly sheared by Galactic rotation. If two differently oriented H i features lie at two velocities along a single line of sight, it is typically not possible to say from the H i data alone where the two features exist along the distance axis. However, since H i orientation traces the plane-of-sky magnetic field orientation, these two H i features imply that two regions with differently oriented magnetic fields lie somewhere along the line of sight. For the purposes of modeling the magnetic field in the neutral medium traced by polarized dust emission, it is sufficient to partition the H i intensity into regions with distinct magnetic field orientations.

Based on the insights outlined in Section 3.1, we posit that H i emission features that are coherent in velocity, plane-of-sky orientation, and spatial extent represent magnetically and spatially coherent clouds. This idea is shown schematically in Figure 2: two clouds along a single line of sight that sit in spatially distinct regions with differently oriented magnetic fields will map to different regions of velocity-orientation space. Observations of the 21 cm line emission (GALFA-H i or H i4PI) map the H i intensity as a function of sky position and line-of-sight velocity, ${I}_{{\rm{H}}{\rm{I}}}(l,b,v)$, where (l, b) represent sky coordinates and v represents the line-of-sight velocity. By measuring the orientation θ of linear structures in each velocity channel, we wish to map this emission into (l, b, v, θ) space. The distribution of H i emission in (l, b, v, θ) traces "magnetic coherence" without cloud-finding or otherwise prescribing dust properties in three-dimensional space.

Zoom In Zoom Out Reset image size

Mapping ${I}_{{\rm{H}}{\rm{I}}}(l,b,v)$ to ${I}_{{\rm{H}}{\rm{I}}}(l,b,v,\theta )$ requires an algorithm that can measure the distribution of H i as a function of orientation on the sky. We use the Rolling Hough Transform (RHT; Clark et al. 2014) to quantify the coherent linearity of H i emission in each channel map. The RHT maps image-plane data into position–position-orientation space in four steps: (1) high-pass filtering and bitmasking the image to highlight structure on small spatial scales; (2) selecting a circular window of data with diameter D_W centered on each image pixel; (3) calculating the fraction of the window that is nonzero in the bitmask as a function of orientation, R(θ); and (4) thresholding R(θ) at a prespecified fraction of the window length, Z. The algorithm used for step 3 is based on the Hough transform (Hough 1962). The basic principle is that straight lines in (x, y) image space are represented by points in (ρ, θ) space, given the mapping

$$\begin{eqnarray}&&\rho =x\,\cos \theta +y\,\sin \theta .\end{eqnarray} \tag{ 3 }$$

In the RHT, each (x, y) image position within a given circular window is mapped to the space of possible lines with ρ = 0: that is, lines that pass through the center of the window.

The high-pass filtering step of the RHT is equivalent to an unsharp mask. In our analysis of the H i4PI data, we use a Gaussian filter with FWHM = 30', following Kalberla et al. (2016). We use D_W = 75' and Z = 0.7 for consistency with Clark et al. (2015). We use the canonical θ-binning of Clark et al. (2014), but note that the alternative binning proposed by Schad (2017) does not substantially change the results. To avoid distortion effects from the map projection, we project each rolling window such that the center of the window is the center of the projection, as was done for the GALFA-H i DR2 data in Peek et al. (2018). The RHT works on any two-dimensional data: apart from H i observations, it has been used to quantify the orientation of structures for such disparate applications as dust emission (Malinen et al. 2016), the solar corona (Asensio Ramos et al. 2017), depolarization canals in radio polarimetric observations (Jelić et al. 2018), X-ray data (Marelli et al. 2019), and magnetohydrodynamic (MHD) simulations (Inoue & Inutsuka 2016).

By performing the RHT on each velocity channel of the H i data, we calculate R(v, θ), the linear intensity as a function of orientation and line-of-sight velocity. We normalize the RHT amplitude at each velocity v such that

$$\begin{eqnarray}&&\displaystyle \sum _{\theta }R\left(v,\theta \right)d\theta =1.\end{eqnarray} \tag{ 4 }$$

We can therefore treat $R\left(v,\theta \right)$ as analogous to a probability distribution function over the possible orientations θ in a given pixel at a given velocity. The orientation bins are uniformly spaced. Note that Equation (4) only applies where the RHT output is nonzero: some velocity channels will have ${\sum }_{\theta }R\left(v,\theta \right)d\theta =0$ if no prominently linear structure is detected.

The mapping to velocity-orientation space is the backbone of the magnetic coherence paradigm illustrated in Figure 2. We show an example of this mapping using the RHT applied to H i4PI data in Figure 3.

Zoom In Zoom Out Reset image size

With measurements of both $I\left(v\right)$ and $R\left(v,\theta \right)$, we can construct synthetic Stokes parameters from the H i data. We assume only that the polarized intensity is proportional to the H i intensity and that the emission is polarized perpendicular to the orientation of H i filaments, as is expected from polarized dust emission. The Stokes parameters of emission with these characteristics are given by

$$\begin{eqnarray}&&{Q}_{{\rm{H}}{\rm{I}}}\left(v\right)=I\left(v\right)\displaystyle \sum _{\theta }R\left(v,\theta \right)\cos \left(2\theta \right)d\theta \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{U}_{{\rm{H}}{\rm{I}}}\left(v\right)=I\left(v\right)\displaystyle \sum _{\theta }R\left(v,\theta \right)\sin \left(2\theta \right)d\theta .\end{eqnarray} \tag{ 6 }$$

Because they are derived from the $I\left(v\right)$ data, ${Q}_{{\rm{H}}{\rm{I}}}\left(v\right)$ and ${U}_{{\rm{H}}{\rm{I}}}\left(v\right)$ have the same units as ${I}_{{\rm{H}}{\rm{I}}}\left(v\right)$, K km s⁻¹.

The orientation ${\theta }_{{\rm{H}}{\rm{I}}}\left(v\right)$ and "polarization fraction" ${p}_{{\rm{H}}{\rm{I}}}\left(v\right)$ of the emission in each velocity channel are computed from the Stokes parameters in the usual way:

$$\begin{eqnarray}&&{\theta }_{{\rm{H}}{\rm{I}}}\left(v\right)=\displaystyle \frac{1}{2}\arctan \displaystyle \frac{{U}_{{\rm{H}}{\rm{I}}}\left(v\right)}{{Q}_{{\rm{H}}{\rm{I}}}\left(v\right)}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{p}_{{\rm{H}}{\rm{I}}}\left(v\right)=\displaystyle \frac{\sqrt{{Q}_{{\rm{H}}{\rm{I}}}{\left(v\right)}^{2}+{U}_{{\rm{H}}{\rm{I}}}{\left(v\right)}^{2}}}{I\left(v\right)}.\end{eqnarray} \tag{ 8 }$$

The additivity of the Stokes parameters permits straightforward integration over velocity channels:

$$\begin{eqnarray}&&{I}_{{\rm{H}}{\rm{I}}}=\displaystyle \sum _{v}I\left(v\right)\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{Q}_{{\rm{H}}{\rm{I}}}=\displaystyle \sum _{v}{Q}_{{\rm{H}}{\rm{I}}}\left(v\right)\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{U}_{{\rm{H}}{\rm{I}}}=\displaystyle \sum _{v}{U}_{{\rm{H}}{\rm{I}}}\left(v\right).\end{eqnarray} \tag{ 11 }$$

In the magnetically coherent cloud paradigm central to this work, the sum over H i velocities and orientations is analogous to the sum over distinct regions along the line of sight.

The orientations ${\theta }_{{\rm{H}}{\rm{I}}}\left(v\right)$ and polarization fractions ${p}_{{\rm{H}}{\rm{I}}}\left(v\right)$ are not additive, and so the velocity integrated orientation and polarization fraction are given by

$$\begin{eqnarray}&&{\theta }_{{\rm{H}}{\rm{I}}}=\displaystyle \frac{1}{2}\arctan \displaystyle \frac{{U}_{{\rm{H}}{\rm{I}}}}{{Q}_{{\rm{H}}{\rm{I}}}}\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{p}_{{\rm{H}}{\rm{I}}}=\displaystyle \frac{\sqrt{{Q}_{{\rm{H}}\,{\rm{I}}}^{2}+{U}_{{\rm{H}}\,{\rm{I}}}^{2}}}{{I}_{{\rm{H}}{\rm{I}}}}.\end{eqnarray} \tag{ 13 }$$

Also by analogy with dust emission, the polarized intensity is defined as

$$\begin{eqnarray}&&{P}_{{\rm{H}}{\rm{I}}}={p}_{{\rm{H}}{\rm{I}}}{I}_{{\rm{H}}{\rm{I}}}.\end{eqnarray} \tag{ 14 }$$

While we demonstrate that these H i-derived quantities are effective predictors of the corresponding quantities in dust emission, a few key differences exist. First, the H i intensity is insensitive to changes in the dust-to-gas ratio and dust temperature, so any variations in those quantities in the Galactic ISM will reduce the correlation with the observed dust emission.

Second, because dust grains align with their rotation axis parallel to the local magnetic field, the polarized dust emission depends on the orientation of the GMF relative to the line of sight. If γ is the angle between the magnetic field and the plane of the sky, then the polarized intensity of the dust emission is proportional to ${\cos }^{2}\gamma$ (Lee & Draine 1985). Although this angle is not directly measurable from the H i data, it is closely related to the dispersion of polarization angles such that more dispersion is expected when the GMF is along the line of sight (γ = π/2) and less when it lies in the plane of the sky (γ = 0). This effect is captured in our maps by the sum over the distribution of orientations (Equations (5) and (6)). Furthermore, if linear H i structures are elongated in the direction of the three-dimensional magnetic field, they will be less well-detected when γ ∼ π/2, as this will diminish their plane-of-sky extent. We thus expect that information about γ is encoded in our maps implicitly.

We first demonstrate in Figure 9 that the velocity-orientation space mapping is predictive of the fractional polarization of each sightline, by showing that the distribution of R(v, θ) is sensitive to the measured Planck 353 GHz polarization fraction, p₃₅₃. We bin the maps into quintiles of p₃₅₃, and observe that the R(v, θ) distribution becomes more narrowly peaked about its mean for sightlines with higher polarization fractions. This confirms our intuition that lines of sight with H i that is less coherent in velocity-orientation space will be more depolarized, as emission from dusty regions with different magnetic field orientations adds vectorially.

Zoom In Zoom Out Reset image size

Figure 8 includes all-sky maps of both p₃₅₃ and p_{H i}, and we show a two-dimensional histogram of the two polarization fractions in Figure 10. The two quantities are strongly linearly correlated, with appreciable scatter. A simple linear regression of p₃₅₃ versus p_{H i} over the whole sky yields a correlation coefficient r ∼ 0.6. As p_{H i} is derived solely from the morphology of I(v), this correlation indicates that much of the variation in p₃₅₃ arises purely from the geometry of the magnetic field rather than, e.g., variable grain alignment, in good agreement with the findings in Planck Collaboration XII (2018) and Clark (2018).

Zoom In Zoom Out Reset image size

The maximum value of the observed polarization fraction is an important constraint on the intrinsic polarizing efficiency of dust grains (Draine & Fraisse 2009; Guillet et al. 2018). Planck Collaboration XII (2018) reported a maximum polarization fraction ${p}_{353}^{\max }\sim 0.22$. This implies that dust grain populations must emit thermal radiation that is intrinsically at least ∼ 22% polarized, before effects like magnetic field tangling induce depolarization. By construction, p_{H i} is sensitive only to these geometrical effects and not to the intrinsic polarization efficiency of dust. One might therefore expect that the most coherent sightlines have p_{H i} ∼ 1, but Figures 8 and 10 show that the range of p_{H i} is similar to p₃₅₃: the 99.9th percentile value of p_{H i} is ∼0.26 at 80' resolution.

If p_{H i} were a numerically accurate predictor of the effects of geometric depolarization, this would imply that dust emission must have an intrinsic polarizing efficiency of nearly unity. However, the alignment of H i gas with magnetic field lines is imperfect, and so we generically expect the H i orientations to be more disordered than the true magnetic field. As a consequence, p_{H i} is an overestimate of the geometric depolarization present in the dust maps. Another reason p_{H i} may overestimate the amount of depolarization is that the denominator of Equation (13) includes all of the H i emission along the line of sight, while Q_{H i} and U_{H i} are weighted only by the emission in linear H i features. Since there are some velocity bins in which R(v, θ) is zero for all θ, these bins contribute to the total but not polarized intensity (Section 3.2). This likely undercounts the fraction of H i emission that is correlated with polarized dust emission, as discussed further in Section 9. The numerical proximity of ${p}_{{\rm{H}}\,{\rm{I}}}^{\max }$ to ${p}_{353}^{\max }$ appears therefore to be a coincidence of the overestimation of depolarization compensating for the absence of an intrinsic polarization fraction in the model.

We can compare the estimate for the magnetic field orientation between two maps by computing the angular difference (Equation (15)) in each pixel, where θ₁ in this case is θ_{H i} computed from the relevant velocity range, and θ₂ is θ₃₅₃. We collapse this δθ distribution to a point statistic for the mean degree of alignment,

$$\begin{eqnarray}&&\xi =\left\langle \cos \phi \right\rangle ,\end{eqnarray} \tag{ 16 }$$

where

$$\begin{eqnarray}&&\phi =2\delta \theta .\end{eqnarray} \tag{ 17 }$$

When θ_{H i} = θ₃₅₃, ϕ = 0, and when θ_{H i} and θ₃₅₃ are orthogonal to one another, ϕ = π. The metric ξ is therefore defined on [−1, 1], such that two perfectly aligned distributions will have ξ = 1, two antialigned distributions will have ξ = −1, and two distributions with no statistical alignment will have ξ = 0. This metric differs from the "projected Rayleigh statistic," (PRS; Jow et al. 2018), by a factor of $\sqrt{2N}$, where N is the number of samples (pixels). That is,

$$\begin{eqnarray}&&\mathrm{PRS}=\sqrt{\displaystyle \frac{2}{N}}\displaystyle \sum _{i}\cos \,{\phi }_{i},\end{eqnarray} \tag{ 18 }$$

and the uncertainty can be estimated as

$$\begin{eqnarray}&&{\sigma }_{\mathrm{PRS}}^{2}=\displaystyle \frac{2}{N}\left[\displaystyle \sum _{i}{\cos }^{2}\,{\phi }_{i}-{(\mathrm{PRS})}^{2}\right].\end{eqnarray} \tag{ 19 }$$

The PRS is equivalent to the modified Rayleigh test for uniformity proposed by Durand & Greenwood (1958), for the specific case where the mean angle of the distribution is ϕ = 0. We define ξ for its convenient range, but will also report the PRS where appropriate. For the full-sky θ_{H i} and θ₃₅₃ maps shown in Figure 8, ξ = 0.71. For the 80' maps pixelated on an N_side = 128 HEALPix grid, PRS = 446 and σ_PRS = 0.5.

One measure of the degree of order in the plane-of-sky magnetic field orientation is the polarization angle dispersion function, ${ \mathcal S }$ (Planck Collaboration Int. XIX 2015). Following the notation used by the Planck collaboration, we define

$$\begin{eqnarray}&&S({\boldsymbol{r}},\delta )=\sqrt{\displaystyle \frac{1}{N}\displaystyle \sum _{i=1}^{N}{\left[\psi ({\boldsymbol{r}}+{\delta }_{i})-\psi ({\boldsymbol{r}})\right]}^{2}},\end{eqnarray} \tag{ 20 }$$

where the sum is over pixels located within an annulus of inner radius = δ/2 and outer radius = 3δ/2. In practice, ${ \mathcal S }$ is computed in terms of the Stokes parameters Q and U, and in order to avoid a spuriously high ${ \mathcal S }$ near the poles from the polarization reference frame, we calculate ${ \mathcal S }$ such that each annulus sits at the equator of its reference frame. Planck Collaboration XII (2018) found that polarization systematics bias the measurement of ${ \mathcal S }$ at high Galactic latitudes unless ${ \mathcal S }$ is computed at a resolution of 160' and lag δ = 80', so we adopt this resolution and lag for both ${{ \mathcal S }}_{{\rm{H}}{\rm{I}}}$ and ${{ \mathcal S }}_{353}$.

Figure 8 includes all-sky maps of ${{ \mathcal S }}_{{\rm{H}}{\rm{I}}}$ and ${{ \mathcal S }}_{353}$. The ${{ \mathcal S }}_{{\rm{H}}{\rm{I}}}$ parameter derived from Q_{H i} and U_{H i} is broadly similar to ${{ \mathcal S }}_{353}$ derived from the Planck observations, particularly on large angular scales. In Figure 11, we examine the relationship between p_{H i} and ${{ \mathcal S }}_{{\rm{H}}{\rm{I}}}$. We find a trend similar to the relationship that Planck Collaboration XII (2018) reported for p₃₅₃ and ${{ \mathcal S }}_{353}$: there exists a negative correlation between ${ \mathcal S }$ and p (see also Fissel et al. 2016). When the polarization angles used to compute a value of ${ \mathcal S }$ are completely randomly oriented, ${ \mathcal S }=\pi /\sqrt{12}$ (Planck Collaboration Int. XIX 2015; Planck Collaboration XII 2018). This is why the values of ${{ \mathcal S }}_{{\rm{H}}{\rm{I}}}$ asymptote near 52°.

Zoom In Zoom Out Reset image size

We also examine the correlation between the magnetic field orientation of our maps, θ_{H i}, and the Planck θ₃₅₃ measurement, as a function of the polarization angle dispersion function (Figure 12). We compute Equation (16) for data binned by ${{ \mathcal S }}_{{\rm{H}}{\rm{I}}}$ value. We find a striking inverse correlation between ${{ \mathcal S }}_{{\rm{H}}{\rm{I}}}$ and ξ, indicating that in regions of sky where the plane-of-sky magnetic field orientation is relatively ordered within a 160' beam (low ${{ \mathcal S }}_{{\rm{H}}{\rm{I}}}$), the magnetic field orientation in our H i-based maps is well-aligned with the magnetic field orientation inferred from the 353 GHz dust measurements (high ξ). Conversely, regions of high polarization angle dispersion see a significant decrease in the statistical alignment of the two sets of angles. In all bins, θ₃₅₃ and θ_{H i} are significantly aligned: PRS = 86, σ_PRS = 0.01 in the lowest ${{ \mathcal S }}_{{\rm{H}}{\rm{I}}}$ bin, and PRS = 14, σ_PRS = 1 in the highest ${{ \mathcal S }}_{{\rm{H}}{\rm{I}}}$ bin, computed at N_side = 128. As discussed further in Section 8.4, this is qualitatively consistent with the hypothesis that at 160', regions of low ${{ \mathcal S }}_{{\rm{H}}{\rm{I}}}$ are preferentially regions where the mean magnetic field orientation lies in the plane of the sky.

Zoom In Zoom Out Reset image size

The polarized emission at 353 GHz is optically thin, so the Planck observations trace dust emission integrated over the entire line of sight. With our velocity-resolved Stokes parameters, we can examine the correlation between the Planck measurements and the H i-based maps as a function of the H i velocity range considered, by changing the bounds of the sums in Equations (10) and (11). The data are summed over an ever-widening velocity range that is always centered at the same velocity bin (v ∼ 0 $\mathrm{km}\,{{\rm{s}}}^{-1}$).

We examine the Spearman correlation coefficient between Q_{H i} and Q₃₅₃ as a function of the velocity range used in this sum. We find a monotonic increase in the correlation, with a steep dependence on the velocity range until about 30 $\mathrm{km}\,{{\rm{s}}}^{-1}$, followed by a gradual improvement of the correlation out to the full velocity range. We repeat this experiment, but compare θ_{H i} and θ₃₅₃ by computing Equation (16). We find the same qualitative dependence on the velocity range for Equation (16) as we do for the Q_{H i}–Q₃₅₃ correlation. Figure 13 shows these correlations over the full sky, but we find a similar trend when restricting the latitude range to $| b| \gt 30^\circ$ or $| b| \gt 60^\circ$.

Zoom In Zoom Out Reset image size

There is, on average, less Galactic H i emission far from the 21 cm rest wavelength, and the distribution of H i intensity along the line of sight is incorporated into Q_{H i}(v) and U_{H i}(v) via the I(v) weighting in Equations (5) and (6). There do exist H i-bright structures at large absolute velocities, notably high-velocity clouds (HVCs; Putman et al. 2012) and nearby galaxies, including the Magellanic Clouds. In general, we expect contamination from these extragalactic objects to decrease the correlation between our H i-based maps and the Planck data, as the component-separated 353 GHz maps in principle trace only Galactic dust. This is particularly true in the case of HVCs, which are H i rich but dust-deficient (Wakker & Boulanger 1986; Planck Collaboration XXIV 2011; Lenz et al. 2017). We have purposefully restricted the velocity range of the H i gas that we consider in this work to −93.9 < v < 96.7 $\mathrm{km}\,{{\rm{s}}}^{-1}$ in order to avoid the bulk of this high-velocity emission. This is consistent with the velocity range identified by Lenz et al. (2017) for which H i column density and dust reddening are most strongly correlated. The monotonic increase of the quantities in Figure 13 suggests that our velocity cutoff is sufficient to avoid significant decorrelation from extragalactic gas, but this statement is only valid as a global average, and individual sightlines may contain significant emission associated with non-Galactic gas. It is straightforward to restrict the velocity range of our 3D maps for studies of such regions.

In the previous sections, we demonstrated that the H i-derived Stokes maps reproduce many key features of the 353 GHz polarized dust emission at the map level. The interest in polarized dust emission as a foreground for CMB science has necessitated characterization of polarized dust emission at the power spectrum level as well, and so in this section we examine the properties of the TE, EE, and BB spectra.

We compute the power spectra ${{ \mathcal D }}_{{\ell }}^{{XY}}\equiv {\ell }({\ell }+1){C}_{{\ell }}/2\pi$ for our H i-derived Stokes maps, where XY are TE, EE, and BB. Here, ${C}_{{\ell }}$ represents the pure-E and pure-B pseudo-C_ℓ estimator, computed to eliminate E–B mixing on our cut-sky maps (Smith 2006). All cross spectra are computed with the NaMaster software (Alonso et al. 2019). We restrict our analysis to $| b| \gt 30^\circ$, but find that our results are stable to more aggressive masking. To facilitate comparison with the Planck 353 GHz power spectra, we mask all pixels in the Planck half-mission I, Q, and U maps that have been set to bad_value. In addition, we mask all pixels for which ${\sum }_{v}{\sum }_{\theta }R(v,\theta )=0$. This mask is apodized with the analytic C² method of Grain et al. (2009), with an apodization scale of 5°. These cuts result in a sky fraction f_sky = 41%.

In Figure 14, we present the three cross spectra binned with Δℓ = 10. The error bars are computed assuming simple Gaussian sample variance and no instrumental noise, i.e.,

$$\begin{eqnarray}&&{\sigma }_{{{ \mathcal D }}_{{\ell }}^{{XY}}}=\sqrt{\displaystyle \frac{{{ \mathcal D }}_{{\ell }}^{{XX}}{{ \mathcal D }}_{{\ell }}^{{YY}}+{\left({{ \mathcal D }}_{{\ell }}^{{XY}}\right)}^{2}}{{f}_{\mathrm{sky}}\left(2{\ell }+1\right){\rm{\Delta }}{\ell }}},\end{eqnarray} \tag{ 21 }$$

and so are somewhat underestimated.

Zoom In Zoom Out Reset image size

On the largest scales (ℓ ≲ 50), the cross spectra have a shallow negative slope. At higher ℓ, and unlike what is observed in the dust emission (Planck Collaboration XI 2018), the cross spectra begin to rise with increasing ℓ, suggesting an excess of small-scale power in the H i maps relative to the dust emission. This behavior is perhaps unsurprising because all polarized emission in this model is explicitly attributed to small-scale, filamentary structures. Indeed, this effect is most pronounced in EE, as expected from this interpretation. We defer a more detailed investigation of the small-scale structure to future work.

One of the most surprising findings about Galactic polarized dust emission is that it has roughly twice as much power in E-mode polarization as B-mode (Planck Collaboration Int. XXX 2016). The prevailing physical explanation is that the preferential elongation of density structures along the local magnetic field direction is responsible for this excess E-mode polarization, as well as for the observed positive TE correlation (Clark et al. 2015; Planck Collaboration Int. XXXVIII 2016; Huffenberger et al. 2019). Indeed, Clark et al. (2015) found that the E/B asymmetry persisted for polarization templates built solely from maps of H i orientation, without including polarized intensity information.

The E/B asymmetry is readily apparent in Figure 14, with BB/EE ≃ 0.6 on large angular scales (ℓ ≲ 120). We also compute the BB/EE ratio derived from power-law fits to the Planck 353 GHz half-mission split cross-power spectra over the same mask and ℓ range, and likewise find BB/EE ∼ 0.6 (see bottom panel of Figure 14). No free parameters were adjusted to achieve this correspondence, strongly suggesting that the observed asymmetry indeed arises from filamentary density structures that are oriented along magnetic field lines. The ratio found here is comparable to—but slightly higher than—the Planck value observed over large sky areas, BB/EE = 0.53 ± 0.01 (Planck Collaboration XI 2018). The sky-variable BB/EE ratio may be a useful probe of the physics of the dusty ISM, perhaps related to local properties of turbulence (e.g., Caldwell et al. 2017).

We quantify the scale-dependent correlation of the H i-based maps with the Planck 353 GHz maps in Figure 15. As described in Section 2.2, we analyze the 353 GHz Planck half-mission maps that have been CMB-subtracted and smoothed to a resolution of 162. We plot the correlation ratios r defined as

$$\begin{eqnarray}&&{r}_{{\rm{H}}\,{\rm{I}}\times 353}^{{XX}}\equiv \displaystyle \frac{{{ \mathcal D }}_{{\ell }}^{{X}_{{\rm{H}}{\rm{I}}}{X}_{353}}}{\sqrt{{{ \mathcal D }}_{{\ell }}^{{X}_{{\rm{H}}{\rm{I}}}{X}_{{\rm{H}}{\rm{I}}}}\times {{ \mathcal D }}_{{\ell }}^{{X}_{353}{X}_{353}}}},\end{eqnarray} \tag{ 22 }$$

where the subscripts H i and 353 denote the H i and Planck 353 GHz maps, respectively, and X denotes either E or B. The correlation between the maps is greatest at large scales, with r > 0.3 for all ℓ < 100 and reaching values as high as 0.75 and 0.60 for EE and BB, respectively. For both cross spectra, r declines roughly monotonically toward small scales where, as was noted in Figure 14, the H i maps appear to have more power than observed in the dust maps.

Zoom In Zoom Out Reset image size

Observations of diffuse synchrotron emission trace the Faraday-rotated magneto-ionic medium. A new generation of radio polarimetric observations is opening a window into magnetic fields in the ionized ISM, with facilities like the Low Frequency Array (LOFAR; van Haarlem et al. 2013), the Murchison Wide-field Array (Tingay et al. 2013), and the multi-instrument project to map the full sky, the Global Magneto-Ionic Medium Survey (GMIMS; Wolleben et al. 2019).

LOFAR observations of a field centered on the quasar 3C 196 revealed strikingly linear, coherent structures in Faraday depth, including one prominent filament several degrees in length (Jelić et al. 2015). These magneto-ionic structures are aligned with H i filaments in the same field (Kalberla & Kerp 2016; Jelić et al. 2018), as well as the Planck 353 GHz magnetic field orientation (Zaroubi et al. 2015). The LOFAR observations are sensitive to magnetic structure in the warm ionized medium, while the dust and H i are components of the neutral phase. Why does a morphological similarity persist across such different phases of the ISM? Jelić et al. (2018) propose that the various ISM phases are confined by a well-ordered magnetic field that lies predominantly along the plane of the sky. Those authors also observe alignment of H i filaments across many velocity channels, indicating coherence of the magnetic field orientation along the line of sight (Clark 2018).

Here, we examine the Jelić et al. (2018) picture in the context of our three-dimensional Stokes maps. Figure 16 shows the sum of I(v)R(v, θ) over a circular region of diameter 5°, centered on 3C 196 at (l, b) = (171°, b = 33°). Our maps show a remarkably coherent orientation of linear H i structures across a range of v_los, with a dominant H i orientation ($\left\langle {\theta }_{{\rm{H}}{\rm{I}}}\right\rangle \sim 4^\circ$ over the full velocity range) that is consistent with the orientation of the Faraday depolarization canals measured by Jelić et al. (2018) (${\left\langle \theta \right\rangle }_{\mathrm{LOFAR}}=10^\circ \pm 6^\circ$).⁵ If the magnetic field is coherent and mostly in the plane of the sky, we expect our maps to have a high p_{H i} and a low ${{ \mathcal S }}_{{\rm{H}}{\rm{I}}}$ in this region. Indeed, averaging our H i-based Stokes parameters over this circular region, we find $\left\langle {p}_{{\rm{H}}{\rm{I}}}\right\rangle =11.5 \%$, higher than the full-sky average $\left\langle {p}_{{\rm{H}}{\rm{I}}}\right\rangle =9.2 \%$. The value of ${{ \mathcal S }}_{{\rm{H}}{\rm{I}}}$ is likewise low (53) over this region. Jelić et al. (2018) found that EBHIS H i filaments were oriented in the same direction in this field over −11.5 < v < +3 $\mathrm{km}\,{{\rm{s}}}^{-1}$. Not surprisingly, Figure 1 shows that the H i4PI data (derived from EBHIS data in this region) are coherent in orientation space over this velocity range, with a mean orientation $\left\langle {\theta }_{{\rm{H}}{\rm{I}}}\right\rangle \sim 12^\circ$. If we restrict our sums in Equations (9)–(11) to approximately this same velocity range, we find for this region $\left\langle {p}_{{\rm{H}}{\rm{I}}}\right\rangle =46.7 \%$ (see discussion of p_{H i} in Section 6.1). The mean Planck polarization fraction in this region is $\left\langle {p}_{353}\right\rangle =9.6 \%$.

Zoom In Zoom Out Reset image size

Faraday depth, like H i velocity, is not a direct probe of distance. However, both this work and Jelić et al. (2018) suggest that velocity- and Faraday depth-resolved probes of the magnetic field structure can be used in tandem to constrain the relative distances to magnetized structures. In conjunction with observations that more directly map to distance, these data can be used to map the magnetic field in three spatial dimensions, as well as to better understand the relationship between different phases of the magnetic ISM. Tomographic analyses combining probes of the neutral ISM with LOFAR observations (Van Eck et al. 2017) and GMIMS observations (Thomson et al. 2019) have already been successful in selected regions of sky. We discuss further possibilities for three-dimensional magnetic tomography below.

Starlight polarization probes the plane-of-sky magnetic field orientation in the dusty medium between observer and star (Davis & Greenstein 1951). As an integral probe, starlight polarization can be used to estimate the incremental polarization as a function of distance, given a sufficient density of polarimetric measurements and accurate stellar distances (e.g., Lloyd & Harwit 1973). Recently, Panopoulou et al. (2019) applied this principle to optical starlight polarimetry obtained with RoboPol (Ramaprakash et al. 2019), for stars with distances estimated from Gaia parallax measurements (Bailer-Jones et al. 2018; Gaia Collaboration et al. 2018).

Panopoulou et al. (2019) selected a circular region of sky of radius 016, centered on (l, b) = (10408, 2231). Dubbed the "two-cloud" region, this sightline passes through two H i emission structures that are distinct in velocity space: an intermediate-velocity cloud (IVC) and a low-velocity cloud (LVC; Figure 17). The polarimetric measurements toward stars in this region are compared to those in a control region of the same angular size centered on (l, b) = (10390, 2197) that contains less H i emission from intermediate velocities. The mean starlight polarization angles for stars that lie in each of these two regions agree within the error bars. However, Panopoulou et al. (2019) use stellar distances and the on-off setup of their selected regions to carefully disentangle the polarimetric properties associated with the near and far clouds. The authors find that the IVC is farther away, at a distance of 1250–2140 pc, and has a mean magnetic field orientation⁶ $\left\langle {\theta }_{\star }\right\rangle =106^\circ \pm 8^\circ$. For the LVC, they find a distance of 346–393 pc with $\left\langle {\theta }_{\star }\right\rangle =42\buildrel{\circ}\over{.} 6\pm 1^\circ$.

Zoom In Zoom Out Reset image size

To compare our three-dimensional Stokes parameter maps with the Panopoulou et al. (2019) polarimetric analysis, we select the velocity range in our maps closest to the ranges associated with the IVC and LVC: these are −51.4 < v < −38.5 $\mathrm{km}\,{{\rm{s}}}^{-1}$ and −3.8 < v < −1.2 $\mathrm{km}\,{{\rm{s}}}^{-1}$, respectively. We compute Q_{H i} and U_{H i} from Equations (10) and (11) over these velocity ranges, and visualize the resulting magnetic field orientation θ_{H i} in the two regions in Figure 17. We compute the mean magnetic field orientation in the two-cloud region defined above, finding $\left\langle {\theta }_{{\rm{H}}{\rm{I}}}\right\rangle =111\buildrel{\circ}\over{.} 6$ for the IVC and $\left\langle {\theta }_{{\rm{H}}{\rm{I}}}\right\rangle =42\buildrel{\circ}\over{.} 6$ for the LVC in our maps at 30'. These values agree well with the polarimetry-derived measurements for the two clouds. We find a higher $\left\langle {p}_{{\rm{H}}{\rm{I}}}\right\rangle$ for the Stokes maps integrated over the LVC velocity range than for the IVC velocity range, just as Panopoulou et al. (2019) derive a higher fractional linear polarization associated with the LVC than with the IVC.

It is worth emphasizing the remarkable correspondence between our three-dimensional Stokes maps and the Panopoulou et al. (2019) measurements. Our maps probe the local magnetic field as a function of velocity. Starlight polarization measures properties of the integrated magnetic field, which can be used in conjunction with stellar distance measurements to probe the magnetic field as a function of distance. The implications for higher-dimensional tomography of the magnetic ISM are discussed further in Section 9.

Mapping magnetic coherence requires some method of quantifying the orientation of H i structure. We use the RHT in this work, but our framework is not algorithm-specific. The RHT is not the only conceivable mapping from image space to orientation space, and there are a range of reasonable parameter choices for the RHT that could be adopted. Such differences would change the numerical values of Q(v) and U(v), although these are circumscribed by I(v), which is directly measured from the sky, and by the fact that the derived orientation in a given H i region (the ratio of Q and U) tends to be fairly insensitive to RHT parameter choices (Clark et al. 2014). We consider the maps provided here to be an excellent estimate of the three-dimensional Stokes parameters derived from H i emission, but not a unique solution. Another way to measure the orientation of image structures is simply to take a spatial gradient and then define the direction of maximum gradient as the local feature orientation. The spatial gradient is widely used in machine vision for edge detection, feature matching, and more complex applications, and has been used in astrophysical contexts to compare ISM structure to the magnetic field orientation (Soler et al. 2013).

Taking a spatial gradient is a computationally simpler procedure than the RHT, but it collapses any information about the distribution of H i emission in orientation space. If this distribution is simply summed along the θ axis as in Equations (5) and (6), then the final orientation angle will be similar as long as the local gradient orientation is similar to the mean RHT angle. A major caveat to this is that the RHT measures linearity as a function of orientation, and pixels will have zero R(v, θ) if there is insufficient linear intensity. This is a desirable property for building three-dimensional Stokes maps if the prominently linear H i structures are the most predictive of the polarized intensity. In contrast, the spatial gradient measures the local direction of maximum change in the H i intensity without regard for the feature that produced it. A circular blob of H i intensity could contribute nothing to R(v, θ) but have strong spatial gradients oriented in all directions around its edge. Without further processing, the spatial gradient is thus very sensitive to small-scale noise and other data features that are not predictive of magnetic field properties. This property is evident in the spatial gradient of H i4PI channel maps. Because the EBHIS and GASS surveys have different angular resolutions, one hemisphere was convolved with a larger beam to create the H i4PI data set at uniform resolution. The noise properties thus differ across the map, and this manifests in starkly different properties of the spatial gradient between the two hemispheres. Furthermore, the mapping to velocity-orientation space enabled by the RHT is a natural description of magnetic coherence that allows further exploration of dust properties (see Section 9).

We create an alternative set of three-dimensional Stokes maps using the spatial gradient of the H i emission I_{H i}(l, b, v). If the H i emission in a given velocity channel is defined on the sphere as I_{H i}(θ, ϕ), we take

$$\begin{eqnarray}&&{\rm{\nabla }}{I}_{{\rm{H}}{\rm{I}}}=\left[\begin{array}{c}{I}_{\theta }\\ {I}_{\phi }\end{array}\right]=\left[\begin{array}{c}{\partial }_{\theta }{I}_{{\rm{H}}{\rm{I}}}\\ \displaystyle \frac{1}{\sin \theta }{\partial }_{\phi }{I}_{{\rm{H}}{\rm{I}}}\end{array}\right]\end{eqnarray} \tag{ 23 }$$

and define the local feature orientation as orthogonal to the direction of the steepest intensity gradient,

$$\begin{eqnarray}&&{\theta }_{\mathrm{grad}}=\arctan \left(\displaystyle \frac{{I}_{\theta }}{{I}_{\phi }}\right).\end{eqnarray} \tag{ 24 }$$

In this approach, as with the RHT, we use an orientation convention that corresponds to the polarization convention used, such that the sense of the H i orientation is a direct proxy for the projected magnetic field orientation, which is orthogonal to the polarization angle of dust emission. We then compute

$$\begin{eqnarray}&&{Q}_{{\rm{H}}\,{\rm{I}}}^{\mathrm{grad}}\left(v\right)=I\left(v\right)\cos \left(2{\theta }_{\mathrm{grad}}\right)\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{U}_{{\rm{H}}\,{\rm{I}}}^{\mathrm{grad}}\left(v\right)=I\left(v\right)\sin \left(2{\theta }_{\mathrm{grad}}\right),\end{eqnarray} \tag{ 26 }$$

and Equations (9)–(11) hold, as in the RHT-based construction of the H i Stokes parameters. We then have Q_{H i}^grad, ${U}_{{\rm{H}}\,{\rm{I}}}^{\mathrm{grad}}$, and the derived quantities ${p}_{{\rm{H}}\,{\rm{I}}}^{\mathrm{grad}}$ and ${\theta }_{{\rm{H}}\,{\rm{I}}}^{\mathrm{grad}}$, which can be compared to both the Planck 353 GHz measurements and the RHT-based Q_{H i} and U_{H i} quantities considered in the rest of this work.

We find that the gradient-based Stokes parameters are similar to the RHT-derived values, but the RHT-derived maps are better correlated with the Planck measurements. For instance, a simple linear regression of Q_{H i} versus ${Q}_{{\rm{H}}\,{\rm{I}}}^{\mathrm{grad}}$ at 80' yields a correlation coefficient r = 0.91. This value is r = 0.87 for Q_{H i} versus Q₃₅₃, and 0.82 for ${Q}_{{\rm{H}}\,{\rm{I}}}^{\mathrm{grad}}$ versus Q₃₅₃. We find similar results when comparing derived quantities across the three sets of maps: ξ = 0.86 between θ_{H i} and ${\theta }_{{\rm{H}}\,{\rm{I}}}^{\mathrm{grad}}$, ξ = 0.71 between θ_{H i} and θ₃₅₃, and ξ = 0.67 between ${\theta }_{{\rm{H}}\,{\rm{I}}}^{\mathrm{grad}}$ and θ₃₅₃. Likewise, p₃₅₃ is better correlated with p_{H i} than with ${p}_{{\rm{H}}\,{\rm{I}}}^{\mathrm{grad}}$. The gradient-based maps have a lower correlation ratio (Equation (22)) with the 353 GHz maps than the RHT-based maps do (${r}_{{\ell }}^{{EE}}$ with the gradient is $\lesssim {r}_{{\ell }}^{{EE}}$ with the RHT for all ℓ).

It is further instructive to compare the small-scale structure of the gradient-based maps to the maps constructed using the RHT. The right-hand panel of Figure 6 shows P_{H i} constructed using Equations (25) and (26). Because the gradient is computed over small scales, the polarized intensity is not concentrated into structures as coherently linear as those in the RHT-based maps. Prominent H i filaments become bright P_{H i} filaments in the RHT-based maps, but it is often the edges of those filaments that contain strong gradients and thus show up as coherent structures in the ${P}_{{\rm{H}}\,{\rm{I}}}^{\mathrm{grad}}$ map.

In the diffuse ISM, the orientation of H i structures traces the plane-of-sky magnetic field with high fidelity. A fundamental assumption of the maps presented here is that the orientation of H i channel map structures is predictive of the magnetic field orientation over the whole sky. In other tracers, notably dust (Planck Collaboration Int. XXXII 2016) and molecular line emission (Goldsmith et al. 2008), there is evidence for a change in the relative orientation between structures and the magnetic field toward high-column sightlines. We find no evidence for a global decrease in the alignment between θ₃₅₃ and θ_{H i} at 80' as a function of ${N}_{{\rm{H}}{\rm{I}}}$ in the H i4PI data. Similar results were obtained for the correlation between θ₃₅₃ and ${\theta }_{{\rm{H}}\,{\rm{I}}}^{\mathrm{grad}}$. We therefore find no motivation to introduce an ${N}_{{\rm{H}}{\rm{I}}}$-dependent expectation for the relative orientation between the H i structures and the local magnetic field. Any trend toward an orthogonal alignment between the H i structures and the magnetic field orientation at high ${N}_{{\rm{H}}{\rm{I}}}$ is likely only detectable on scales smaller than the 80' Planck resolution that we have used in this analysis. In contrast, we note that the dependence of ξ on ${{ \mathcal S }}_{{\rm{H}}{\rm{I}}}$ is a significant effect (see Figure 12).

Galactic H i exists in multiple phases, and any given sightline may contain emission from the CNM, the warm neutral medium (WNM), and gas in a thermally unstable phase (e.g., Cox 2005). As discussed above, the H i structures that are most predictive of the magnetic field orientation are predominantly CNM. In fact, the small-scale structure in H i channel maps in general is preferentially cold gas (Clark et al. 2019). By construction, our maps assign polarized intensity to these structures in a given velocity bin, but only total intensity to the more diffuse, interstructure gas. In reality, all H i intensity is probably correlated with polarized dust emission at some level, and it may be that our maps currently overweight the contribution of the CNM to the polarized emission from the neutral medium. Because the WNM-associated H i emission does not strongly contribute to the H i orientation measurements, the empirical correlation between our H i-based maps and the 353 GHz dust polarization constrains the magnetic coherence between the WNM and the CNM.

Phase-decomposed maps of H i emission can be used to explicitly assign different emission properties to different H i phases or velocity ranges, a strategy employed in the dust models of Ghosh et al. (2017) and Adak et al. (2019), and in the foreground cleaning for CIB maps of Lenz et al. (2019). The relationship between the phase structure of H i and the ambient magnetic field is worth further exploration, especially as sophisticated Gaussian fitting algorithms for phase-decomposing hyperspectral data are being developed (e.g., Marchal et al. 2019; Riener et al. 2019).

We establish in this work that regions of high polarization angle dispersion (${ \mathcal S }$) are regions where θ_{H i} is correlated less well with θ₃₅₃ (Figure 12). We find a relationship between ${{ \mathcal S }}_{{\rm{H}}{\rm{I}}}$ and p_{H i} similar to the one found between these quantities in the Planck 353 GHz maps (Figure 11). These observations are qualitatively consistent with a picture in which ${ \mathcal S }$ is sensitive to the angle γ between the mean magnetic field and the plane of the sky (Falceta-Goncalves et al. 2008; Planck Collaboration Int. XX 2015; King et al. 2018). When the magnetic field is along the line of sight, both small perturbations to its 3D orientation, e.g., from interstellar turbulence, and instrumental noise in the Q and U maps induce large changes to its projected 2D orientation, thereby leading to large values of ${ \mathcal S }$. In contrast, when the magnetic field lies mostly in the plane of the sky, perturbations to its 3D orientation have little effect on the projected 2D orientation, leading to small values of ${ \mathcal S }$. Hensley et al. (2019) recently provided empirical evidence for this picture by finding that the dust emission per N_{H i} is positively correlated with ${ \mathcal S }$, as expected from orientation effects.

The dust polarization fraction depends in large part on the value of γ. When the magnetic field is along the line of sight, grain rotation about the field results in no net polarization. When the magnetic field lies in the plane of the sky, the polarization fraction is maximal. In the modeling framework presented here, we do not explicitly account for the effects of γ on the polarized emission. Yet we find that p_{H i} is correlated remarkably well with p₃₅₃ (see Figure 10). To understand this, we note that Q_{H i} and U_{H i} are constructed by integrating R(v, θ) over orientation. In regions where the dispersion in polarization angles is high, this results in strong depolarization. These regions of high ${ \mathcal S }$ are also the ones where, on average, the magnetic field is expected to be oriented more along the line of sight, and the dust polarization fraction is consequently low. On the other hand, regions with low ${ \mathcal S }$ are not significantly depolarized when integrating R(v, θ) over θ, resulting in high polarization fractions. Thus, the effect of γ is implicitly encoded in p_{H i}. Explicitly modeling γ from the H i, perhaps by introducing a factor inversely proportional to ${{ \mathcal S }}_{{\rm{H}}{\rm{I}}}$ to Equations (5) and (6), could improve the correlation between these maps and the 353 GHz polarized dust maps.

On a related note, this work suggests that the velocity structure of H i can be used to decompose the two geometric effects that contribute to depolarization: plane-of-sky depolarization, arising from dispersion of polarization angles within the beam; and line-of-sight depolarization, from multiply oriented dust polarization angles along the line of sight. These effects can, in principle, be isolated and quantified with our 3D Stokes maps.

The present maps are constructed using a fixed RHT window size, meaning that the minimum angular length of H i structures that contribute to θ_{H i} does not vary as a function of velocity. Because a fixed angular size corresponds to different physical sizes for structures at different distances, maps integrated over a large velocity range likely incorporate structures of different physical scales: specifically, more distant structures would need to be physically longer in order to be detected as linear features. This may not be optimal for some use cases, and could in principle be adjusted, given distance information to the gas. We note, however, that any scale-mixing from this effect is likely unimportant at high Galactic latitudes, where the dust is concentrated no farther than a few hundred parsecs away (Capitanio et al. 2017).

We do not find any significant loss of correlation between our maps and the Planck 353 GHz dust Q and U maps as a function of Galactic latitude. This is perhaps surprising, as the previous work that motivated our magnetically coherent mapping (Section 3.1) was mostly focused on high Galactic latitudes. Evidently, H i orientation remains a reliable proxy for magnetic field orientation along sightlines closer to the Galactic plane, although it would be interesting to investigate whether the structures that dominate this measurement are substantially different physically from the "H i fibers" (Clark et al. 2014) that pattern the diffuse sky.

In addition to the study of the ISM, our model has significant relevance to cosmology. One of the major goals of modern observational cosmology is the search for primordial B-mode polarization in the CMB, a signature of inflation-era gravitational waves (Kamionkowski et al. 1997; Seljak & Zaldarriaga 1997). A paramount obstacle to successful detection is the polarized foreground signal from Galactic dust emission, which at CMB frequencies is several orders of magnitude brighter in polarization than the sought-after primordial signal (Dunkley et al. 2009; Flauger et al. 2014; BICEP2/Keck and Planck Collaborations et al. 2015; Planck Collaboration XI 2018).

Realistic models of Galactic polarized dust emission are therefore a crucial ingredient for experiment design and forecasting, component separation, and validation of any claimed detection. Polarized CMB foregrounds arise from the three-dimensional ISM, which is shaped by nonlinear MHD processes that are fundamentally non-Gaussian, particularly on smaller scales (e.g., Burkhart et al. 2009; Allys et al. 2019). While data-driven sky models like the Planck Sky Model (Delabrouille et al. 2013) or the Python Sky Model (Thorne et al. 2017) are able to incorporate this non-Gaussianity on scales that have already been measured well, they typically resort to simple Gaussian parameterizations at small scales.

A number of approaches to provide this small-scale information are being actively pursued. MHD simulations have successfully reproduced many of the observed properties of the polarized dust emission (Kritsuk et al. 2018; Kim et al. 2019), and if pushed to sufficiently high resolution, can generate small-scale, non-Gaussian spatial fluctuations in a physically motivated way. However, MHD simulations are limited in being able to reproduce only the statistical properties of the Galaxy, not the specific morphology of the Galaxy itself. High-resolution ancillary data like Galactic H i also naturally contain physical, non-Gaussian spatial structure while also corresponding directly to the particular structure of our Galaxy. Past efforts have largely focused on two-dimensional data with a parametric, phenomenological perscription for the line-of-sight structure (e.g., Planck Collaboration Int. XLIV 2016; Ghosh et al. 2017; Vansyngel et al. 2017; Adak et al. 2019). While simple, these models have been effective in reproducing many of the spatial statistics of the polarized dust emission, such as the E/B asymmetry and polarization cross spectra. Rather than employ a phenomenological prescription for line-of-sight structure, our maps incorporate the three-dimensional distribution of magnetic coherence entirely from data.

The maps presented in this work also provide a data-driven framework for modeling spatial variability in the frequency-dependent dust emission. The superposition of multiple emitting regions along the line of sight, each with their own dust properties (such as temperature) and magnetic field orientation, hinders the ability of a map of dust emission at one frequency to be extrapolated to another frequency, i.e., "frequency decorrelation" (Tassis & Pavlidou 2015). While currently not detected in the Planck dust polarization data (Planck Collaboration XI 2018), even modest levels of decorrelation can impact B-mode science (Poh & Dodelson 2017; Hensley & Bull 2018).

Line-of-sight variability of the dust emission has been modeled both with simple parametric prescriptions that reproduce the 2D observations as well as directly incorporating 3D data sets such as dust extinction maps (Poh & Dodelson 2017; Martínez-Solaeche et al. 2018). Our maps can improve considerably on these methods by incorporating not only line-of-sight information of the density distribution, but also the magnetic field orientation. Mapping different dust emission models onto the three-dimensional maps derived in this work will provide new data-driven models and predictions for the expected levels of frequency decorrelation.

Finally, we note that even as-is, the H i-based Stokes maps are an ideal template for assessing the residual contamination in a foreground-cleaned CMB observation (e.g., Thorne et al. 2019). Any cross-correlation signal between such maps and our model is certain to arise from the Galactic ISM. Because we restrict the H i data to velocities $| v| \lt 90$ km s⁻¹, our H i-based maps are free from CIB contamination (Chiang & Ménard 2019). Further, as an entirely independent data set, the H i-derived maps will not contain correlated systematics with any CMB observations.

The global structure of the GMF is an outstanding problem in astrophysics, with wide-ranging consequences for understanding cosmic ray propagation, dynamo theory, star formation, and galactic evolution (Haverkorn 2015; Jaffe 2019). At present, the state-of-the-art data-driven models for the GMF are parametric fits to a few observables, typically some subset of extragalactic and pulsar rotation measures (RMs), synchrotron emission, and polarized dust emission (e.g., Jansson & Farrar 2012; Jaffe et al. 2013; Terral & Ferrière 2017). Each model relies on a set of assumptions and contains large uncertainties, such that there is a wide variation in the GMF structure obtained, even between models based on the same observational data. New observational constraints and more robust statistical frameworks are two avenues for improvement (Boulanger et al. 2018; Haverkorn et al. 2019). The H i-derived Stokes maps constitute a new observational constraint, and here we detail a few specific ideas for their use in GMF modeling.

The influence of nearby (≲few hundred pc) ISM structure on observational tracers presents a major challenge for GMF models. Our proximity to these local features means that they may be observed on large angular scales, but may be unimportant to the Galactic-scale field structure. For instance, the Sun sits in the Local Bubble (also called the Local cavity), an underdense region of the nearby ISM thought to have been carved out by a series of supernovae (Cox & Reynolds 1987; Lallement et al. 2003). The expansion of the Bubble likely distorted the local GMF, measurably impacting the observed dust emission, and yet the Local Bubble is not included in GMF models (Alves et al. 2018). The 3D Stokes maps presented here can be used to disentangle the influence of these local features from the rest of the emission in observations like the Planck 353 GHz map. Our maps may be particularly interesting to compare with detailed 3D maps of the local ISM (Lallement et al. 2014; Capitanio et al. 2017), as well as with the Alves et al. (2018) model for the magnetic structure of the Local Bubble.

Faraday rotation is an important and widely used probe of the GMF (e.g., Hutschenreuter & Enßlin 2019). RM is an integral quantity, proportional to the product of the thermal electron density (n_e) and the line-of-sight component of the magnetic field integrated between the observer and the source. The largest uncertainty in estimating the magnetic field strength from measurements of RM is thus the distribution of n_e along the line of sight (Beck et al. 2003). A better understanding of the magnetic field and its relationship to gas density between phases of the ISM is crucial to making progress (Boulanger et al. 2018). If our maps primarily trace the CNM, they may be useful for mapping the three-dimensional phase distribution of the ISM, or for identifying sharp changes along the line of sight in either n_e or the magnetic field, in the vein of Van Eck et al. (2017).

Ultimately, we see this work as an important step toward the goal of mapping the three-dimensional structure of the GMF in three spatial dimensions. Our 3D Stokes maps constitute a novel measurement: a local probe of the magnetic field orientation as a function of velocity. Other measurements of the magnetic field orientation in the neutral ISM are not local,⁷ but integrated: polarized dust emission traces the field integrated along the entire line of sight, and starlight polarization traces the field between the observer and the star.

The next step is to map these data to distance. This calls for the synthesis of maps like those presented in this work with other data sets. Stellar reddening and other measurements enable maps of the dust distribution in three spatial dimensions (e.g., Marshall et al. 2006; Green et al. 2019). Large-scale starlight polarization surveys like PASIPHAE at high Galactic latitudes (Tassis et al. 2018) will soon provide optical polarimetry toward stars with distances measured by Gaia (Gaia Collaboration et al. 2016), improving the density of starlight polarization measurements at high Galactic latitude by orders of magnitude over existing catalogs (e.g., Heiles 2000; Berdyugin et al. 2014). These data sets can be combined with our velocity-resolved maps to produce maps in four dimensions—three spatial dimensions and radial velocity. Some of these data have already been used to map the four-dimensional density distribution in the Galactic plane (Tchernyshyov & Peek 2017).