The Physical Nature of Neutral Hydrogen Intensity Structure

TL;DR

This study demonstrates that small-scale structures in neutral hydrogen emission maps are actual density features, primarily cold neutral medium, aligned with magnetic fields, and not artifacts of turbulent velocity fluctuations.

Contribution

It provides evidence that HI intensity structures are real density features, challenging the velocity caustics hypothesis, and links small-scale HI structures to the cold neutral medium phase.

Findings

HI structures are real density features, not velocity caustics.

Small-scale HI structures are associated with colder, denser gas.

The shallow power spectrum of narrow HI channels is due to CNM properties.

Abstract

We investigate the physical properties of structures seen in channel map observations of 21-cm neutral hydrogen (HI) emission. HI intensity maps display prominent linear structures that are well aligned with the ambient magnetic field in the diffuse interstellar medium (ISM). Some literature hold that these structures are "velocity caustics", fluctuations imprinted by the turbulent velocity field, and are not three-dimensional density structures in the ISM. We test this hypothesis by stacking probes of the density field -- broadband far infrared (FIR) emission and the integrated HI column density (NHI) -- at the locations of linear HI intensity structures. We find that the HI intensity features are real density structures and not velocity caustics. We generalize the investigation to all small-scale structure in HI channel maps, and analyze this correlation as a function of velocity channel width, finding no measurable contribution from velocity caustics to the HI channel map emission. Further, we find that small-scale HI channel maps structures have elevated FIR/NHI, implying that this emission originates from a colder, denser phase of the ISM than the surrounding material. The data are consistent with a multi-phase diffuse ISM in which small-scale structures in narrow HI channel maps are preferentially cold neutral medium (CNM) that is anisotropically distributed and aligned with the local magnetic field. The shallow spatial power spectrum (SPS) of narrow HI channels is often attributed to velocity caustics. We conjecture instead that the small-scale structure and narrow linewidths typical of CNM explain the observed relationship between the SPS and channel width.

Equations

$$I_{\mathrm{H}\text{sc}\mathrm{I}}(\vec{x})={\int }_{v_1}^{v_2}{T}_b(\vec{x},v)\delta v,$$

$${\bar{I}}_{857}^w=\frac{{\sum }_{i=1}^n{I}_i{w}_i^{\delta v}}{{\sum }_{i=1}^n{w}_i^{\delta v}},$$

$$\Delta I_{857}={\bar{I}}_{857}^w-{\bar{I}}_{857}.$$

Full text

dust, extinction — ISM: magnetic fields — ISM: structure — magnetic fields — radio lines: ISM — turbulence

††software: astropy (Astropy Collaboration et al., 2013, 2018), cmocean (Thyng et al., 2016), matplotlib (Hunter, 2007), numpy (Oliphant, 2015)

1 Introduction

The diffuse interstellar medium (ISM) fills most of the volume of the Milky Way. This material, the progenitor of dense structures that eventually form stars, is multi-phase, magnetized, and turbulent. The ISM is broadly partitioned into ionized, molecular, and atomic components, with the atomic gas composed of a cold neutral medium (CNM), warm neutral medium (WNM), and a thermally unstable component (for reviews, see Ferrière, 2001; Cox, 2005). Understanding the distribution of matter and energy within and between these phases is a major goal of ISM research. Since its initial detection in 1951 (Ewen & Purcell, 1951; Muller & Oort, 1951), the 21-cm line from the hyperfine transition of neutral hydrogen (Hi) has been an extremely fruitful probe of the distribution of interstellar gas (e.g. Kalberla & Kerp, 2009).

The ISM occupies three spatial dimensions, but our direct observations are limited to the two-dimensional plane of the sky. A third dimension is available to spectroscopic observations, which record intensity as a function of frequency. Emission line observations like 21-cm Hi measurements are typically reported in position-position-velocity (PPV) space, where the third dimension represents the line-of-sight velocity implied by the frequency shift from the line rest frequency. Such shifts are not caused exclusively by bulk gas motions, but also by thermal and turbulent line broadening. This PPV information is therefore distinct from position-position-position (PPP) space, but is also not merely a projection of the six-dimensional spatio-kinematic information that describes flows of interstellar material. This observational restriction complicates the interpretation of ISM structures, particularly the determination of how observed PPV intensity structures relate to “real” PPP density structures (Beaumont et al., 2013). This is important for understanding molecular clouds, which are commonly probed by PPV molecular line observations (Goodman, 2011).

In simple flows, there exists a one-to-one mapping between PPP and PPV, e.g. in a pure Hubble flow. At least in theory, it is possible to construct a circumstance in which the structure seen in PPV space can be dominated by the velocity field, and thus not a simply-interpretable probe of the density field. This can be understood by considering the extreme case: incompressible, constant-density turbulence. Such an idealized flow will lack any PPP density structure, but can still exhibit PPV intensity structure. A featureless incompressible flow is far from a realistic description of the ISM (e.g. Elmegreen & Scalo, 2004). Still, it is reasonable to interrogate to what extent intensity structures in Hi PPV observations are influenced by the true density or velocity fields.

Recent high-dynamic range observations of the diffuse ISM have renewed interest in this question. Hi surveys such as the Galactic Arecibo $L$-Band Feed Array Survey (GALFA-Hi; Peek et al., 2018) and Hi4PI (HI4PI Collaboration et al., 2016) show a wealth of linear structure, particularly prominent at high Galactic latitudes and in narrow ($\sim$few km ${\rm s^{-1}}$) velocity channels. Clark et al. (2014) developed a machine vision algorithm called the Rolling Hough Transform (RHT) to measure the linearity of image structure, and used the RHT to show that Hi intensity structure is well aligned with the ambient magnetic field as traced by optical starlight polarization. With the advent of all-sky maps of polarized dust emission at 353 GHz (Ade et al., 2015), Clark et al. (2015) measured a robust, tight correlation between linear Hi intensity structures and the magnetic field. Similar structures have been seen in other Hi emission data observed at high angular resolution (Martin et al., 2015; Kalberla et al., 2016, 2017; Blagrave et al., 2017). These Hi intensity observations are reminiscent of linear, magnetically aligned structures in other media, including Hi absorption (McClure-Griffiths et al., 2006) and Faraday depth structures in low-frequency radio polarimetric data (Jelić et al., 2015).

Do these results indicate that the true (PPP) density field is striated, and that these striations are preferentially elongated along the ambient magnetic field? Or is this an effect of the velocity field? A number of explanations have been proposed for the origin of this magnetically aligned structure. Most theoretical explanations have focused on physical mechanisms for a correlation between the orientation of elongated density structures and the magnetic field (e.g. Hennebelle & André, 2013; Inoue & Inutsuka, 2016; Soler & Hennebelle, 2017).

However, the striking magnetic alignment of Hi structures in velocity channel maps (Clark et al., 2014, 2015) could in theory be an imprint of the velocity field. This interpretation has been put forward by Lazarian & Yuen (2018), based on the theory of Lazarian & Pogosyan (2000): that the structures seen in these Hi channel maps are “velocity caustics”: intensity structures created by turbulent velocity fluctuations coincident at a given $v_{lsr}$. In this picture, the measured alignment between the magnetic field orientation and the orientation of Hi intensity structures has little or no relation to the underlying density field, but is instead caused by the correlation between velocity and magnetic field orientation that is predicted from magnetohydrodynamic (MHD) turbulence (Goldreich & Sridhar, 1995; Elmegreen & Scalo, 2004). Shearing turbulent eddies that are elongated in the direction of the local magnetic field have a velocity gradient that is steepest in the direction perpendicular to the magnetic field. If Hi intensity structures are dominated by velocity caustics, the correlation between the gradient of the intensity in thin velocity channel maps and the dust polarization angle would be a measurement of the velocity-magnetic field alignment (Lazarian & Yuen, 2018).

In this work, we apply a direct test to determine whether the linear structures identified in narrow spectral channels are velocity caustics or real density structures. We measure the correlation between linear channel map Hi structures and the density field. We further generalize our investigation to all small-scale structure in thin-channel Hi intensity maps. This paper is structured as follows. In Section 2 we detail the data products and numerical simulations used in this work. In Section 3 we test whether the structure seen in Hi intensity maps is consistent with being velocity caustics, and investigate the physical nature of structures observed in Hi intensity. In Section 4 we discuss how our findings fit into a broader picture of the neutral ISM, and suggest several testable predictions. In Section 5 we summarize and conclude.

2 Data and Simulations

In this work we analyze 21-cm line observations of Galactic Hi and measurements of the far infrared (FIR) intensity at 857 GHz. We make use of several publicly available Planck legacy data products. Much of this work uses the 857 GHz frequency map from the Planck collaboration data release R3.01 (Planck Collaboration et al., 2018). We also use the High Frequency Instrument (HFI) point source mask to exclude regions of the data that are identified as point sources in the FIR (Planck Collaboration et al., 2016b). The zero-point of the Planck 857 GHz emission is significantly uncertain, and the data must be corrected for the cosmic infrared background monopole. In this work we apply a monopole offset correction of $0.64$ MJy/sr (Planck Collaboration et al., 2016a). However, we intentionally design the experiments in this paper such that our results are not qualitatively dependent on this value.

We use the publicly available Data Release 2 (DR2) of GALFA-Hi (Peek et al., 2018). GALFA-Hi is a large-area ($\sim 4$ sr) survey with high spatial ($4^{\prime}$) and spectral ($0.18$ km ${\rm s^{-1}}$) resolution. We use three of the data products provided in DR2: PPV cubes of the Hi brightness temperature ($T_{b}(v)$), the RHT output of GALFA-Hi velocity channel map data, and the stray radiation-corrected column density map ($\mathrm{N}_{\mathrm{H{\sc I}}}$), derived from $T_{b}(v)$ integrated over $\left|v\right|<90$ km ${\rm s^{-1}}$. While Hi is susceptible to self absorption effects in very cold, dense media, this phase is relatively uncommon in the low-density ISM under investigation here (Murray et al., 2018a). In this work we ignore this effect and treat velocity-integrated brightness temperature and column density interchangeably. In what follows, an Hi “channel map” refers to GALFA-Hi $T_{b}(v)$ integrated over a narrow range in velocity:

[TABLE]

where $\vec{x}$ represents the plane-of-sky spatial dimensions. The velocity “channel width” is $\delta v=v_{2}-v_{1}$, and the channel is centered on $v_{0}=(v1+v2)/2$. $\mathrm{N}_{\mathrm{H{\sc I}}}~{}$is calculated via Equation 1 with $v_{1}=-90$ km ${\rm s^{-1}}$ and $v_{2}=+90$ km ${\rm s^{-1}}$, and reported in units of $\mathrm{cm}^{-2}$.

Unless otherwise noted, the analysis in this work is carried out on a large region of the GALFA-Hi sky. We select data in the GALFA-Hi footprint at Galactic latitudes $|b|>30^{\circ}$. We exclude sightlines in the HFI point source mask, as well as those with \mathrm{N}_{\mathrm{H{\sc I}}}~{}$$>8\times 10^{20}~{}\mathrm{cm}^{-2}, to focus on the diffuse ISM. The resulting sky area considered comprises $\sim 50\%$ of the GALFA-Hi sky, or $\sim 2$ sr. A portion of the Planck 857 GHz emission and the GALFA-Hi channel map data is shown in Figure 1.

Finally, to illustrate the relationship between the channel map structure and the column density structure of turbulent flows, we run a series of simple isothermal simulations. These are three-dimensional hydrodynamic simulations produced with RAMSES (Teyssier, 2002). We run simulations at three sonic Mach numbers ($\mathcal{M}_{s}$=0.5, 1, 5), each with size 1283, and translate these to PPV space with a velocity resolution $\delta v_{0}=0.1$ km ${\rm s^{-1}}$.

3 The physical nature of thin-channel Hi intensity structures

Broadband observations of thermal dust emission are sensitive to the density, but not velocity, field. No velocity caustics exist in dust emission maps. The observed correlation between the orientation of filamentary structures in dust emission and the dust polarization field (Planck Collaboration et al., 2016c; Planck Collaboration Int. XXXV et al., 2016; Malinen et al., 2016; Soler et al., 2017) is thus unequivocally a correlation between density structures and the projected magnetic field. $\mathrm{N}_{\mathrm{H{\sc I}}}~{}$is likewise sensitive to the density field, but not the velocity field, because of the integration over a broad velocity range (see Section 2).

3.1 Theoretical considerations

Before undertaking any data analysis, we ask under what physical conditions we might expect the thin-slice Hi intensity to be dominated by velocity caustics. Here we should emphasize an important distinction between the velocity caustics picture and other questions of PPV structure. Velocity caustics occur when the dominant intensity contours – the spatial structures – in a velocity channel are imprinted by the turbulent velocity field rather than the density field. The absence of velocity caustics does not imply that the distribution of emission along the line of sight is unaffected by the radial velocity. The classic example of an expanding shell, for instance, will cause emission from a single PPP structure to contribute to multiple line-of-sight velocities. This type of velocity-induced emission structure, and its inverse, whereby emission at one velocity originates from multiple density structures along the line of sight, certainly complicates the relationship between PPP and PPV space (see Figure 1 of Beaumont et al., 2013). The difficulty of reverse-engineering PPP density structures from PPV observations is a long-standing problem (e.g. Ostriker et al., 2001; Clarke et al., 2018). The velocity caustics picture tested here is a distinct and specific claim: that the intensity structure in narrow velocity channels is an image of the turbulent velocity field. Put another way, the intensity structure observed at line-of-sight velocity $v_{0}$ converges to the constant-density turbulent velocity field at $v_{0}$ as the channel width $\delta v$ becomes small enough. This turbulent velocity field can be uncorrelated with the density field, as explicitly assumed in the Lazarian & Pogosyan (2000) treatment of incompressible turbulence.

A significant velocity caustics contribution to synthetic channel map emission has been reported in a number of simulation studies. These results are often interpreted as meaning that channel map observations of Hi contain little to no structure that can be associated with density structures. The relevance of this picture for PPV observations of Hi is partly a question of how well these simulations capture the salient ISM physics, and partly a question of whether the transformation from the simulated PPP space to mock PPV observations correctly mimics the process by which 21-cm emission maps to line-of-sight velocity. Unphysical transformations between PPP and PPV space can generate misleading effects.

One early exploration relating the morphology of synthetic channel maps of MHD simulations to the underlying velocity and density fields was undertaken by Pichardo et al. (2000). The authors compared the structure of slices through the PPP velocity field along the line of sight ($v_{z}$) to synthetic channel map data. Pichardo et al. (2000) concluded that their channel maps were more morphologically similar to the $v_{z}$ slices than to slices through the PPP density field, and that the channel maps exhibited more small-scale structure than the density slices. This result was in qualitative agreement with the contemporary work of Lazarian & Pogosyan (2000).

However, as the authors note, Pichardo et al. (2000) neglect thermal broadening in their computation of channel map structure. Translating PPP simulations into PPV space requires the inclusion of thermal line broadening, which is a convolution along the velocity axis by a Gaussian of $\sigma=\sqrt{k_{B}T/m}$, where $T$ is the gas temperature, $m$ is the particle mass, and $k_{B}$ is the Boltzmann constant. Thermal line broadening is not a radiative transfer effect, but rather a fundamental property of thermodynamics. If this temperature-dependent convolution along the velocity axis is ignored, the influence of the turbulent velocity on channel map structures is artificially inflated. The same omission of thermal broadening is made in Lazarian et al. (2001), another work that finds prominent velocity caustics in simulations and concludes that features in observational channel maps must also be velocity structures. Miville-Deschênes et al. (2003) noted the severe limitations that thermal broadening places on the feasibility of using the Lazarian & Pogosyan (2000) method to derive properties of the turbulent velocity field from the structure of narrow channel maps.

The velocity caustics picture and the effects of thermal broadening are illustrated in Figure 2. Velocity caustics manifest as channel map structures induced by the three-dimensional velocity field, but weakly or not at all correlated with the three-dimensional density field (Figure 2, panel C2). This will not be observable if thermal broadening washes out velocity fluctuations (panel D2) or if three-dimensional density structures dominate the emission contours (panels C4, D4). Omitting thermal broadening from synthetic observations is particularly pernicious because it implies that channel maps directly measure the mass distribution as a function of radial velocity (panels B2 and B4). This is, in reality, unobservable, and thus cannot be compared to any real data.

Burkhart et al. (2013) examined the effect of the sonic Mach number ($\mathcal{M}_{s}$) on the relationship between structures in PPP and PPV. High $\mathcal{M}_{s}$ flows can contain shocks, which sweep up matter and thus strongly correlate the turbulent velocity field with the density field. Low $\mathcal{M}_{s}$ flows do not strongly correlate the velocity and density structures. However, $\mathcal{M}_{s}$ is the ratio of the turbulent velocity to the sound speed in the medium, which is set by the kinetic temperature. Burkhart et al. (2013) neglected thermal broadening in much of their analysis, and so the PPV cubes they analyze are not the PPV-space representation of turbulent flows, but instead are histograms of the turbulent velocity. Assuming that these are the PPV-space representation of flows with a given $\mathcal{M}_{s}$ leads to the conclusion that PPP and PPV space are more closely related for higher $\mathcal{M}_{s}$ flows, summarized in the abstract of Burkhart et al. (2013) as “the dominant emission contours in PPP and PPV are related for supersonic gas but not for subsonic”.

Missing from this analysis is the fact that low $\mathcal{M}_{s}$ flows have higher temperatures relative to their turbulent velocities, and thus broader emission profiles. PPV cubes of low $\mathcal{M}_{s}$ flows that correctly incorporate the thermal line broadening have features that are smooth in velocity space and highly correlated with the PPP density structure. These lack velocity caustics: any given velocity slice of a low $\mathcal{M}_{s}$ PPV cube is highly correlated with the integrated line emission, especially for lightweight carriers such as Hi where this effect is strong. Thus, while PPP and PPV structures are correlated dynamically for high $\mathcal{M}_{s}$ flows, they are also correlated thermodynamically for low $\mathcal{M}_{s}$ flows.

We demonstrate this in Figure 3 using the isothermal turbulent simulations described in Section 2. These are hydrodynamic simulations that are not representative of the ISM, particularly for their lack of magnetic fields and thermal structure. This is intentional, as we illustrate here a simple point about $\mathcal{M}_{s}$ and thermal broadening, that holds irrespective of those physics or of the particular simulation scheme used. Figure 3 shows narrow channel maps of each synthetic 21-cm PPV cube for $v_{0}=0$, along with its corresponding column density map. In each simulation, the channel map structure is well correlated with the integrated intensity (column density) structure, but the reason for this correlation changes as function of $\mathcal{M}_{s}$. The supersonic simulation has correlated density and velocity fields that cause structures in the channel map to be well correlated with column density structures. The subsonic simulation has channel map structures that are correlated with the column density by thermal broadening. Between the dynamic and thermodynamic effects, none of the synthetic channel maps display prominent velocity caustics.

If 21-cm PPV data contains prominent velocity caustics, then, our thermodynamic expectation is that these structures would have to be produced by the cold, supersonic gas in the diffuse Hi. As noted, the expectation from theoretical studies of turbulence is that this medium will have correlated three-dimensional density and velocity fields. One might therefore not expect that intensity channels show significant velocity-induced structure that is uncorrelated with the density field. Still, the idea that thin Hi channel maps are dominantly structured by the turbulent velocity field is pervasive in the literature. We have devised a direct observational test of this picture.

3.2 Magnetically aligned Hi intensity structures are not velocity caustics

We conduct an experiment to test whether thin-channel Hi intensity structures are an imprint of the velocity field. We stack FIR emission on the locations of the Hi intensity structures. As noted above, FIR emission probes the total density, but does not at all trace the velocity field. If the Hi intensity structures are velocity caustics, there should be no enhancement in the FIR emission strength at the location of the intensity structures, relative to the surrounding medium. Conversely, an enhancement in the FIR emission indicates that the structures are true density features.

To highlight the linear structures in Hi emission we use the RHT (Clark et al., 2014). The RHT is a machine vision algorithm that quantifies linearity as a function of orientation, and has been used to trace the orientation of Hi emission structures that are prominently aligned with the local magnetic field (Clark et al., 2014, 2015; Jelić et al., 2018), as well as to measure linearity in other data. The RHT output for $I_{v_{0}}(\vec{x})$, an Hi channel map centered at velocity $v_{0}$, is $R_{v_{0}}(\vec{x},\theta)$, the linear intensity as a function of orientation ($\theta$) for every image pixel ($\vec{x}$). Clark et al. (2014) referred to the narrow Hi channel map structures that are well aligned with the magnetic field as “Hi fibers”.

We measure the correlation between structures in RHT intensity and in the $I_{857}$ and $\mathrm{N}_{\mathrm{H{\sc I}}}~{}$maps, by stacking $I_{857}$ and $\mathrm{N}_{\mathrm{H{\sc I}}}~{}$on $R_{v_{0}}(\vec{x},\theta)$. Here, $R_{v_{0}}(\vec{x},\theta)$ is the RHT intensity for a channel map with $\delta v=2.94$ km ${\rm s^{-1}}$, centered at $v_{0}\sim 1.5$ km ${\rm s^{-1}}$. We stack a $101^{\prime}\times 101^{\prime}$ region of sky around each nonzero pixel in the RHT output. We perform a separate stack on the RHT intensity maps at each orientation angle. We weight each of these sky regions by the $R_{v_{0}}(\vec{x},\theta)$ value of the central pixel. The result of this process, for each stacked map, is a $101^{\prime}\times 101^{\prime}$ stacked map for each orientation bin of the RHT output. We collapse the orientation information by rotating each stacked map by its orientation, such that each map is aligned with $\theta=0^{\circ}$. The GALFA-Hi data contains telescope scan artifacts at specific orientations because of the survey’s “meridian nodding” or “basketweave” scan strategy (see Peek et al., 2011b). We try cutting orientations near $180^{\circ}$ that are known to be most contaminated by these artifacts, and verify that even conservative cuts do not meaningfully alter our results.

Figure 4 shows the result of this stacking experiment. We emphasize that if the thin-channel Hi intensity structures were velocity caustics, they would have no correlation with the total density tracers, and neither of these stacks would show an enhancement at the location of the $R_{v_{0}}(\vec{x},\theta)$ emission. Clearly, there is a significant enhancement both in 857 GHz emission and in $\mathrm{N}_{\mathrm{H{\sc I}}}~{}$at the location of the Hi intensity structures. A description of the magnetically aligned Hi fibers as velocity caustics is inconsistent with this result.

3.3 Thin-channel Hi intensity structures are not velocity caustics

In Section 3.2 we showed that elongated linear structures in thin-channel Hi intensity images correspond to dust-bearing density structures in the ISM. We made this measurement using the RHT, but other methods have been used to study the intensity structure of thin-channel Hi observations, such as the spatial gradient (e.g. Lazarian & Yuen, 2018). Could it be that the structures that are highlighted by the RHT are special in terms of their physical nature – that is, are Hi intensity features with strong RHT signatures more likely than other features to be real density structures?

We test this by applying a different algorithm to the thin-channel Hi intensity data: the unsharp mask (USM). The USM is a high-pass filter, an operation that filters out low-frequency structure. Many popular algorithms for quantifying ISM spatial structure share this high-pass filter property, including the spatial gradient (Gaensler et al., 2012; Planck Collaboration Int. XXXV et al., 2016) and the Hessian (e.g. Planck Collaboration et al., 2016c). The USM is used on its own to study ISM structure (e.g. Lee et al., 2009; Kalberla et al., 2016), and is also the first step in the RHT. Applying the USM to the Hi channel map data highlights the small-scale intensity structure.

We repeat the experiment of Section 3.2, but stack the FIR and $\mathrm{N}_{\mathrm{H{\sc I}}}~{}$data on the USM of thin-channel Hi intensity rather than the RHT output. We find a similar result to the one summarized in Figure 4: there is a clear FIR and $\mathrm{N}_{\mathrm{H{\sc I}}}~{}$excess at the location of the stack. In Figure 5 we show small regions of Hi velocity channel data and their corresponding USM intensity. The Hi channel maps shown are all centered at $v_{lsr}\sim 0.4$ km ${\rm s^{-1}}$, but vary in channel width from $0.18$ km ${\rm s^{-1}}$ to $7.9$ km ${\rm s^{-1}}$. Beneath each data snapshot in Figure 5 we show the corresponding stack of $I_{857}$ on the USM over the region described in Section 2. Each shows a clear excess of FIR emission at the center of the stack. Thin-channel Hi intensity structures are not velocity caustics.

In the picture set forth by Lazarian & Pogosyan (2000), the degree to which the intensity structure is an imprint of the turbulent velocity field depends on the velocity channel width. The claim made by Lazarian & Yuen (2018), that intensity structures in channels of width $\delta v\sim 3$ km ${\rm s^{-1}}$ are velocity caustics, is already ruled out. Indeed, the correlation with the FIR shown in Figure 5 reveals that the structures in Hi channel maps are correlated with the density field down to the narrowest measured channel map width.

In Figure 6 we show the results of an experiment to test the Lazarian & Pogosyan (2000) prediction that narrower Hi velocity channels will be more dominated by velocity caustics. We compute the mean of the Planck 857 GHz data, $\bar{I}_{857}$. This mean is computed over the region of diffuse, high Galactic latitude sky described in Section 2. We then compute the weighted mean of the $I_{857}$ data, where the weights are the USM of the Hi intensity channel data. That is, we compute

[TABLE]

where $I_{i}$ is $I_{857}$ in a given pixel, and $w_{\delta v}$ is the USM of an Hi intensity channel with velocity width $\delta v$.

We then compute

[TABLE]

The quantity $\Delta I_{857}$ represents the average FIR enhancement correlated with small-scale Hi intensity channel structure. If the Hi intensity structure is caused entirely by velocity caustics, $w^{\delta v}$ will be uncorrelated with the density field, and $\Delta I_{857}$ will equal 0.

As velocity caustics are expected to contribute more strongly to the Hi emission in narrower velocity channels, the caustics picture predicts a decrease of $\Delta I_{857}$ as $\delta v$ decreases. This caustics expectation is illustrated in Figure 6 by a pale green wedge. While the slope and amplitude of the expected decline depends on details of the ISM density and velocity fields, the predicted trend is toward a more prominent contribution of the turbulent velocity field to the Hi intensity in narrower velocity channels. The gray dashed vertical line in Figure 6 indicates the velocity channel width used in the Clark et al. (2014, 2015) papers and in the gradient-based analyses of Lazarian & Yuen (2018) and Lazarian et al. (2018), to reflect the claim in the latter papers that structures in Hi intensity at these channel widths are velocity caustics.

We calculate $\Delta I_{857}$ down to the narrowest $\delta v$ possible: the GALFA-Hi velocity resolution of $0.18$ km ${\rm s^{-1}}$, well below the full width at half maximum (FWHM) linewidth of the coldest known Hi structures (Peek et al., 2011a). Even with $\delta v=0.18$ km ${\rm s^{-1}}$, $\Delta I_{857}$ is significantly inconsistent with [math]. This is a qualitatively identical result to the stacking experiment shown in Figure 5: there is a significant $I_{857}$ GHz excess at the location of the USM structures for all $\delta v$. We find the same trend if we instead measure \Delta$$\mathrm{N}_{\mathrm{H{\sc I}}}, the average $\mathrm{N}_{\mathrm{H{\sc I}}}~{}$enhancement correlated with small-scale Hi channel map structure, computed in analogy to Equations 2 and 3. The channel maps used in this analysis are all centered at $v=0.4$ km ${\rm s^{-1}}$, but the measured trend is qualitatively unchanged if we instead center the channel maps near the peak of the Hi line integrated over our region, $v\sim-3.3$ km ${\rm s^{-1}}$.

To estimate the variance in $\Delta I_{857}$, we block bootstrap the data by dividing the selected sky area into 40 regions of roughly equal numbers of pixels. We sample these regions with replacement $10^{7}$ times and compute the $16^{\mathrm{th}}$ and $84^{\mathrm{th}}$ percentile ($1\sigma$) $\Delta I_{857}$ values from this distribution. These are visualized as the light gray error bars in Figure 6. As error bars, these are conservative and covariant, because they represent the sky variation rather than a statistical error.

Small-scale noise and telescope scan artifacts, more prominent in narrow velocity channels, will be uncorrelated with $I_{857}$ and will drive $\Delta I_{857}$ toward 0. To measure the effect of noise as a function of velocity channel width, we repeat the experiment shown in Figure 6, but prior to taking the USM of each velocity channel we add a “pure noise” component. This noise map is the $\delta v=0.72$ km ${\rm s^{-1}}$ GALFA-Hi emission centered at $v_{lsr}=400$ km ${\rm s^{-1}}$, multiplied by a factor $1/\sqrt{N}$, where $N=\delta_{v}/0.72$. We repeat the same experiment multiplying the noise component by $2$. The results are shown in Figure 7. The expected noise contribution has the same shape as the slight downturn in $\Delta I_{857}$ for the narrowest velocity channels in Figure 6. The amplitude of the slight downturn is likewise consistent with being solely caused by noise. Hi channel map noise similarly affects the stacked FIR emission in the narrowest velocity channels of Figure 5.

Neither $\Delta I_{857}=0$ nor a decline in $\Delta I_{857}$ with decreasing velocity channel width is consistent with the data. The FIR associated with small-scale channel map structure is significantly higher than the sky average for all velocity channel widths, and remarkably flat as a function of channel width. Only for the narrowest channel widths, $\delta v<1$km ${\rm s^{-1}}$, do we measure even a marginal decrease in $\Delta I_{857}$. This $\sim 15\%$ decrement in $\Delta I_{857}$ is consistent with noise in the Hi intensity maps. Evidently, not only are thin-channel Hi intensity structures not dominated by velocity caustics, but there is no measurable contribution to the Hi intensity from velocity caustics at all.

3.4 Thin-channel Hi intensity structures are preferentially CNM

We have so far established that high-spatial frequency structures in narrow Hi velocity channels are density structures, not velocity caustics. Here we investigate the physical nature of these density structures, and suggest a physical picture of the ISM that is consistent with these observations.

We examine the distribution of FIR/$\mathrm{N}_{\mathrm{H{\sc I}}}~{}$for structures at different spatial scales. Figure 8 shows the (normalized) histogram of $I_{857}$/$\mathrm{N}_{\mathrm{H{\sc I}}}~{}$as a function of the intensity of the USM ($w^{\delta v}$) of an Hi intensity channel with width $\delta v=2.2$ km ${\rm s^{-1}}$. Higher values in the USM correspond to channel map features that are higher contrast with respect to the local low-frequency spatial structure. The gray histogram shows the $I_{857}$/$\mathrm{N}_{\mathrm{H{\sc I}}}~{}$for the entire region of sky considered. The histogram of $I_{857}$/$\mathrm{N}_{\mathrm{H{\sc I}}}~{}$clearly shifts toward higher values as sky pixels with higher USM values are selected. The high spatial frequency structures that the USM highlights have preferentially higher $I_{857}$/$\mathrm{N}_{\mathrm{H{\sc I}}}~{}$than the sky as a whole. The inset in Figure 8 shows the USM of a small region of sky. The structures highlighted tend to be elongated, a fact discussed further below.

The FIR/$\mathrm{N}_{\mathrm{H{\sc I}}}~{}$associated with the small-scale structures is consistently higher than the sky-averaged value. This suggests that high spatial frequency Hi intensity structures are physically distinct from the surrounding medium, beyond simply having higher density. There are three physical effects that are expected to raise the FIR/$\mathrm{N}_{\mathrm{H{\sc I}}}~{}$ratio. These are (e.g. Ysard et al., 2015; Nguyen et al., 2018):

Optically thick Hi emission, which will lower the Hi column without affecting the associated FIR emission. 2. 2.

Spatially correlated molecular hydrogen (H2), which depletes the population of atomic hydrogen, thereby lowering $\mathrm{N}_{\mathrm{H{\sc I}}}$. 3. 3.

Increased dust emissivity or dust-to-gas ratio associated with dense gas, which will raise the FIR emission relative to $\mathrm{N}_{\mathrm{H{\sc I}}}$.

These three possible contributing effects are notoriously difficult to disentangle (Burstein & Heiles, 1978; Boulanger et al., 1996; Lenz et al., 2017). Nevertheless, all three of these physical situations are associated with the cold phase of the ISM (referring both to CNM and to the diffuse molecular phase).

The enhanced FIR/$\mathrm{N}_{\mathrm{H{\sc I}}}~{}$cannot be due to noise bias from weighting the data by the USM of an Hi channel map. The primary effect of radiometer noise in a narrow velocity channel will be to add spurious power both in that narrow channel map and in the $\mathrm{N}_{\mathrm{H{\sc I}}}~{}$map. This correlated noise will correspond to an enhancement in $\mathrm{N}_{\mathrm{H{\sc I}}}~{}$but not in the dust tracer. This bias will thus tend to suppress, rather than enhance, the FIR/$\mathrm{N}_{\mathrm{H{\sc I}}}~{}$ratio, and so noise bias cannot account for our result.

We therefore expect that the observed enhancement of the FIR/$\mathrm{N}_{\mathrm{H{\sc I}}}~{}$ratio correlated with small-scale Hi intensity structure is because these channel map structures are preferentially CNM. Although this is, to our knowledge, the first measurement showing that high spatial frequency Hi channel map structures are associated with higher FIR/$\mathrm{N}_{\mathrm{H{\sc I}}}$, the picture that CNM structures are preferentially small-scale is consistent with the current understanding of the ISM phase distribution. Indeed, while CNM accounts for roughly half of the Hi mass fraction, the CNM volume filling fraction is of order 1% of the WNM volume (Heiles & Troland, 2003; Murray et al., 2018b). This alone implies that the CNM is preferentially distributed into small-scale structures: the CNM is clumpy compared to the extended, diffuse WNM.

4 Discussion

We have shown that the small-scale structure in narrow Hi channel maps is predominantly cold phase material, and the preferential elongation of these density structures in the direction of the ambient magnetic field is responsible for the measured alignment between Hi intensity structure and the plane-of-sky magnetic field orientation. Other observations bolster the picture that the magnetically aligned Hi intensity structures are CNM. Filamentary structures in narrow Hi channel maps have FWHM linewidths consistent with CNM temperatures (Clark et al., 2014; Kalberla et al., 2016). Most strikingly, observations of Hi self-absorption in the Riegel-Crutcher cloud show that the CNM is organized into slender, linear structures that are very well aligned with the local magnetic field (McClure-Griffiths et al., 2006). Our results suggest that, if Hi self-absorption imaging were possible for more of the sky, such a cloud structure would not be uncommon. Indeed, taken together with the results of Clark et al. (2014, 2015), we suggest that the CNM is generically organized not into spherical clumps, but into magnetically aligned filamentary structures. The prediction of this work, that high spatial frequency Hi intensity structures are preferentially associated with cold gas, should be tested with absorption line measurements of the Hi spin temperature. The Square Kilometer Array (SKA) and the Australian SKA Pathfinder (ASKAP) will allow this picture to be tested in many different Galactic environments (McClure-Griffiths et al., 2015; Dickey et al., 2013).

How do the magnetically aligned CNM features in Hi channel maps compare to magnetically aligned FIR intensity structures? Planck Collaboration Int. XXXV et al. (2016) used the spatial gradient of the column density to characterize the local orientation of dust structures. The authors reported an alignment between the dust intensity structures and the magnetic field for low column densities, and quantified the correlation statistically using the Histogram of Relative Orientations (Soler et al., 2013). Planck Collaboration et al. (2016c) similarly find that dust structures are well aligned with the projected magnetic field at lower column densities. Lazarian & Yuen (2018) refer to the spatial gradient of column density maps as “intensity gradients”, and argue that they do not trace the magnetic field as well as gradients of velocity channels. The authors find that gradients of thin channel maps (“velocity channel gradients”) are better correlated with the local magnetic field orientation than gradients of either the column density map (“intensity gradients”) or of the first moment of the velocity (“velocity centroid gradients”) (Yuen & Lazarian, 2017; González-Casanova & Lazarian, 2017). Lazarian & Yuen (2018) (and references within) argue that this is because thin channel maps trace the velocity-magnetic field alignment, and so the higher-fidelity alignment of narrow velocity channels confirms the theoretical expectation that the density-magnetic field correlation is weak compared to the velocity-magnetic field correlation (e.g. Cho & Lazarian, 2003; Passot & Vázquez-Semadeni, 2003). The present work overturns that interpretation.

Why, then, are the Hi channel map structures better aligned with the magnetic field than structures in $\mathrm{N}_{\mathrm{H{\sc I}}}$? In their analysis of the relative orientation of the magnetic field and ridges in the FIR emission, Planck Collaboration et al. (2016c) estimate volume densities for the dust filaments that are consistent with CNM or diffuse molecular gas, although they primarily analyze higher column density features (median $\mathrm{N}_{\mathrm{H}}=1.2\times 10^{21}$ cm*-2*) than considered here. Planck Collaboration et al. (2016c) find that dust structures are better aligned with the projected magnetic field at lower column densities. To the extent that the magnetic alignment is stronger for structures in the Hi channel maps, it may be because diffuse CNM structures have a higher contrast relative to the surrounding medium when the Hi intensity is integrated over a narrow velocity range. Projection effects and line-of-sight effects can cause the observed correlation between the plane-of-sky magnetic field and intensity structures to vary from the true relative orientation between the three-dimensional magnetic field and density structures. 21-cm observations may thus be a useful tracer of the three-dimensional magnetic field, because line-of-sight effects can be probed by analyzing multiple channel maps along the line of sight (Clark, 2018).

This work calls into question the interpretation of other analyses that are premised on a significant contribution of the turbulent velocity field to structures in narrow Hi velocity channels. Based on this expectation, Lazarian & Pogosyan (2000) proposed that the three-dimensional velocity spectral index can be extracted from spectral observations by measuring the power spectral index of channel maps of varying $\delta v$ (see also Pogosyan & Lazarian, 2005). Lazarian & Pogosyan (2000) predict a steepening of the power spectrum slope with increasing $\delta v$, as a result of a decreasing contribution of the velocity field to the intensity fluctuations. One observational test of this picture in Hi is given in Pingel et al. (2013), where indeed the slope of the two-dimensional spatial power spectrum (SPS) steepens with increasing velocity slice thickness. Our work suggests that an alternative explanation for this result should be considered: that the intensity contribution of CNM structures increases with decreasing $\delta v$. If the CNM has a shallower power spectrum than the warmer-phase gas, this would also reproduce the qualitative trend of steeper SPS for larger $\delta v$. The preferential distribution of CNM into small-scale structures is consistent with a shallower SPS for the CNM. Other observational analyses that interpret shallower channel map spectral indices as a signature of velocity fluctuations should be reevaluated in light of this work (e.g. Dickey et al., 2001), and should include a careful treatment of the noise (e.g. Blagrave et al., 2017).

5 Conclusions

We draw a number of conclusions about the physical nature of Hi intensity maps, the structure of the multiphase ISM, and the interpretation of synthetic channel maps. These are summarized as follows:

• •

Small-scale structures seen in narrow velocity channels are dust-bearing density structures. A description of Hi intensity structures as being entirely or predominantly velocity caustics is incorrect.

• •

Techniques based on the gradient of intensity structures in narrow spectral channels successfully trace the orientation of the interstellar magnetic field because they measure the orientation of density structures that are preferentially aligned with the magnetic field, not because they measure an imprint of the turbulent velocity field.

• •

Small-scale Hi channel map structures have higher FIR/$\mathrm{N}_{\mathrm{H{\sc I}}}~{}$than their surroundings, which implies that they are associated with colder, denser phases of the ISM.

• •

The data presented are consistent with a picture of the neutral ISM in which the colder-phase gas is measurably distributed into small-scale, anisotropic structures, preferentially aligned with the magnetic field. The WNM, by contrast, is more extended and isotropic.

• •

A number of numerical studies predict a significant influence of the turbulent velocity field on the structure of Hi channel maps, in conflict with the data. We discuss the origin of some misleading theoretical expectations. Neglecting thermal broadening is unphysical and artificially inflates the intensity contribution from velocity caustics. Thermal broadening in subsonic Hi smooths out velocity-induced intensity fluctuations and correlates channel maps with the density field. Supersonic motions dynamically correlate the three-dimensional density and velocity fields. We demonstrate that these effects correlate synthetic Hi channel maps near the peak of the line with the synthetic column density map regardless of sonic Mach number.

• •

A shallower spatial power spectrum slope measured in narrower velocity channels is qualitatively consistent with a higher contribution to the intensity in narrow channels from CNM gas. Previous work has attributed such a change in the spatial power spectrum to velocity fluctuations, and used this to derive the slope of the turbulent velocity power spectrum in the gas. This work calls those analyses into question.

This work calls for a significant reassessment of many observational and theoretical studies of turbulence in Hi.

We thank François Boulanger for many enlightening discussions over the course of this investigation, and for reading an early version of this manuscript. We thank Claire Murray and Ludovic Montier for illuminating conversations on this work and related topics. We thank Cédric Colling and Patrick Hennebelle for the numerical simulations. This work took place in part under the program Milky-Way-Gaia of the PSI2 project funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02. We thank the other participants of that program for stimulating discussion. We thank the Flatiron Institute and Chang-Goo Kim for hospitality at the Center for Computational Astrophysics, where some of this analysis took place. We thank the anonymous referee for a thoughtful and prompt review. S.E.C. is supported by NASA through Hubble Fellowship grant #HST-HF2-51389.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. This publication utilizes data from Galactic ALFA Hi (GALFA-Hi) survey data set obtained with the Arecibo L-band Feed Array (ALFA) on the Arecibo 305m telescope. The Arecibo Observatory is operated by SRI International under a cooperative agreement with the National Science Foundation (AST-1100968), and in alliance with Ana G. Méndez-Universidad Metropolitana, and the Universities Space Research Association. The GALFA-Hi surveys have been funded by the NSF through grants to Columbia University, the University of Wisconsin, and the University of California.