© ESO 2021

1 Introduction

The interstellar medium (ISM) is shaped by many processes. Energy is injected by sources such as supernovae and stellar winds. This energy affects the dynamics of the ISM: it stimulates turbulence that cascades down from large to small scales. These processes have an imprint on the dynamics and the distribution of diffuse neutral hydrogen (H I) as well as on the physical state of the H I, which is characterized by the composition in different phases. Gas traced by H I emission is not completely neutral; as a consequence, the neutral phase is coupled to ions and the dynamics of the H I is linked to the magnetic field. This process is usually described as flux freezing (Heiles & Crutcher 2005). Gas and dust are well mixed (Clark et al. 2019; Kalberla et al. 2020). One of the most important tracers of the magnetic field is polarization from aligned dust grains, which emit the absorbed starlight in the far infrared (FIR; Heiles & Crutcher 2005). A close relation between H I and FIR emission and the magnetic field is expected (Clark & Hensley 2019).

Recent investigations have shown that the H I at high Galactic latitudes is organized in thin filamentary structures, sometimes described as fibers (Clark et al. 2014). These filaments are associated with FIR filaments (Kalberla et al. 2016). The morphology of these filaments suggests that the structures have been shaped by the magnetic field, and it has been shown that they are indeed extremely well aligned with the plane-of-sky magnetic field as probed by both starlight polarization (Clark et al. 2014) and polarized dust emission (Clark et al. 2014; Kalberla et al. 2016). Thus, H I and FIR structures are coupled to the magnetic field, and the morphology of these structures may be used to study magnetism in this environment. The most recent publications in this field have shown that H I filaments are coherent structures, indicating ordered polarized emission analogous to FIR polarization (Clark 2018; Clark & Hensley 2019). In turn, these coherent FIR and H I structures can be used to constrain models for the three-dimensional orientation of the interstellar magnetic field and the twisting of the magnetic flux tubes.

Magnetohydrodynamics (MHD) simulations by Planck Collaboration Int. XX (2015) are based on a large-scale anisotropic component of the magnetic field, with an additional turbulent component due to velocity perturbations imposed on converging flows that drive turbulence. A different assumption is that there is a uniform magnetic field with isotropic turbulent fields in different three-dimensional layers (Planck Collaboration XII 2020). We observe a systematic bending of the filaments that can be used to probe models on turbulence. Magnetic fields exert a tension force that opposes bending; they behave like elastic strings threading the fluid. Accordingly, the field strength and the field-line curvature are affected. Field tangling in the case of a small-scale turbulent dynamo, also called a fluctuation dynamo (St-Onge & Kunz 2018), leads to distinct predictions for the curvature distribution (Schekochihin et al. 2002, 2004) that can be verified.

Previous investigations of the coupling between the magnetic field and FIR and H I emission by Clark & Hensley (2019) made use of the Rolling Hough Transform (RHT), a technique from machine vision for detecting and parameterizing linear structures (Clark et al. 2014). This method, as applied by Clark & Hensley (2019), measures the linearity of structures in a region of diameter 75′ around each pixel, and so the effective resolution is rather limited. These authors considered angular resolutions between 80′ and 160′. Here we study the relationsbetween H I and FIR structures at the highest possible resolution. We use HI4PI survey data (HI4PI Collaboration 2016) that have beam widths between 10.′8 and 14.′5, the best currently available resolutions for all-sky data. However, our Hessian analysis is pixel based and therefore by construction limited to a resolution of18′. This is still an improvement by about a factor of four to nine, and we present, for the first time, an all-sky high resolution study of H I filaments. In Sect. 2 we explain our data reduction, the Hessian analysis, and the comparison of filamentary structures in H I and FIR at 857, 545, and 353 GHz. We find coherence between H I and FIR and explore the physical conditions of these H I filaments in Sect. 3. Curvatures in the case of magnetized filaments imply a magnetic tension with back reactions on the filament shape. We determine the distribution of filament curvatures in Sect. 4 and relate our results to predictions for the small-scale turbulent dynamo. Changes in curvatures between the central parts of the filaments and their environment are discussed in Sect. 5. We summarize our results in Sect. 6.

2 Observations and data reduction

We use HI4PI H I observations (HI4PI Collaboration 2016), combining data from the Galactic all sky survey (GASS; Kalberla & Haud 2015), measured with the Parkes radio telescope and the Effelsberg-Bonn H I Survey (EBHIS; Winkel et al. 2016) with data from the 100 m telescope. The H I data from both surveys have been gridded to independent nside=° 1024 HEALPix databases, retaining the original resolution for each of the telescopes (14.′5 for the GASS and 10.′8 for the EBHIS). Subsequently these H I profiles were decomposed into Gaussian components (Haud 2000; Kalberla & Haud 2015, 2018).

To compare the H I data with the FIR emission we used the most recent maps from the Public Data Release 4 (PR4; Planck Collaboration Int. LVII 2020) at frequencies of 857, 545, and 353 GHz¹. These maps have been generated using the NPIPE data processing pipeline, a natural evolution of previous Planck analysis efforts. PR4 data need to be corrected for the Solar dipole. We used parameters from Planck Collaboration Int. LVII (2020) with longitude 263. °986 ± 0. °035 and latitude 48. °247 ± 0. °023 and amplitude of 3366.6 ± 2.7 μK. The maps are published on an nside = 2048 HEALPix grid. We downgraded the maps to an nside = 1024 grid using the ud_grade software from the HEALPix software distribution. As detailed later, we smoothed some of the maps using the HEALPix smoothing program.

Our aim is to characterize filamentary structures that are common to the H I and FIR emission. Several different methods have previously been used to detect filamentary structures. Clark et al. (2014) applied RHT to data from the Galactic Arecibo L-Band Feed Array H I survey (GALFA-H I Peek et al. 2018) and to GASS data. Clark & Hensley (2019) applied the same approach to HI4PI data. Kalberla et al. (2016) applied unsharp masking (USM) to suppress large-scale structures in HI4PI data on scales below 30′. In addition they used the Hessian matrix to work out the all-sky distribution of filamentary structures in analogy to similar investigations by Schisano et al. (2014) and Planck Collaboration Int. XXXVIII (2016) in the FIR and submillimeter range. Recently Soler et al. (2020) used the Hessian analysis to analyze the filamentary structure in the disk of the Milky Way. These authors have shown that results from alternative algorithms, FilFinder (Koch & Rosolowsky 2015) and RHT (Clark et al. 2014), are consistent with those from the Hessian matrix. In the following we describe the Hessian analysis and the data processing as applied by us.

2.1 Hessian analysis

We classified filamentary structures as enhancements in intensity that are organized as elongated regions with a preferred spatial orientation. Such features on scales of arc minutes up to more than 20 degrees are often approximated as one-dimensional structures and termed fibers (Clark et al. 2014), but they may also be bent, resembling worms (Heiles 1984; Kalberla et al. 2016). Filaments may also be caused by projection effects from two-dimensional structures, such as sheets or walls, seen edge on.

As a tool to classify structures as filament-like, we used the Hessian operator H, which is based on partial derivatives of the intensity distribution $$H(x,y)\equiv \left(\begin{array}{cc}{H}_{xx}& H_{xy}\\ H_{yx}& H_{yy}\end{array}\right).$$(1)

Here x and y refer to true angles in longitude x = l cosb and latitude y = b. The second-order partial derivatives are H_xx = ∂²I∕∂x², H_xy = ∂²I∕∂x∂y, H_yx = ∂²I∕∂y∂x, and H_yy = ∂²I∕∂y².

The eigenvalues of H, $${\lambda }_{\pm }=\frac{(H_{xx}+H_{yy})\pm \sqrt{{(H_{xx}-H_{yy})}^2+4H_{xy}{H}_{yx}}}{2},$$(2)

describe the local curvature of the features; λ₋ < 0 is in direction of least curvature and indicates filamentary structures or ridges. Alternatively λ₊ < 0 highlights blobs. Equation (2) with λ₋ < λ₊ implies that the blobs are not independent but located as enhancements along the filaments.

The Hessian operator H has first been used for the analysis of H I data by Kalberla et al. (2016) but for consistency with Clark et al. (2014), Clark & Hensley (2019), Jow et al. (2018) and Soler et al. (2020) we follow here the definition x = l cosb. The derivatives are calculated by us from symmetric differential quotients. Due to the interleaved structure of the HEALPix database it is necessary to interpolate the HEALPix data in Galactic longitudes and we use a cubic spline interpolation for this purpose. The Hessian matrix is symmetrical, H_xy = H_yx but to minimize possible numerical artifacts we calculate both ways independently. All calculations are done in double precision. The numerical procedure takes only a few seconds for an nside = 1024 data slice.

The local orientation of filamentary structures relative to the Galactic plane is given by the angle $$\theta =\frac{1}{2}\arctan \left[\frac{H_{xy}+H_{yx}}{H_{xx}-H_{yy}}\right],$$(3)

in analogy to the relation $${\theta }_S=\frac{1}{2}\arctan \frac{U}{Q},$$(4)

that can be derived from polarimetric observations that provide the Stokes parameters U and Q.

This analogy between these two relations is intriguing but before we use angular orientations θ from Eq. (3) we need to understand the meaning of this measure in comparison to θ_S. Results according to Eq. (4) need independent observations of U and Q, two of the four Stokes parameters. In the case of Eq. (3) we are limited to a single observation of the Stokes parameter I. The orientation angle θ is derived merely from the morphology of the intensity distribution. The advantage in this case is however the high signal-to-noise ratio (S/N) of I, while the polarization information entering Eq. (4) is typically much closer to the noise level.

The derivatives in Eq. (1) are calculated over a fixed grid of 5 × 5 pixels, centered on x, y. This implies that θ from Eq. (3) is a measure for the local intensity fluctuations within this grid only. The Stokes parameters U and Q from polarimetric measurements are integrated along the line of sight and depend on the telescope beam. The 5 × 5 pixel filter of the Hessian operator is only meaningful if it is matched to the spatial scale under consideration, in our case defined by the resolution of the data product. The sensitivity to structures on scales that are large compared to the 5 × 5 pixel filter is low. For a more general application a multi-scale morphology approach may be necessary (Aragón-Calvo et al. 2007).

We analyzed HEALPix databases with nside = 1024; therefore, the effective pixel resolution is 3.′44 (Górski et al. 2005) and, accordingly, the H(x, y) filter acts on a scale of 18′. The FIR data, in comparison to the H I channel maps, are sensitivity limited (see, e.g., Figs. 3 and 4 of Kalberla et al. 2016). To avoid unnecessary noise amplification when deriving first and second derivatives, we smoothed the Planck 857 GHz data to an effective resolution of 18′. Thus we matched the full width at half maximum (FWHM) of the smoothing filter to the 5 pixels of the grid used by the Hessian analysis. Such an adapted smoothing leads to significant reduction of the uncertainties (see Sect. 2.3). In Fig. 1 we display the resulting distribution of eigenvalues that characterize filaments and blobs for FIR structures at 857 GHz.

For the H I data we used Gaussian components from Kalberla & Haud (2019) with a FWHM beamwidth of 14.′5 for the GASS and 10.′8 for the EBHIS. In comparison to the FIR databases there is a slight mismatch in spatial resolution but we intend to analyze correlations between FIR emission and H I at the highest possible resolution. A smoothing would disable the use of the Gaussian database since the fitting is a nonlinear process. Our analysis is hierarchical in the sense that we always use the 857 GHz FIR data as a template to search for structures at other frequencies. A slightly smaller FWHM width for the H I data does not violate the Nyquist theorem since both databases are on the same HEALPix grid. Tests with several different smoothing kernels indicated anyway that the derived eigenvalues do not depend critically on the smoothing kernel used.

2.2 Comparison of filamentary structures in different databases

To compare different data sets we apply simultaneously a Hessian analysis to both data sets. According to our definition of filaments as intensity enhancements that are organized and have a preferred spatial orientation we demand that both databases have at the same position significant eigenvalues λ₋ < 0. If this is the case we calculate and compare the orientation angles θ.

Clark & Hensley (2019) and Jow et al. (2018) have previously considered similar comparisons and developed different strategies to measure the alignment of filamentary structures. For a deeper discussion we refer to the original publications; we use these measures but define them below only briefly.

According to Planck Collaboration Int. XXXII (2016) and Clark & Hensley (2019), the angular difference between the orientation angles θ₁ and θ₂ at each position of two data sets can be calculated by $$\delta \theta =\frac{1}{2}\text{arctan}\left[\frac{\text{sin}(2{\theta }_1)\text{cos}(2{\theta }_2)-\text{cos}(2{\theta }_1)\text{sin}(2{\theta }_2)}{\text{cos}(2{\theta }_1)\text{cos}(2{\theta }_2)+\text{sin}(2{\theta }_1)\text{sin}(2{\theta }_2)}\right],$$(5)

and a mean degree of alignment, $$\xi =\text{ left}<\text{cos}\varphi \text{ right}>,$$(6)

can be defined for $$\varphi =2\delta \theta .$$(7)

A different metric, called projected Rayleigh statistic (PRS), was advocated by Jow et al. (2018), $$\text{PRS}=\sqrt{\frac{2}{N}}{\sum }_i\text{cos}{\varphi }_i,$$(8)

and the uncertainty of this measure can be estimated from $${\sigma }_{\text{PRS}}^2=\frac{2{\sum }_i{\text{cos}}^2{\varphi }_i-{(\text{PRS})}^2}{N}.$$(9)

2.3 Proof of concept

To test our data processing and the performance of the different alignment measures we use two neighboring Planck frequency channels. Both probe essentially the same dust foreground spectral energy distribution at slightly different frequencies. Data collected close in frequency should have the highest degree of similarity and thus the best spatial correlation of two independently measured data sets. For a quantitative comparison we compare first FIR data at 857 and 545 GHz. These data have been processed with the NPIPE pipeline and are claimed to have currently the best internal consistency between the various frequency channels (Planck Collaboration Int. LVII 2020). However, despite noise and differences in instrumental performance, we expect deviations when comparing 857 and 545 GHz data. For frequencies ν > 217 GHz the contribution of the Galactic foreground is dominated by the ISM, at lower frequencies the Galactic synchrotron and free-free radiation contribute significantly. The cosmic microwave background (CMB) contribution is increasing for lower frequencies and the relative contribution from the cosmic infrared background (CIB) can be considered as a significantly fainter and spatially unrelated offset to the H I emission (Odegard et al. 2019, their Fig. 3). In general, for Planck 545 GHz and 857 GHz frequency channels the Galactic foreground dominates across the sky. The total emission for the different grain models is highest at 857 GHz (Draine & Fraisse 2009, Fig. 7). Despite the fact that across the sky some fluctuations in the spectral energy distributions are expected, a good agreement is observed at 545 and 857 GHz (Irfan & Bobin 2018, Fig. A1). These two data sets are best suitable to test how far the application of the Hessian operator H and the subsequent determination of the angular orientation of filamentary structures (Eq. (3)) lead to coherent results at both frequencies.

For a basic test of our data analysis we use first unsmoothed 857 and 545 GHz data, downgradedto a HEALPix grid with nside = 1024. We select structures with eigenvalues λ₋ < −1.5 K deg⁻² at 857 GHz. To determine coherence we apply the measures according Eqs. (5)–(9), the results are listed in Table 1. We repeat analysis and comparison for 857, 545, and 353 GHz data after Gaussian smoothing of all data sets to an effective resolution of 18′ FWHM.

The distribution of angular differences δθ (Eq. (5)) for all sky 857 and 545 GHz data is shown in Fig. 2. From this plot it is obvious that for the original unsmoothed data the angular differences δθ are limited by noise. However, applying to the data a Gaussian smoothing with a kernel that is matched to the 5 × 5 pixel Hessian operator H leads to significant improvements with a lower spread for the angular alignment deviations.

The probability density distribution for angular differences δθ is a circular distribution in the range − π to + π. A normal distribution is in this case approximated by the von Mises distribution (e.g., Jow et al. 2018) $$\mathcal{F}(\delta \theta |\text{upmu},\kappa )\equiv \frac{e^{\kappa \cos \delta \theta -\text{upmu}}}{2\pi I_0(\kappa )},$$(10)

with κ approaching for narrow distributions the reciprocal of the dispersion of a normal distribution, κ ~ 1∕σ². I₀ is the modified Bessel function of order 0. Fitting $\mathcal{F}$ we found in cases with narrow δθ distributions virtually identical distributions (within the thickness of the lines plotted) in comparison to a Gaussian fit. In the case of broad distributions the $\mathcal{F}$ fit was found to diverge while the Gaussian fit was still stable. For this reason we plot Gaussian approximations and cite Gaussian parameters only. The shapes of the δθ distribution in Figs. 2 and 3 are difficult to fit. We also considered a simple estimate by measuring the width of the δθ distributions at half maximum without fitting, but such a measure seriously underestimates the width in the case of a strong central peak (central panels of Fig. 3). Central peaks and extended wings are best taken into account by a Voigt function and we found that the Voigt function approximates noisy data best².

To compare the distribution of angular differences between different data sets in a quantitative way we list in Table 1 the surface filling factor f for common filamentary structures, the dispersions σ_Gauss and σ_Voigt and further alignment indicators ξ, PRS, and σ_PRS according to Eqs. (5)–(9).

Fig. 1 All-sky Mollweide displays of derived FIR structures at 857 GHz. Top: distribution of eigenvalues λ− that characterize filamentary structures. Bottom: distribution of λ+, characterizing local enhancements (blobs) along filaments.

Fig. 2 Histograms of angular alignment deviations for filamentary structures according to Eq. (5). We compare 857 and 545 GHz data, unsmoothed as observed and smoothed to a resolution of 18′ FWHM.

2.4 Alignment between FIR and H I column density structures

We aim to determine the degree of alignment between FIR and H I structures. In the following we use only smoothed 857 GHz data for the FIR since these have among the Planck observations the best S/N (Planck Collaboration Int. XVII 2014, Fig. C.1). In contrast to the H I distribution the 857 GHz data are noise limited (Kalberla et al. 2016) but, as shown in Fig. 2, smoothing helps to improve our analysis significantly. For the H I we chose first a velocity range of − 50 < v_LSR < 50 km s⁻¹. This range comprises about 75% of the total H I emission. It also accounts for the phenomenon of dusty intermediate-velocity clouds (Röhser et al. 2016) and thus is complete in the sense that it samples all H I emission that is associated with FIR emission. For a first comparison we calculate the integral across this radial velocity range (the moment zero map) and derive column densities N_HI. As mentioned before we do not smooth H I data.

Table 1 shows the resulting alignment measures for two cases, using all sky data and restricting the analysis to high Galactic latitudes, |b| > 20°. The high latitude comparison shows the best agreement, the distribution function for angular deviations is best fit with a Voigt profile with a dispersion of about 20°. We display this distribution in Fig. 3 (top right) for comparison with alignment deviations between FIR at 857 and 353 GHz (Fig. 3 top left). There is some alignment between FIR filaments and H I but apparently the N_HI filaments fail a detailed correlation with the FIR.

Fig. 3 Histograms of angular alignment deviations according to Eq. (5) for filamentary structures. Top left: Planck 857 GHz compared with Planck 353 GHz and angular orientations according Eq. (4), |b| > 20°. Top right: Planck 857 GHz compared with H I column densities for |vLSR| < 50 km s−1, |b| > 20°. Center left: Planck 857 GHz compared with best fit single channel H I filaments, all sky. Center right: Planck 857 GHz compared with best fit single channel H I filaments, |b| > 20°. Bottom left: Planck 857 GHz compared with best fit H I filaments at a channel width of 16 km s−1, all sky. Bottom right: Planck 857 GHz compared with best fit H I filaments at a channel width of 16 km s−1, |b| > 20°.

2.5 Alignment between FIR and H I in single channels

Cold H I filaments are observed to have average FWHM line widths of 3 km s⁻¹ (e.g., Clark et al. 2014 or Kalberla et al. 2016). This empirical result, typical for a cold neutral medium (CNM), may imply that coherence between FIR and H I can exist only for a restricted range in radial velocity. The filamentary CNM structures are embedded in a more diffuse and warmer environment, the warm neutral medium (WNM). Dust and gas are well mixed and the CNM is in any case found to be associated with WNM (e.g., Kalberla & Haud 2018). The question arises whether we can distinguish between different H I phases that can be related to the FIR filaments. Hence it is mandatory to consider the line-width of possible H I counterparts of FIR filaments. For a turbulent medium it is in addition a matter of debate whether filaments are real density structures or just velocity caustics, fluctuations imprinted by the turbulent velocity field (Clark et al. 2019). In the case of velocity caustics it is assumed that density and velocity fields are completely uncorrelated. Sufficient large velocity channel widths of ~ 17 km s⁻¹, typical for a WNM, need accordingly to be used to characterize the density field.

To verify the different hypotheses we modify our analysis and probe the alignment between FIR and H I structures depending on both, radial velocities and velocity channel widths. In five independent runs we generate separate HEALPix databases with 101, 51, 25, 13, and 7 velocity slices in the range − 50≲v_LSR≲50 km s⁻¹, each slice with a width of 1, 2, 4, 8, and 16 km s⁻¹, respectively. We apply the Hessian operator individually to each of these slices and calculate λ_± (l, b, v_LSR) and associated angles θ(l, b, v_LSR). Next we compare for all HEALPix positions the angular alignment between FIR at 857 GHz and H I at each of the velocity channels with different channel widths. The Hessian analysis is performed in all cases over a 5 × 5 pixel region centered on (l, b). For each of the five different setups independently the velocity with the minimum scatter in the δθ distribution (or alternatively the best alignment) defines the best fit filament velocity v_fil at this position. Selecting then λ_±(l, b, v_fil) we reduce the three-dimensional λ_±(l, b, v_LSR) distributions to two dimensions. For each of the five velocity channel widths we obtain this way an H I distribution that can be tested for internal coherence in velocity (discussed in Sect. 2.6) but also for consistency with FIR at 857 GHz.

The derived alignment measures are summarized in Table 1. During our investigations we preferred to inspect the distributions of the angular alignment deviations, characterized by σ_Voigt, but found none of the listed alignment measures to be superior over the others. However, all of these parameters show a common trend. To our understanding the collapsed filamentary H I distribution that was extracted for the velocity slice with a width of 1 km s⁻¹ results in the best fit. However, we also note opinions that the 8 and 16 km s⁻¹ integrated maps might be most appropriate for this kind of analysis and should be preferentially analyzed. In the following we consider both solutions, for 1 and 8 or 16 km s⁻¹ channel widths in some detail.

Figure 3 shows examples of the derived distribution functions for angular alignment deviations, all sky and restricted to high latitudes |b| > 20°. Approximating the alignment distributions with Voigt functions we obtain σ_Vogit < 3° at high latitudes and σ_Voigt < 2° all sky in the case of 1 km s⁻¹ broad channels (Fig. 3 center). This comes even close to the correlation of the two neighboring Planck frequency bands, inspected in Sect. 2.3. We conclude that there is a very well-defined coherence between FIR and H I filaments at the narrowest velocity interval. Figure 3 shows at the bottom the distribution of angular alignment deviations for a channel width of δv_LSR = 16 km s⁻¹ with significantly larger dispersions.

2.6 Velocity dispersion along H I filaments

In this subsection we consider the problem of whether the correlation between Planck and H I data might be manufactured by the data processing. The correlation could be an artifact of selecting only the best of many noisy channels. In other words, the question is whether one could generate a map with a high correlation coefficient by using just N realizations of random noise given the number N is large enough.

To investigate this question we consider the general case of a three-dimensional distribution of H I filaments that may be correlated with FIR structures that are observable only in two dimensions. To establish a correlation, we strictly demand that such H I filaments are not only coincident and aligned with the FIR in position but are also homogeneous in the third dimension, thus uniform without irregularities in velocity. A distribution is usually considered to be uniform if its parameters are not arbitrary but lie between certain narrow bounds. Irregularities are then measurable as deviations.

To determine irregularities in velocity we define a velocity dispersion along the H I filaments $$\mathcal{V}\left(r,\delta \right)=\sqrt{\frac{1}{N}{\sum }_{i=1}^N{\left[v_{\text{LSR}}(r+{\delta }_i)-v_{\text{LSR}}(r)\right]}^2}.$$(11)

The sum extends over all pixels along the filament with positions (r + δ_i) within an annulus centered on r and having inner and outer radii δ_inner and δ_outer, respectively. According to our definition of the Hessian operator with adapted Gaussian smoothing over five pixels we select an inner radius δ_inner =9′. For the outer radius we chose δ_outer = 1°, a value that appears to be appropriate as discussed below in Sect. 5.2. We also tested smaller δ_outer values and found no significant biases when changing this parameter. We count only pixels along filaments and determine thus the local velocity dispersion $\mathcal{V}$ only along the filaments.

Figure 4 shows the distributions of velocity dispersions $\mathcal{V}$ for three different settings with channel widths of 1, 8, and 16 km s⁻¹. For 1 km s⁻¹ the filament velocities have a typical scatter $\mathcal{V}~6$ km s⁻¹ all sky and $\mathcal{V}~5.5$ at high latitudes. Close to the Galactic plane there is evidence for some confusion, causing the bump in the distribution function at $\mathcal{V}~15$ km s⁻¹. $\mathcal{V}$ includes systematicaleffects from linear velocity gradients. The measured velocity dispersions contain also statistical uncertainties in the determination of v_fil caused by the limited velocity resolution of 1 km s⁻¹ (Sect. 2.5). $\mathcal{V}$ can therefore only partly be caused by fluctuations in the velocity field, the FIR and H I filaments are rather well defined and coherent in velocity space.

Numerical studies of the condensation of the WNM into CNM structures under the effect of turbulence and thermal instability by Saury et al. (2014) indicate that the velocity field of the CNM reflects the velocity dispersion of the WNM. These authors estimate a CNM cloud-to-cloud velocity dispersion of 5.9 km s⁻¹, in good agreement with Fig. 4 in the case of a 1 km s⁻¹ channel width. The velocity dispersion distributionsfor the 8 and 16 km s⁻¹ integrated maps deviate by construction significantly from these theoretical estimates.

It remains to be discussed whether the correlation between FIR and H I filaments at 1 km s⁻¹ channel width could be flawed because the S/N in individual H I channel maps is much lower than in the integrated map. The Hessian analysis of individual channel maps may increase the uncertainties even more and with 101 opportunities a channel with a concordant angle can nearly always be found.

To verify this assertion we consider how far noise degrades the H I data (Kalberla & Haud 2019, Sect. A2.3). The system noise T_sys contains several independent contributions. Most important is the thermal noise from the receiver system and the elevation-dependent ground radiation including spill-over. These components are variable, but we use here an average thermal contribution T_sys = 30 K. The line signal T_B (v_LSR) adds to the noise contributionT_Noise(v_LSR) that can be approximated as (Haud 2000) $$T_{\text{Noise}}(v_{\text{LSR}})={\sigma }_{\text{av}}[T_{\text{sys}}+T_{\text{B}}(v_{\text{LSR}})]/T_{\text{sys}},$$(12)

where σ_av is the average noise level in the baseline, determined at velocities without H I line emission. We use this approximation of the radiometer equation to increase artificially the noise for all observed H I line profiles step by stepby a factor of two. For each of these cases all steps of the data processing pipeline are then repeated. Figure 5 shows the results. The velocity dispersions increase significantly and $\mathcal{V}$ is already unacceptable large for a noise enhancement by a factor N_f = 2. Likewise we are unable to get lucky imaging with higher additional noise injection.

Figures 4 and 5, using the dispersions $\mathcal{V}$, are somewhat abstract and it may be difficult to realize the consequences. In Fig. 6 we demonstrate therefore that the derived coherent velocity field in the upper panel is degraded significantly when enhancing the telescope noise level by a factor of N_f = 4. In the lower panel of Fig. 6 it is just possible to recognize a few H I filaments as coherent structures. Considering a noise amplification by a factor N_f = 8 (not shown), the image is completely decorrelated. Likewise it is not possible to generate or improve coherent structures for larger values of N_f. This result reflects a well known physical experience; the additional noise that we consider here is incoherent (independent in position and velocity) and by adding incoherent events it is not possible to generate coherent structures with a high correlation coefficient. Generating a coherent result in N random realizations as a sum of incoherent events demands N to be virtually infinite. In other words: the H I structures that are derived from Hessians at 1 kms resolution are not random but highly correlated, not only in position but also in velocity. Observers are working hard to get significant data. The average noise term in Eq. (12) scales with the integration time τ for individual positions as σ_av ∝ τ^−1∕2. The velocity field presented at bottom of Fig. 6 simulates observations that are therefore unacceptable short in observing time by a factor of 16.

Observers common experience is also that it is necessary to resolve structures to be analyzed. In the case of H I lines, CNM structures with an approximate FWHM width of 3 km s⁻¹ are common and according to the Nyquist sampling theorem (Shannon & Weaver 1975) a spectral resolution of at least 1/2 of this value is necessary, independent of any other requirements concerning the noise level. The possibility that averaging over broader channels would increase the correlation by reducing noise violates the sampling theorem. Applying such an averaging causes a loss of information about the H I structures and does according to Fig. 4 not improve the velocity dispersion along the filaments. Also from Fig. 3 and Table 1 we find for none of the different alignment measures any evidence for an improved correlation with increasing channel width.

This subsection on the velocity dispersion summarizes one of the most important results from our investigations. Comparing local orientation angles θ in the FIR at 857 GHz along filaments with H I orientation angles for 101 channels it is no surprise to find some agreement at one of these velocity channels. For structures that are unrelated this should result in a more or less random velocity distribution. Along the filaments we find however a surprisingly close agreement of the orientation angles with a well-defined coherence in velocity space (see Fig. 6), a necessary condition for alignment between FIR and H I structures along magnetic field lines as conjectured by Clark (2018) and Clark & Hensley (2019).

Fig. 4 Distribution of H I velocity dispersions $\mathcal{V}$ along FIR filaments at channel widths of 1, 8, and 16 km s−1.

Fig. 5 H I velocity dispersions $\mathcal{V}$ along FIR filaments at a channel width of 1 km s−1 depending on simulatednoise enhancements Nf = 2, 4, and 8 of thetypical observed noise level.

2.7 FIR and H I alignment conditions

To clarify conditions that lead to a good alignment between FIR and H I filaments, we need to sort out several effects that influence the data analysis.

First of all we need to consider data processing and observational uncertainties. The application of the Hessian operator, thus the calculation of first- and second-order derivatives, leads unambiguously to an amplification of uncertainties. In addition the definition of the 5 × 5 Hessianmatrix, applied to a HEALPix data cube introduces a spatial filtering, depending on the dimension defined by the nside parameter.

To minimize the uncertainties, also to minimize the calculation efforts, it is reasonable to use Gaussian derivatives for the calculation of the Hessian matrix (e.g., Soler et al. 2020). We do not follow this track but apply smoothing and calculation of the Hessian matrix independently, allowing H(x, y) to be determined for smoothed FIR and unsmoothed H I data with the same software. We adapt a Gaussian smoothing to the FIR data only, but with a FWHM smoothing kernel that matches the 5 × 5 pixel size of the Hessian operator. Assuming that the 857 and 545 GHz Planck FIR data are essentially a representation of similar ISM structures we reproduce the angular orientation of filamentary structures in 857 and 545 GHz FIR data according to Eq. (3) within σ_Voigt ~ 4. °3 all sky and σ_Voigt ~7. °5 toward high latitudes. The scatter in δθ decreases byabout a factor of four if we apply a Gaussian smoothing. Applying the same alignment test to smoothed 857 and 353 GHz FIR data, we find that thedispersions σ_Voigt between both FIR data sets increase roughly by a factor of 3.5; we explain this by sensitivity limitations, caused by the significant total intensity difference of the CMB and CIB in both frequency bands (Odegard et al. 2019, Fig. 3). Using angular orientations from Stokes parameters at 353 GHz in comparison to those from 857 GHz we find also an increased mismatch. As mentioned above, this mismatch is only partly due to sensitivity limitations. The application of the Hessian matrix causes a reduced sensitivity to low spatial frequencies but in Sect. 4.1 it will become clear that spatial filtering does not affect our analysis severely.

Considering similarities between 857 GHz FIR and H I we find some alignment for column densities but the correlation is weak. Searching for an angular alignment between FIR and H I in single channels we obtain a far better agreement but only for a velocity resolution of δv_LSR = 1 km s⁻¹ the agreement is perfect. The alignment measures are close to that of our test case at 857 and 545 GHz. Such an alignment can only be meaningful if the coherence between FIR and H I can be verified also in velocity space. Continuity in velocity along the filaments is discussed in Sect. 2.6 and later in Sect. 3.6. For the rest of this paper we consider only the best fit H I filaments with a velocity resolution of δv_LSR = 1 km s⁻¹.

To demonstrate the small-scale structure of orientation angles θ in filaments we show in Fig. 7 in the top and middle panels the distributions of θ for FIR at 857 GHz and for the H I in narrow velocity intervals. Both distributions are almost indistinguishable. At the bottom we display θ_S from Eq. (4), calculated from 353 GHz U and Q Stokes maps. Still, there is some correlation of this with maps in the upper panels (see also Table 1) but it is obvious that parameters derived from the Stokes maps differ from those derived by a Hessian analysis. Large-scale structures are suppressed by the Hessian analysis despite the fact that all maps in Fig. 7 have a similar spatial resolution. Some of the narrow structures resemble chromospheric fibrils that tend to be aligned with the magnetic field (Asensio Ramos et al. 2017). For the interpretation we need to be careful with respect to spatial filtering caused by the Hessian operator (Aragón-Calvo et al. 2007). Only at the very end of our paper in Sect. 4.1 it will become clear that neither sensitivity limitations nor biases caused by spatial filtering are affecting our conclusions severely.

The focus of our discussions in this section went from data processing to more physical influences like velocity widths and scale sizes. In the following we intend to explore astrophysical issues.

Fig. 6 Velocity field in the case of H I filaments derived for data with a 1 km s−1 channel width at an arbitrarily selected central position, l = 160°, b = 30° in gnomonic projection. The field size is 27. °7; only velocities in the range − 30 < vLSR < 30 km s−1 are displayed. Top: velocities from original telescope data. Bottom: results derived after simulating a noise amplification by a factor of Nf = 4.

Fig. 7 Orientation angles according to Eq. (3). Top: from the FIR distribution at 857 GHz.Middle: from the best fit H I distribution. Bottom: distribution of orientation angles according to Eq. (4) calculated from Stokes parameters U and Q of the 353 GHz FIR distribution.

Fig. 8 Distribution of harmonic mean H I Doppler temperatures in coherent FIR and H I filaments.

3 Coherence conditions

The close angular alignment between FIR emission and H I filaments at narrow velocity intervals implies that the FIR emission of such filaments has to come from regions where gas and dust are well mixed and in dynamical equilibrium. The best fit local H I velocities v_fil (l, b) are accordingly characteristic for the filamentary FIR emission structures. In this section we explore conditions and limitations for FIR and H I coherence in the ISM.

3.1 H I Doppler temperatures

H I filaments are cold with median Doppler temperatures of 220 K (e.g., Clark et al. 2014; Kalberla et al. 2016). For single dish observationswe cannot determine excitation temperatures, we can only use a Gaussian decomposition to derive Doppler temperatures T_D as upper limitsof the true exitation temperatures of individual H I components. In the case of several such components, we need to use harmonic mean Doppler temperatures defined as the brightness temperature (T_B) weighted average Doppler temperature for all Gaussian components along the line of sight, $$\langle T_{\text{D}}(v_{\text{LSR}})\rangle =\frac{{\displaystyle \int T_{\text{B}}}(v_{\text{LSR}})T_{\text{D}}\text{d}{v}_{\text{LSR}}}{{\displaystyle \int T_{\text{B}}}(v_{\text{LSR}})\text{d}{v}_{\text{LSR}}}.$$(13)

This ⟨T_D(v_LSR)⟩ definition is analogous to the definition of optical depth weighted average spin temperatures ⟨T_S (v_LSR)⟩ (e.g., Murray et al. 2020). Figure 8 shows the log-normal distribution of ⟨T_D (v_fil)⟩ for H I filaments at high latitudes |b| > 20° and all sky. While the all sky distribution may suffer from confusion caused by the Galactic plane, the high latitude data indicate that most of the gas in the filaments belongs to the CNM with low Doppler temperatures.

3.1.1 H I phase composition in filaments

To determine the phase composition of the coherent FIR and H I filaments we search at each position along the filaments for the Gaussian component with the center velocity v_Gauss with the least deviation |v_fil − v_Gauss| from the filament velocity. We determine the Doppler temperature T_D = 21.86 Δv² of this component, here Δv is the FWHM velocity width of this Gaussian.

Figure 9 shows the derived log-normal distribution of H I Doppler temperatures in filaments. We also apply a three-component fit to the T_D distribution.At high Galactic latitudes we find the CNM as the dominant component with a geometric mean T_D,CNM = 196 K. We get for the unstable lukewarm neutral medium (LNM) T_D,LNM = 567 K and T_D,WNM = 2156 K for the WNM. These values can be compared to characteristic values at high Galactic latitudes T_D,CNM = 283 K, T_D,LNM = 2014 K, and T_D,WNM = 11879 K, as determined by Kalberla & Haud (2018) from a Gaussian analysis, unconstrained to filamentary structures. The formal errors of these values are low, around 1%. We conclude that H I in coherent FIR and H I filaments is for all phases significantly colder than the H I outside such filaments.

Fig. 9 Distribution of H I Doppler temperatures of coherent FIR and H I filaments compared with a Gaussian three-component fit. Top: high latitude data, |b| > 20°. Bottom: all sky.

3.1.2 H I phase composition in blobs

The Hessian operator allows two eigenvalues, representing filaments and local enhancements (blobs) along the filaments, to be distinguished with Eq. (2). Figure 1 shows that there are numerous local enhancements at a resolution of 18′. Such structures are well defined at the resolution of 14.′5 for the GASS and 10.′8 for the EBHIS. We repeat the analysis from Sect. 3.1.1 and determine the phase distribution for these H I blobs.

At high Galactic latitudes we find with Fig. 10 a T_D distribution that is in its shape similar to Fig. 9 for the filaments. The CNM with a geometric mean T_D,CNM = 183 K is cold but the LNM with T_D,LNM = 282 K dominates the low temperature wing. This part is colder than the LNM of the filaments. Opposite for the WNM, here we get with T_D,WNM = 4766 K somewhat higher temperatures in comparison to the filaments. All sky we obtain T_D,CNM = 207 K, T_D,LNM = 360 K and T_D,WNM = 5528 K. In summary, local enhancements along filaments contain also predominantly cold H I gas. The formal parameters for the phase decompositions in Figs. 9 and 10 may not be very well defined but we clearly find in all phases significantly lower Doppler Temperatures than the averages determined previously by Kalberla & Haud (2018). Comparing these blobs with sources from the Planck catalog of Galactic cold clumps (Planck Collaboration Int. XLVIII 2016) we find that only 2397 out of 13242 sources listed in this catalog belong to the population of blobs. We conclude that most of the blobs shown in Fig. 1 share the properties of the filamentary structures; they are just condensations along the filaments.

Fig. 10 Distribution of H I Doppler temperatures of local enhancements (blobs) along filaments and their approximations by Gaussian distributions (green). To distinguish the distributions from H I phases, we separate the CNM (blue) from the warmer phase, LNM and WNM (red).

3.2 Velocity distribution

Figure 11 shows the velocity distributions for H I in coherent FIR and H I filaments. Filaments are local phenomena, they do not share high rotational velocities from the rotation curve in the Galactic plane. Using unconstrained all sky data we find that filaments are mostly in the range |v_LSR|≲30 km s⁻¹. At high Galactic latitudes we find barely filaments with |v_LSR| > 20 km s⁻¹. Our analysiscovers the range |v_LSR| < 50 km s⁻¹ but this limited range is complete and allows an unbiased description of filamentary structures. A Gaussian fit to the velocity distribution leads to a FWHM width of 16.6 km s⁻¹ at high latitudes and to 19.1 km s⁻¹ for the all sky case and confirms the narrow velocity distribution of small-scale filamentary structures derived previously from USM by Kalberla et al. (2016).

To check whether filament velocities may depend on the phase composition we distinguish velocity distributions according to contributions from the CNM and the warmer phases. Figure 11 shows that the warmer contributions tend to be asymmetric,indicating dynamical interactions.

Fig. 11 Histograms of the velocity distributions for filamentary structures and their approximations by Gaussian distributions (green). To distinguish the distributions from H I phases, we separate the CNM (blue) from the warmer phase, LNMand WNM (red).

3.3 Polarization: H I coherence in velocity

Clark (2018) investigated the alignment of H I features with the plane-of-sky magnetic field orientation. To quantify the correlation between H I structures and the magnetic field as characterized by the polarization fraction of 353 GHz dust emission, she proposed characterizing H I structures by their H I coherence or H I polarization. This metric is defined for the degree of coherence of the H I orientation as a function of velocity by Stokes-like parameters U_HI and Q_HI. Clark et al. (2014) and Clark & Hensley (2019) use the orientation of H I structures with angles derived from a RHT. Here we use angles θ according to Eq. (3), derived after application of the Hessian operator. These are weighted by the local H I brightness temperature T_B and integrated over the line of sight velocity v_LSR, $$U_{\text{HI}}={\displaystyle \int T_{\text{B}}}(v_{\text{LSR}})\text{cos}(2\theta (v_{\text{LSR}}))\text{d}{v}_{\text{LSR}}$$(14)

and $$Q_{\text{HI}}={\displaystyle \int T_{\text{B}}}(v_{\text{LSR}})\text{sin}(2\theta (v_{\text{LSR}}))\text{d}{v}_{\text{LSR}},$$(15)

and the H I coherence (Clark et al. 2014) or H I polarization (Clark & Hensley 2019) is defined as $$p_{\text{HI}}=\frac{\sqrt{U_{\text{HI}}^2+Q_{\text{HI}}^2}}{{\displaystyle \int T_{\text{B}}}(v_{\text{LSR}})\text{d}{v}_{\text{LSR}}}.$$(16)

Clark & Hensley (2019) applied this concept to map magnetically coherent regions of space. They have proven the existence of H I regions that are highly correlated with the 353 GHz Q and U maps of polarized dust emission observed by Planck on scales of 80′. Comparing 857 GHz data on scales of 18′ with H I emission from GASS and EBHIS we use this concept for our high resolution analysis and display the distribution of derived H I polarization data p_HI in Fig. 12. These distributions are close to Gaussians with peaks at p_HI ~ 0.3 but with extended wings up to p_HI ~ 0.8. As discussed below in Sect. 3.5 in more detail, the dust polarization fraction is limited and p_HI is not representative for the dust.

The derived distributions for the H I polarization p_HI in Fig. 12 are remarkably well defined and the peaks appear only little affected by confusion effects in the Galactic plane. On scales of 18′ p_HI is significantly enhanced in comparison to previous low resolution investigations. Planck Collaboration Int. XXXVIII (2016) determine from 353 GHz data on multipole scales 30≲l≲300 a mean polarization fraction of magnetized filaments of 11%. Clark (2018) and Clark & Hensley (2019) report H I polarizations around 15% on scales of 90′ and 80′, respectively, with a maximum near 30%. From 353 GHz Planck data at 1° resolution Planck Collaboration Int. XIX (2015) determine a maximum polarization of 19.8% and Planck Collaboration XII (2020) derive a maximum of 22% at 80′ resolution.

For all these publications the spatial filters were chosen in such a way that they highlight all the bright filaments. Our aim is to compare FIR and H I filaments at the highest possible resolution with the best sensitivity. Our data are matched to the 18′ resolution of the Hessian operator and can only be consistent with previous results if there are significant small-scale effects along the filaments and within the analyzed beam. The CNM as the coldest and densest part of the H I distribution dominates the structure of the filaments. The densest portions are in the very centers of the filaments and the magnetic field strengths correlates with volume density. These cold and dense parts cause the local bending for both, FIR and H I filaments with fluctuations in θ. Such a tangling must cause a de-correlation of the polarization observed at larger scales.

Fig. 12 Distribution of the H I polarization pHI according to Eq. (16), all sky and for high latitudes at |b| > 20°.

Fig. 13 Distribution of the polarization angle dispersion $\mathcal{S}$ along FIR filaments at a spatial resolution of 18′. The peak at $\mathcal{S}={52}^{°}$, annotated with a vertical line, indicates that many positions have significant deviations of polarization angles along the filaments. The black and red lines mark filaments with eigenvalues λ− < −1.5 K deg−2, and blue andorange lines represent filaments limited to λ− < −3 K deg−2.

3.4 Polarization angle dispersion

To derive the mean scatter of polarization angles along the FIR filaments at 857 GHz we use the polarization angle dispersion function $\mathcal{S}$, introduced in Planck Collaboration Int. XIX (2015) $$\mathcal{S}\left(r,\delta \right)=\sqrt{\frac{1}{N}{\sum }_{i=1}^N{\left[\theta (r+{\delta }_i)-\theta (r)\right]}^2}.$$(17)

The sum extends over all pixels along the filament with positions (r + δ_i) within an annulus centered on r and having inner and outer radii δ∕2 and 3δ∕2, respectively.According to our definition of the Hessian operator with adopted Gaussian smoothing over five pixels we select a lag of δ =18′. Contrary to the previous usage of $\mathcal{S}$ we count only pixels along filaments with eigenvalues λ₋ < −1.5 K deg⁻².

Figure 13 shows that polarization angles along filaments have considerable dispersions. We consider here for the moment only the black and red lines in that figure. The peaks of these distribution functions are close to 52°. Planck Collaboration Int. XIX (2015) and Planck Collaboration XII (2020) use $\mathcal{S}$ to quantify the regularity of the magnetic field. They conclude that a distribution peaking around $\mathcal{S}=\pi /\sqrt{12}(~{52}^{°})$ is characteristic for a chaotic and spatially completely uncorrelated distribution of polarization angles θ(r + δ_i) − θ(r). Opposite to this conclusion we demonstrate in Sect. 4 that the H I polarization angles follow a distinct nonrandom pattern if one takes the relation between orientation angles θ at neighboring positions into account. The observed orientation angels for neighboring positions are correlated due to systematic field curvatures with increased tangling on small scales. The definition for $\mathcal{S}$ (Eq. (17)) disregards such effects. As shown in Sect. 4.1, the distribution of curvatures is highly nonrandom, indicating that the filaments follow a systematical bending of the magnetic field lines. The peak at $\mathcal{S}~{52}^{°}$ is consistent with an uniform distribution of orientation angles without preferred orientation (e.g., Naghizadeh-Khouei & Clarke 1993, Appendix A) but may not be mistaken as a proof for a purely random distribution of filament structures. Here we consider only the question whether the extended wings for $\mathcal{S}\text{ ga}{90}^{°}$ could be caused by blending of unrelated filaments. To check for effects from such a blending we excluded in Eq. (17) positions with velocity deviations exceeding 10 km s⁻¹. We find no significant change in the histograms compared to the unbiased distribution shown in Fig. 13. In presence of the narrow velocity distribution derived in Sect. 3.2 blending may be unavoidable. At high Galactic latitudes on average 2.5–3 H I clouds are expected along the line of sight (Panopoulou & Lenz 2020). In Planck Collaboration XII (2020) it is assumed that the polarization angle dispersion measure is characteristic for a turbulence in layered structure of the magnetic field along the line of sight.

The unexpected strong peak of the polarization angle dispersion at $\mathcal{S}~{52}^{°}$ in Fig. 13 is an intrinsic property of the small-scale structures we observe along the filaments as a response to the Hessian operator. To check the dependence of $\mathcal{S}$ on the chosen confidence level λ₋ < −1.5 K deg⁻² that is used throughout this paper, we exacerbate this constraint by a factor of two. In response only the most prominent filamentary structures are considered. Toward high Galactic latitudes this constraint affects 50% of all filaments and all-sky 30% are discarded. However, the peak of the polarization angle dispersion distribution in Fig. 13 remains at $\mathcal{S}~{52}^{°}$ (blue and orange lines). On average about 109 positions along the filaments contribute to the measurement of $\mathcal{S},$ while in the case of λ₋ < −3 K deg⁻² only 103 positions are used.

3.5 Polarization angle dispersion versus polarization fraction

Planck Collaboration Int. XIX (2015) report that at 353 GHz $\mathcal{S}$ is in general anticorrelated to the observed polarization fraction. This result was confirmed by Planck Collaboration XII (2020) and Clark & Hensley (2019) report a similar relation also for H I polarization with a resolution of 160′. In the following we focus on the polarization angle dispersion along individual filaments at considerable smaller scales of 18′.

The $\mathcal{S}$ from Eq. (17) essentially measures fluctuations of the polarization angle perpendicular to the line of sight on scales of 3δ, three times the selected resolution and in our case on a scale of ~54′. Tangling causes fluctuations in the plane-of-sky magnetic field orientation. Figure 14 shows that, in this case, $\mathcal{S}$ does not depend on the H I polarization p_HI according to Eq. (16). The H I polarization p_HI is defined as a measure of coherence along the line of sight, hence across the H I fibers. The dispersion $\mathcal{S}$ is along the filament and both are apparently not correlated. Our results clearly contradict Clark & Hensley (2019, Fig. 11), who observe an anticorrelation between $\mathcal{S}$ and p_HI similar to the anticorrelation in the case of polarization fractions p₃₅₃ from Stokes parameters. H I polarization fractions and polarization angle dispersions on scales of 160′ and 18′ are not compatible. The scale differs by a factor of ~9 and we conclude that just on angular scales of ~18′ we are able to resolve details of the local magnetic field structure manifested in the bending or tangling of fibers perpendicular to the line of sight.

In Fig. 15 we display dependences between polarization angle dispersion along the same filaments as in Fig. 14, but now for polarization fractions determined from Stokes parameters at 353 GHz. These two-dimensional distributions confirm the previously published anticorrelations between polarization angle dispersion and polarization of FIR filaments. We also used the molecular cloud sample as considered by Planck Collaboration Int. XX (2015) for an additional comparison between polarization angle dispersions and polarization fractions. The result is shown in Fig. 16. The only difference from the previously published two-dimensional distributions for the diffuse ISM is the larger fractional polarization at 353 GHz for the denser molecular clouds. In the case of the molecular clouds, a decrease in the maximum polarization fraction with increasing column density is reported by Planck Collaboration Int. XX (2015). So we would expect for the molecular cloud sample a trend to lower p₃₅₃ values but we find the opposite. Our investigations in Sect. 2.2 show that there is in general only a weak correlation for orientation angles θ between FIR filaments and H I structures in column densities. Accordingly we find only a weak and probably insignificant correlation between H I column densities and polarization dispersions $\mathcal{S}$ and fractions p_HI.

Additionally, in simulations of MHD turbulence of polarized thermal emission from Galactic dust by Planck Collaboration Int. XX (2015), the polarization angle dispersion $\mathcal{S}$ is found to anticorrelate with the polarization fraction. These authors considered a lag of δ =16′, comparable to our resolution. The model assumption is that there is a large-scale anisotropic component of the magnetic field. Cold dense filaments and clumps are condensing of the magnetized WNM. In this case a turbulent magnetic field component is linked to turbulent velocity perturbations caused by a converging flow. The high H I polarization p_HI, as well as the absence of a correlation with $\mathcal{S}$ from our investigations are inconsistent with these studies. Also, the case for a strong turbulence, advocated by Planck Collaboration XII (2020, Appendix F.6) on basis of an anticorrelation between $\mathcal{S}$ an p at 353 GHz, isnot supported by our investigations if we consider filaments with coherence in FIR and H I.

The discrepancies between the polarization fractions derived from the H I Hessians and the Stokes parameters at 353 GHz indicate systematical differences caused by the data processing. The derived polarization is certainly affected by the relative orientation of the local magnetic field with respect to the line of sight. When comparing gas and dust we assume that both are affected in a similar way. But we need to take into account that different parts of the ISM are sampled. Figure 17 shows the relation between polarization measures p₃₅₃ and p_HI. For the full sample of filamentary structures considered by us the polarization derived from Stokes parameters is on average only 10% of the H I polarization. As detailed by Clark & Hensley (2019, Sect. 6.1), the observed polarization fraction is an important constraint on the intrinsic polarizing efficiency. In the case of dust grains, the size distributions and alignment functions are important, as is the local orientation of the magnetic field (Draine & Fraisse 2009). This leads to a limited observable polarization of ~ 22% (Planck Collaboration XII 2020). The H I polarization p_HI reflects onlythe internal coherence of the gas distribution and not the intrinsic polarization efficiency of dust. Clark & Hensley (2019) speculate that the most coherent sightlines may have p_HI ~ 1. At a resolution of 80′ they find p_HI ≤ 0.26, only slightly more than the observed maximum dust polarization of ~22% at 353 GHz. Our analysis results in observed dust polarizations that are comparable to previous investigations; however, we get far stronger H I polarizations, up to 80% and on average ~ 30%. It is not possible to explain high H I polarization fractions with noise or systematic effects. We interpret therefore the lower H I polarization fractions observed by Clark & Hensley (2019) with beam depolarization due to tangling within their larger beam. In a similar way the high polarization angle dispersions $\mathcal{S}$, derived by us from the Hessians, can be explained. The H I filaments are well resolved structures in narrow velocity intervals of 1 km s⁻¹, with a previously unreachable spatial resolution. Our smaller beam allows fluctuations in orientation angles to be observed on small scales.We emphasize that the Hessian analysis is most sensitive on such small scales.

Fig. 14 Two-dimensional density distribution functions for the polarization angle dispersion $\mathcal{S}$ along FIR filaments at a spatial resolution of 18′ as a function of the H I polarization fraction. Top: high Galactic latitudes. Bottom: all sky.

Fig. 15 Two-dimensional density distribution functions for the polarization angle dispersion $\mathcal{S}$ along FIR filaments at a spatial resolution of 18′ as a function of the polarization fraction derived from Stokes parameters at 353 GHz. Top: high Galactic latitudes. Bottom: all sky.

Fig. 16 Two-dimensional density distribution for the polarization angle dispersion $\mathcal{S}$ along FIR filaments in molecular clouds at a spatial resolution of 18′ as function of polarization fractions in H I filaments (top) and at 353 GHz (bottom). We used the same molecular cloud sample as considered by Planck Collaboration Int. XX (2015).

Fig. 17 Two-dimensional density distribution showing the relation between polarization fractions of filaments in H I and at 353 GHz.

3.6 Spatial distribution of coherent structures

Figure 18 shows on top the spatial distribution of the velocity field that we derive for the coherent FIR and H I filaments for the dominant velocities − 50 < v_LSR < 50 km s⁻¹, but to improve the presentation we display only − 30 < v_LSR < 30 km s⁻¹. Contributions outside this velocity range are statistically unimportant (Fig. 11). In the middle we present the distribution of the velocity dispersions $\mathcal{V}$ along filaments and below the harmonic mean Doppler temperatures in logarithmic scaling.

H I filaments are local structures that are embedded in the diffuse H I distribution. Derived harmonic mean Doppler temperatures according to Eq. (13) are therefore upper limits to Doppler temperatures for the H I gas that is directly associated with the FIR filaments. Figure 18 bottom shows a general trend that coherent FIR and H I filaments are cold, at least for high latitudes. Many structures close to the Galactic plane may be severely affected by confusion. The velocity dispersions $\mathcal{V}$ for the most prominent filaments are low. In some regions we find filaments with significant velocity gradients, indicating either dynamical interactions or overlap of unrelated filaments from different velocity regimes.

4 Filament curvature diagnostics

In Sect. 2.5 we determined FIR and H I filaments as coherent structures with nearly perfect angular alignment. This concordance is confined to the narrowest velocity interval and the continuity of the velocity field (Figs. 6 and 18) implies spatial coherence for FIR and H I structures in the plane of the sky. We adopt in the following the paradigm by Clark et al. (2014) and Clark & Hensley (2019) that such coherent H I structures can be used to probe the properties of the associated magnetic field, in particular the field tangling. We find that the FIR and H I coherence is confined to particular cold gas and is best defined for small-scale structures. The key for our understanding of the relations between gas, dust and magnetic fields is therefore the small-scale structure in narrow velocity intervals on arcmin scales despite the fact that such structures are part of larger and more prominent filaments (Kalberla et al. 2016) that appear to be more interesting upon visual inspection. Clark et al. (2014) and Clark & Hensley (2019) analyzed filamentary structures with a resolution of 60′ –80′ but we focus here on scales of 18′.

4.1 The filament curvature distribution

The small-scale structure of MHD turbulence is characterized by a folded structure of the fields, observable as transverse spatial oscillation of the field direction, while the field lines remain mainly unbent on large scales. Clark et al. (2014) and Clark & Hensley (2019) analyzed the line of sight field tangling. Our high resolution analysis allows us to focus on tangling as observed in projection perpendicular to the line of sight, hence changes of the orientation angle from θ₀ at the beam center to θ₁ at an offset position.

To parameterize the bending we need to consider that the magnetic field points in the direction of the flux tube’s tangent vector while a curvature is the response to a local perturbation perpendicular to the field line. We determine for each position with orientation angle θ₀ along the filaments the local curvature C = 1∕R for a structure with a radius determined from R = δ∕(2sin((θ₀ − θ₁)∕2)). We select offsets δ= 9′ with orientation angle θ₁ on both sides of the central position and average the radii. We verified that our results do not depend critically on the particular choice of δ. Figure 19 displays the distribution of curvatures calculated from FIR filaments derived at 857 GHz. Narrow channel H I filaments (see Figs. 7 top and middle for comparison) lead within the uncertainties to indistinguishable results. Our observational findings are quantitatively consistent with a power tail ∝ C^−13∕7 predicted by Schekochihin et al. (2002) for curvatures of field line structures generated by a small-scale turbulent dynamo. According to this model the small-scale magnetic turbulence is caused by folded structures of the magnetic field B and the geometry of the tangled magnetic fields can be studied in terms of the statistics of their magnetic tension forces or curvatures $K=\widehat{b}\cdot \nabla \widehat{b}$ with $\widehat{b}=B/B$.

Schekochihin et al. (2002) considered several model variants, the power tail ∝ C^−13∕7 remains valid in all cases for large curvatures. After an initially exponential growth the dynamo saturates and the fully developed, forced, isotropic MHD turbulence is expected to be the final saturated state of the small-scale dynamo. These authors derive with their Eq. (9) a curvature distribution $$P(K)\propto \frac{K}{{(1+K^2)}^{10/7}}$$(18)

with two asymptotes, P(K) ∝ K¹ for low curvatures and P(K) ∝ K^−13∕7 for high curvatures.

In Fig. 20 we compare our data after fitting the scaling parameter with this relation³. On top we display the curvature distributions at 857 GHz derived after application of the Hessian operator, using Eq. (3). Below we plot the curvature distributions calculated from 353 GHz Stokes parameters U and Q along the same filament positions using Eq. (4) to calculate orientation angles. Figure 20 shows systematical differences in the curvature distributions that are related to different analysis methods. Inherent to the application of the Hessian operator is a low sensitivity to structures larger than the applied pixel matrix, here 5 × 5 pixels or 1. °5 × 1. °5. The processing according to Eqs. (1)–(3) does however not bias significantly the calculation of the orientation angels compared to the straight forward calculation from Stokes parameters using Eq. (4). The shapes of the curvature distributions remain consistent with the Eq. (18).

To estimate limitations from observational uncertainties we repeated the calculations of the curvature distribution. To test the hypothesis that the $\mathcal{S}$ map could be noise dominated, we replace the orientation angles derived from Eqs. (3) and (4) with a cyclic uniform distribution. Repeating the complete data analysis results in a flat curvature distribution. For a random uniform distribution of orientation angles without any correlation of angles between neighboring positions all curvatures have the same probability. This is clearly inconsistent with observations. Our simulations include beam effects, resulting in a decreased sensitivity at high spatial frequencies. Accordingly all slopes in Fig. 20 at high curvatures are affected by sensitivity limitations. The predicted tail ∝ C^−13∕7 shown in Fig. 19 is not useful to constrain the model in this range. Comparing the 353 GHz and 857 GHz curvature distributions, we do not find any obvious noise biases. This is on the first glance surprising because the 353 GHz data are affected by noise much more severely than the 857 GHz ones but we show below that the S/N for the selected 353 GHz filaments is high enough.

We like to point out here first that the observed orientation angles for neighboring positions are correlated due to systematic field curvatures. The definition for $\mathcal{S}$ (Eq. (17)) does not account for such a systematic behavior. The distribution of curvatures is highly nonrandom, indicating that the filaments follow a systematical bending of the magnetic field lines. The peak at $\mathcal{S}~{52}^{°}$ in Fig. 13 is consistent with strongly fluctuating orientation angles on small scales without directional preferences as expected in the case of a tangled magnetic field. This peak does however not imply a completely random and uncorrelated distribution of orientation angles. Such a case requires the curvatures (derived from orientation angels at adjacent positions along the filaments) to also have a flat distribution.

We emphasize further that the derived distribution of observed local H I filament curvatures is far from noise dominated. Using Eq. (12) we determine the S/N at each of the selected filament positions for all H I channels with δv_LSR = 1 km s⁻¹ considered by us throughout this paper as related to the FIR. These data exceed an average S/N level of 87 all sky and 68 at high latitudes. The Hessian operator is selective and most sensitive to H I data with local S/N maxima (constant multiple rule). These are positions that are dominated by the CNM. The structures disclosed this way are predominantly observed close to zero velocity. Fortypical applications of a Hessian analysis a S/N of better than five is considered to be sufficient (e.g., Polychroni et al. 2013 or Soler et al. 2020, we refer also to the discussion in Sect. 2 of Schisano et al. 2014).

In a similar way we estimate the S/N for filaments at 353 GHz from the smoothed Planck intensity map. Following Planck Collaboration Int. XLVIII (2016) we use a 12. °5 × 12. °5 field at l = 90°, b = −80° to determine the background noise. Accordingly the 353 GHz FIR filaments used in our analysis exceed an average S/N level of 30 all sky and 18 at high latitudes. This is only about 1/4 of the S/N in H I but still high enough that we do not need to be worried about hidden noise biases in Fig. 20 and other parts of our analysis. Montier et al. (2015) consider uncertainties on polarization fraction and angle measurements of Planck polarization data and conclude that these measurements are little affected by uncertainties for an intensity S/N level above 10. Skalidis & Pelgrims (2019) determined the dependences between beam depolarization and smoothing radius of 353 GHz Planck data. They describe in their Appendix A that a smoothing radius of 20′ results in the best compromise to mitigate the beam depolarization and work with sufficiently high quality data.

Previous investigations of filamentary structures by Planck Collaboration Int. XIX (2015), Planck Collaboration Int. XXXVIII (2016), Clark (2018), Clark & Hensley (2019), and Planck Collaboration XII (2020) were focused mainly to angular resolutions between 80′ and 160′. On these scales only the very tip |C|≲1 deg⁻¹ of the curvature distribution discussed in the previous subsections is accessible but we extended this range by an order of magnitude. The agreement between our data and relation (18) is in general excellent for |C|≳1 deg⁻¹. We select only the range |C| > 1 deg⁻¹ to fit the observed curvatures to the distribution Eq. (18). The small-scale turbulent dynamo considered by Schekochihin et al. (2002, 2004) dominates the curvature distribution only for such scales.

For Kolmogorov turbulence the energy that feeds the turbulence is injected at an outer scale of ~ 100 pc (Kalberla & Haud 2019). The energy cascades then down to the viscous dissipation scale (called Kolmogorov inner scale in purely hydrodynamic systems). It is important to realize that the small-scale turbulent dynamo can feed energy in the opposite way. It acts initially below the viscous scale of ~ 10⁻² pc (Schekochihin et al. 2004) but perturbations in the magnetic field can propagate to larger scales. The dynamo leads initially to an exponential magnetic energy growth. A sketch of scale ranges and energy spectra is given in Sect. 1.2 of Schekochihin et al. (2004) with their Fig. 1. Folding of the magnetic fields causes magnetic tension forces with a back reaction on the turbulent flow. The system becomes increasingly nonlinear. Geometrically, one observes that stretched structures and flux sheets are converted to flux ribbons. In this so-called saturated stage, low curvatures are more dominant and the peak of the curvature distribution has shifted up and to small curvatures (see Fig. 25 of Schekochihin et al. 2004). For 353 GHz Stokes parameters (unaffected by biases from spatial filtering) the peak is at |C| ~ 0.4 deg⁻¹, corresponding to radii R≳2.5°. The kinetic energy on larger scales has to be fed in by external sources. It is expected that the dynamo saturates with an equilibrium between magnetic and kinetic energy. According to Schekochihin (2020, Sect. 12) details about the final state of the saturated small-scale dynamo remain however an unsolved problem, both numerically (due to lack of resolution) and theoretically (due to lack of theoreticians). Beresnyak (2012, his Fig. 1) argues that most of the magnetic energy may be concentrated at the scale corresponding to the peak in the curvature distribution. Supernovae are the dominant sources of turbulent energy in the ISM (Mac Low & Klessen 2004) and may play a major role in this case. The size distribution of H I shells in the outer Galaxy as potential sources for diameters larger than 1. °5 was studied by Ehlerová & Palouš (2013) and is in good agreement with the observed distribution of filaments at low and intermediate latitudes.

According to Schekochihin et al. (2004) the timescale for the growth of the small-scale magnetic energy is characterized by the viscous eddy turnover time in the order of 10⁵ yr. The mean-field dynamo theory predicts a large-scale galactic field exponentiating timescale in a rotation time of ~ 10⁸ yr. Thus the small-scale dynamo is very fast. These timescales may be compared to the typical cooling time of ~ 1.4 × 10⁵ yr for the CNM and ~5 × 10⁷ yr for the WNM. Shocks due to supernovae are expected to occur every ~ 5 × 10⁶ yr (Kalberla & Kerp 2009). In this situation the WNM can hardly be in equilibrium while phase transitions to the CNM can easily be triggered by a small-scale dynamo. In consequence, cold structures in the diffuse ISM that are driven by a small-scale dynamo, must evolve from small to large scales. Balsara & Kim (2005) studied the magnetic field amplification in the ISM by supernova-driven turbulence. They find a field amplification that increases with volume density as |B| ∝ ρ^0.386 and conclude that the field amplification takes place more vigorously in the lower temperature, denser gas.

Several authors have studied the role of the magnetic field in molecular cloud formation and evolution. These investigations have been reviewed by Hennebelle & Inutsuka (2019). From MHD simulations the consensus is that the magnetic field is strongly shaping the interstellar gas by generating a lot of filaments. For the diffuse H I distribution with column densities below N_H= 10^21.7 cm⁻² the filaments are aligned with the magnetic field. Recent investigations of phase transitions by Falle et al. (2020) emphasize that cooling behind a shock is affected by the magnetic field. Even a small initial magnetic field can lead to a magnetically dominated state on the unstable part of the equilibrium curve. These authors conclude that the magnetic field must dominate in the final state of shocks even for an implausibly small initial magnetic field strength and derive filament widths of ~ 0.52 pc (Figs. 11 and 12 of Falle et al. 2020).

The coherence between cold CNM and magnetized filamentary structures on small scales is explainable in the framework of these investigations. For the case of a magnetic pressure confinement of CNM filaments Kalberla et al. (2016) derived a median filament thickness of 0.09 pc assuming that the distance to the filaments is on average 100 pc. This estimate is consistent with Clark et al. (2014) who use better resolved Arecibo telescope observations and find largely unresolved structures, corresponding to a scale of 0.12 pc or below. Such structures have so far not been related to the small-scale dynamo. Shukurov & Sokoloff (2007) and Rincon (2019) consider the turbulent magnetic field as presumably confined to the WNM. If we, however, think about the simplified case of a flux tube compression in ideal MHD, we find that, in the case of phase transitions, the gas density and the magnetic field strength in such tubes should be amplified (Heiles & Crutcher 2005, Sect. 4.4). Phase transitions coupled to flux tubes may eventually lead to a balance between magnetic and gaseous pressure with typical parameters as used by Kalberla et al. (2016) in their Sect. 5.12 (see also Heiles & Troland 2005, Sect. 7.2, for discussion). For a magnetically confined CNM it is necessary to have a sufficient large magnetic field outside. Observational evidence for a flux tube scenario was reported by Kalberla et al. (2020, Sect. 3.3). In this context it may be worth to reconsider the equilibrium arguments given by Beresnyak (2012) by including local interactions from phase transitions.

Power spectra for the local CNM distribution show a remarkable power excess at multipoles l≳100, corresponding to scales below ~180°∕100 = 1.8° (Kalberla & Haud 2019). This is the range that we find to be occupied by the small-scale dynamo in Fig. 20. Characteristic for this excess is a notable shift to higher values of l when decreasing the width of the analyzed velocity slices or decreasing the Doppler temperature of the H I distribution under consideration (Fig. 15 of Kalberla et al. 2020). The power excess for the coldest structures is at multipoles l≳500, corresponding to scales of 180°∕500 ~ 20′ (see Fig. 15 of Kalberla et al. 2020).

The so far unexplained excess in turbulent CNM power spectra (Kalberla et al. 2020) is compatible with the transition between large and small-scale turbulent driving from our current analysis. On multipoles l≳100 the small-scale dynamo is dominant. Last but not least we need to take into account that λ₊, the second eigenvalue of the Hessian operator H from Eq. (2), discloses numerous unresolved structures along filaments at similar scales (see Fig. 1). These sources are also cold (Sect. 3.1.2) and belong therefore to the excess part of the exponential power distribution for the CNM. The coherence between FIR and H I filaments exists only for cold H I (Sect. 3.1.1) and turbulence is in this case affected by the small-scale dynamo. The injection of turbulent energy on small scales implies additional power at large multipoles and correspondingly a flattening of the CNM power spectra in comparison to multiphase H I power spectra. Changes of the spectral index depending on the phase composition of the H I (Kalberla & Haud 2019) are expected.

Fig. 18 All-sky Mollweide displays of derived H I parameters. Top: velocity distribution for coherent FIR and H I filaments in the velocity range − 30 < vLSR < 30 km s−1. Middle: distribution of the velocity dispersions along filaments for $0<\mathcal{V}<25$ km s−1. Bottom: spatial distribution of harmonic mean H I Doppler temperatures for 2 < log(TD) < 4.

Fig. 19 Distribution of local filament curvatures C from H I filaments and FIR at 857 GHz compared with a power tail ∝ C−13∕7 according to Schekochihin et al. (2002).

Fig. 20 Distribution of observed local filament curvatures |C| at a resolution of 18′ in comparison to curvatures |K| from Eq. (18) (Eq. (9) of Schekochihin et al. 2002). Top: curvatures for FIR filaments at 857 GHz and in H I derived from Eq. (3). Bottom: curvatures from 353 GHz Stokes parameters Q and U, using Eq. (4). As indicated by the vertical line, the fit was only applied to |C| > 1 deg−1.

4.2 Curvature versus filament strength

In the saturated state of the small-scale dynamo the magnetic field strength and the curvature are anticorrelated, |B| ∝ C^−1∕2 (Fig. 17 of Schekochihin et al. 2004, Figs. 5, 15, and 17 of Schekochihin et al. 2004, and Fig. 5 of St-Onge & Kunz 2018). Sharply curved fields (high C values or small radii) imply a high field tension and the field strength is reduced. The anticorrelation between field strength and curvature implies that the most prominent straight filaments popping up in FIR or H I maps must be indicating high values for the field strength and polarization fraction (little depolarization by tangling). For a prominent example we refer to a display of the Riegel–Crutcher cloud, Figs. 10 and 12 of Clark et al. (2014), as well as McClure-Griffiths et al. (2006).

To allow a general comparison between curvatures C and strength of filamentary structures we use the eigenvalues λ₋ as proxies for the strength of the magnetic field. The relation between λ₋ and the intensity of a filamentary structure, hence the strength of the magnetic field, is nonlinear. λ₋ depends on the second-order partial derivatives but we expect that the prominence of filaments is related to λ₋ (see Fig. 5 of Kalberla et al. 2016). Figure 21 displays the two-dimensional distribution of observed local filament curvatures C and eigenvalues − λ₋ on a logarithmic scale, verifying the expected relation. The most prominent filaments have low curvatures, leading to their description as fibers. Tangled worm-like structures have low intensities, are less prominent and have therefore previously mostly been discarded when discussing filamentary structures (e.g., by Kalberla et al. 2016). However, these less prominent worms belong definitively to the structures caused by a small-scale Galactic dynamo and should not be excluded from the scientific discussion of the ISM.

Fig. 21 Two-dimensional density distribution of observed local filament curvatures C at 857 GHz in comparison to eigenvalues λ− from the Hessian analysis as a measure of the filament strength.

Fig. 22 Two-dimensional density distribution of observed local filament curvatures C in comparison to H I polarization fraction. Top: high latitude sky. Middle: selected fields with strong H I anisotropies and associated radio-polarimetric filaments (Kalberla & Kerp 2016; Kalberla et al. 2017). Bottom: two-dimensional density distribution of spurious H I polarization derived by modeling a random distribution of orientation angles along the line of sight.

4.3 Curvature versus H I polarization fraction

A general relation between the curvature |C| and the strength of the magnetic field implies that also the H I polarization fraction p_HI must be related to C. Here we must take geometrical effects into account. The curvature |C| is derived from bending of filaments in the plane of the sky while the H I polarization fraction p_HI is defined along the line of sight. Small curvatures imply high magnetic field strengths with high internal coherence and H I polarization fractions. Figure 22 shows on top the observed relation between C and p_HI at latitudes |b| > 20°; all sky results are similar but not displayed. The highest polarization fractions are reached for filaments with the lowest curvatures as expected for a small-scale dynamo.

We repeat the calculations by modeling as in Sect. 4.1 a random distribution of orientation angles. We consider the same sample of positions along the filaments except that now the angles θ from Eq. (3) are replaced by a random distribution in velocity. The result is shown in the bottom panel of Fig. 22. We find nowhere an indication that a spurious polarization signal of this kind could affect our results, the orientation angles θ in our analysis are well defined

The FIR and H I coherence is not the only case that hints to relations between H I filaments and magnetism. Three fields that allow a detailed comparison of filamentary structures observable with the Low Frequency Array (LOFAR) as polarimetric filaments in close connection to cold H I filaments have been studied previously by Kalberla & Kerp (2016) and Kalberla et al. (2017). The LOFAR observations revealed strikingly linear coherent structures in Faraday depth, including some prominent filaments several degrees in length. The orientations of the LOFAR structures are affected by magnetic fields in the warm ionized medium. The observed alignment between magneto-ionic structures and H I filaments in these fields implies that the H I filaments are shaped by the magnetic field. The data from these polarimetric studies are displayed in the middle of Fig. 22 for comparison. This plot shows statistical uncertainties from the restricted number of analyzed filaments but there are no indications for any noise biases as documented in the lower panel. The two-dimensional density distribution of C against p_HI shows a similar trend as the distribution from FIR and H I structures at high latitudes, except that the polarization fractions in the case of the polarimetric filaments are somewhat lower. We question that the lower polarization is significant; a part of the filaments are close to the Galactic plane and may suffer from confusion. The investigations on these radio-polarimetric structures indicate anisotropies and systematic changes of turbulent power spectra for those parts of the H I velocity distribution that are associated with the magneto-ionic medium, we refer to Kalberla & Kerp (2016) and Kalberla et al. (2017) for details.

5 Curvature changes across flux tube envelopes

The determination of filament properties and curvatures in the previous sections was limited either to FIR filaments or to H I structures in narrow channels with best fit velocities v_fil as determined in Sect. 2.5. We adopt here the working hypothesis that these filaments stand for structures of magnetized flux tubes in the diffuse ISM. Coherent FIR and H I structures at positions with the best fit velocities v_fil (Fig. 18) represent in such a model the central parts, the bones of the flux tubes. We intend in this section to determine the properties of filament envelopes along the flux tubes at positions deviating from the central parts. We distinguish structures with positional offsets in the plane of the sky but also entities along the line of sight with velocities that differ from v_fil.

5.1 Structural changes at offset velocities

Curvatures that we discussed in the previous sections depend on positions and orientation angles that are defined only along the filaments at velocities v_fil. Considering deviating velocities Δv_LSR = |v_LSR − v_fil| we may expect that on average structures with increasing offsets Δv_LSR are also offset increasingly in distance. For a turbulent medium the characteristic velocities v_l and v₀ at scales l and l₀ are related according to $v_l/v_0={(l/l_0)}^{1/3}$ (e.g., Frisch 1996, Sect. 7.6). Considering H I channels at velocities offsets from v_fil would accordingly allow the flux tube properties far from the central parts to be determined, though an exact scaling between velocity and distance is missing.

We extend the calculation of curvatures to all individual H I channels and consider at each filament position also the velocities v_fil ± Δv_LSR for |Δv_LSR| < 25 km s⁻¹. To avoid any noise biases we use identical criteria as in Sect. 2.5 for the significance of filamentary structures. We characterize systematical changes in curvature between envelope and center of the flux tube by the ratio R_C (Δv_LSR) = C(v_fil + Δv_LSR)∕C(v_fil). Figure 23 displays two-dimensional histograms of curvatures C and curvature ratios R_C at velocity offsets of Δv_LSR = 1, 4, and 16 km s⁻¹. We observe the general trend that for increasing velocity offsets Δv_LSR the curvature ratios R_C increase. In other words, the filament curvatures increase systematically when considering outer filament layers that are characterized by increasing velocity offsets |Δv_LSR|. Considering the prognosis that the magnetic field strength and curvature are for a small-scale dynamo anticorrelated, |B| ∝ C^−1∕2 (Schekochihin et al. 2004; St-Onge & Kunz 2018), we conclude that the magnetic field of the low curvature segments in the centers of the flux tubes must be enhanced relative to the more curled environments. These central parts are most significant (see Sect. 4.2). The eigenvalues λ₋ increase systematically with increasing velocity offsets and at |Δv_LSR|~ 4 km s⁻¹ the number of significant filamentary structures that can be analyzed drop by about 50%. At |Δv_LSR|≳16 km s⁻¹ only a few significant filamentary H I structures remain.

5.2 Structures offset in position

Changes in filamentary structures for narrow H I channels can also be observed at positions offset from the central parts of the filaments. However, the disadvantage of such an analysis is that we can no longer track curvatures along the bones of the filaments. Curvatures in the previous sections were determined along the filament ridges but for offset positions this orientation is lost. We can only determine average local curvatures that can be related to the curvature at the closest filament position with distance $D~\sqrt{{((l-l_f)\cos (b))}^2+{(b-b_f)}^2}$. It also can happen that the most significant structure with the lowest λ₋ value at the position l, b is found at a velocity that deviates from v_fil at the nearest position l_f, b_f. Accordingly we need to change definition for the curvature ratios to R_C(D) = C(l, b, v_LSR)∕C(l_f, b_f, v_fil). In Fig. 24 we display these ratios for H I layers in the velocity ranges |v_LSR − v_fil|≤ 1, 4, and 16 km s⁻¹. R_C (D) can in most cases only be traced to D≲1°. Searching for related filamentary structures up to distances D ~ 10° needed lots of CPU time but remained without success.

The plots in Fig. 24 are not directly comparable to those in Fig. 23 but they tell the same. Curvatures for structures offset in velocity-position space from the central parts of flux tubes increase significantly in comparison to curvatures along central bones of the filaments. Increased curvatures imply increased curling of magnetic field lines with decreasing field strengths according to predictions for the small-scale dynamo. Such structures have lower H I intensities and the λ₋ values are less significant.

5.3 Structures in H I column densities and FIR

We repeated the analysis from the previous subsection but substituted the individual channel maps with a single H I column density map for |v_LSR| < 50. The result is shown in Fig. 25 on the left. In this case we find no evidence for systematic changes of the curvature distributions with significant preferences in curvature radius as observed in Figs. 23 and 24. Averaging across numerous filaments smears out any information on the curvature of the ISM. Structures along the line of sight but outside the filaments are weaker, accordingly tangling causes de-correlation when integrating in velocity. The same effect occurs if we consider the distribution of FIR intensities. Figure 25 shows that the curvature distributions in FIR (right) replicate the distributions for N_HI on the left hand side.

As discussed in Sect. 2.6, a coherent velocity field along the filaments is only detectable if we recognize that the FIR filaments are associated with H I structures on small scales. We demonstrated in Sect. 3 that these features are cold. Accordingly they need to be analyzed with narrow channel spacings. Integrating in velocity smears out structures caused by the CNM. The differences between Figs. 23 and 24 on one side and Fig. 25 on the other side demonstrate impressively observer’s obligation to obey the sampling theorem (Blackman & Tukey 1958; Shannon & Weaver 1975). Otherwise sensitive information about cold turbulent structures on small scales, imprinted on the velocity field, gets lost. This information is available in full resolution three-dimensional spectral H I data with sufficient sensitivity but not in two-dimensional Planck FIR observations or H I data that are missing the third dimension. The necessity to consider the three-dimensional distribution in the ISM was also demonstrated recently by Pelgrims et al. (2021) when studying the line of sight frequency decorrelation of polarized dust emission along the line of sight.

Fig. 23 Two-dimensional histograms of filament curvatures C in central parts of the filaments versus changes in curvature radius, characterized by the ratio RC for channels offset from vfil by 1, 4, and 16 km s−1 (top to bottom). Displayed are the all sky data (left) and at the high latitude data (right).

6 Summary and discussion

Filamentary H I structures are known to be well aligned with the plane-of-sky magnetic field orientation as measured with optical starlight polarization and polarized thermal dust emission (e.g., Clark 2018). This tight correlation exists for narrow velocity intervals only, and the orientation of these features is found to be coherent in velocity space. Clark (2018) and Clark & Hensley (2019) propose using this H I coherence to define for such structures an H I polarization along the line of sight. Analogous to optical polarization that can be produced by transition through aligned structures in refracting media, this H I polarization is indicative of coherence and alignment in H I filaments. It was shown by Clark & Hensley (2019) that H I Stokes parameter maps are comparable to the Planck 353 GHz Q and U maps of polarized dust emission. At a resolution of 80′, many of the large-scale features are reproduced in H I.

Here we study the FIR and H I coherence in more detail with a better spatial resolution. We consider the more sensitive Planck 857 GHz FIR data smoothed to 18′ resolution. For the H I we use unsmoothed HI4PI data at resolutions of 10.′8 for EBHIS and 14.′5 for the GASS survey. Using the Hessian operator, we extract filaments simultaneously for both data sets and determine alignment angles along the filaments. We consider several measures for the alignment of the filamentary structures, and these indicate in all cases that the alignment on small scales improves significantly with respect to previous determinations on scales of around one degree. For H I in narrow velocity intervals of 1 km s⁻¹, we obtain an excellent alignment between FIR and H I filaments, indicating a close coherence between FIR and CNM. These filamentary structures are also coherent in velocity space and, in general, are cold. Following Clark (2018) and Clark & Hensley (2019), we determine in Sect. 3.3 the H I polarization along the line of sight and find for filamentary structures an average H I polarization fraction of 30% with peak values around 80%. This is a factor of two to three more than reported in previous determinations. Since we did not apply any smoothing to the H I data, we probe, at an assumed distance of 100 pc (e.g., Sfeir et al. 1999), regions of about 0.3 pc.

We probed the velocity structure of the H I counterparts for the FIR filaments and discovered a well-defined coherent velocity field. Filaments are local phenomena with center velocities around 0 km s⁻¹. The FWHM width of this distribution is 16.6 km s⁻¹. Along the filaments the typical velocity dispersion is $\mathcal{V}=5.5$ km s⁻¹. Most of the H I filaments are associated with the CNM. Toward the structures investigated here we also observe Doppler temperatures that are lower by up to a factor of six for the accompanied LNM and WNM gas, implying a CNM in cold cores of a multiphase medium in phase transition. Low Doppler temperatures imply that FIR and H I coherence is best defined in narrow velocity intervals and probably also in small and dense regions in the centers of the filaments. The spatial distribution of the velocity field in comparison to velocity dispersions $\mathcal{V}$ and harmonic mean Doppler temperatures shows common structures for most of the prominent filaments, indicating that velocity dispersions and Doppler temperatures are correlated along these filaments.

The polarization angle dispersion $\mathcal{S}$ introduced in Planck Collaboration Int. XIX (2015) measures fluctuations of the polarization angle along FIR filaments. The total angular scale probed by us is four times the 18′ FIR resolution, hence five to six times the H I beam, or about 2 pc at an assumed distance of 100 pc. At such a scale we find, for a large fraction of the positions, significant fluctuations in the position angles relative to the center position. Contrary to previous investigations by Clark & Hensley (2019), we find no indications for a negative correlation between $\mathcal{S}$ and the polarization fraction p_HI. We confirm, however, the inverse correlation between $\mathcal{S}$ and p₃₅₃ at 353 GHzfrom Stokes parameters, derived previously mostly on angular scales between 80′ and 160′ (Planck Collaboration Int. XIX 2015; Planck Collaboration Int. XXXVIII 2016; Clark 2018; Clark & Hensley 2019; Planck Collaboration XII 2020).

The high spatial resolution of our analysis allows us to parameterize the bending of filamentary structures in terms of a curvature distribution P(C). We find a distribution that is well approximated by curvatures, as predicted in the framework of a small-scale turbulent dynamo by Schekochihin et al. (2002). From our data we derive an excellent agreement with the predicted curvature distribution (Schekochihin et al. 2002, 2004). This agreement is, however, limited to curvatures |C|≳1 deg⁻¹, corresponding to tangling radii of R≲1°. While the small-scale turbulent dynamo acts below this scale, the turbulent energy on the largest scales must be fed by external sources, most probably by supernovae. We confirm the prediction that for a small-scale dynamo field strengths and curvatures are anticorrelated (Schekochihin et al. 2002, 2004). The most prominent parts of the filaments have the lowest curvatures, and, accordingly, these structures should have the strongest fields. The curling in envelopes around the filaments increases significantly, implying a decay of the field outside the central parts of the flux tubes. Furthermore, we find that curvature is in general correlated with filament coherence along the line of sight, characterized by the H I polarization fraction.

Our results on coherent FIR and H I structures are for |C|≳1 deg⁻¹ incompatible with simulations of anisotropic MHD turbulence in several of the fields selected by Planck Collaboration Int. XX (2015). Throughout those simulations there is a large-scale anisotropic component of the magnetic field, as well as a turbulent component linked to the velocity perturbations imposed on converging flows. From our results we conclude that the spatial coherence of magnetized FIR and H I filaments is shaped by magnetic fields that were amplified by a small-scale turbulent dynamo. Turbulence emerging from small scales can change the spectral index of the turbulent flow at those velocities, which are characteristic for FIR and H I coherence. In the framework of a small-scale turbulent dynamo, it is also easy to support, on a Galactic scale, phase transitions in a magnetically dominated state, as considered recently by Falle et al. (2020).

One of the main conclusions by Schekochihin et al. (2004) is that the fully developed, forced, isotropic MHD turbulence is the saturated state of the small-scale dynamo. In this case the kinetic energy spectrum is dominated by the outer scale and has a steeper-than-Kolmogorov scaling in the inertial range, while the magnetic energy is dominated by small scales, at which it substantially exceeds the kinetic energy. The strongest part of the local Galactic H I emission extends over large fractions of the sky at velocities close to zero. Assuming that this emission is representative of a fully developed MHD turbulence driven by a large-scale dynamo, we expect steep power spectra for this part of the multiphase H I. Focusing our attention on the coldest part of the CNM, with volume filling factors at a 5% level (e.g., Murray et al. 2020), we find coherence between H I structures and the distribution of magnetic field lines. We find evidence that this range is dominated by the magnetic energy and the small-scale dynamo. Cold neutral medium power spectra are expected to be shallow in comparison to the large-scale, dominating part of the warmer H I distribution with excess power from structures on small scales.

In addition to filamentary structures, we find local condensations – numerous blobs – located along the filaments. These are cold coherent FIR and H I structures, point-like at 18′ resolution, but they share the properties of the filaments. This additional cold small-scale emission is consistent with the excess power observed in power spectra of the CNM at the highest spatial frequencies (Kalberla & Haud 2019). Since filaments are in general cold and restricted locally to small scales, it is also plausible that the CNM power spectra are shallower compared to multiphase power spectra since the power at high spatial frequencies is increased by the small-scale dynamo. Cold structures, driven by a small-scale dynamo, must evolve from small to large scales. Turbulence in the local diffuse ISM has to be related to phase transitions, but observational evidence currently appears to be disregarded in many theoretical investigations (Clark et al. 2019).

Fig. 24 Two-dimensional histograms of filament curvatures C in central parts of the filaments versus changes in curvature radius, characterized by the ratio RC for curvatures at the closest positions off the filaments. We integrate RC for velocitiesat channels |vLSR − vfil| by 1, 4, and 16 km s−1 (top to bottom). Displayed are the all sky data (left) and the high latitude data (right).

Fig. 25 Two-dimensional histograms of filament curvatures C in central parts of the filaments versus changes in curvature radius, characterized by the ratio RC for curvatures at the closest positions off the filaments. To the left we display ratios RC derived fromH I column densities, and to the right RC ratios are given for FIR intensities at 857 GHz. Top panels display distributions restricted to high Galactic latitudes. Bottom panels: from all sky data.

Acknowledgements

We acknowledge the referee for a careful reading and critical comments that helped to improve the quality of the manuscript. HI4PI is based on observations with the 100-m telescope of the MPIfR (Max-Planck-Institut für Radioastronomie) at Effelsberg and the Parkes Radio Telescope, which is part of the Australia Telescope and is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO. This research has made use of NASA’s Astrophysics Data System. Some of the results in this paper have been derived using the HEALPix package.

References

All Tables

All Figures

	Fig. 1 All-sky Mollweide displays of derived FIR structures at 857 GHz. Top: distribution of eigenvalues λ− that characterize filamentary structures. Bottom: distribution of λ+, characterizing local enhancements (blobs) along filaments.
In the text

	Fig. 2 Histograms of angular alignment deviations for filamentary structures according to Eq. (5). We compare 857 and 545 GHz data, unsmoothed as observed and smoothed to a resolution of 18′ FWHM.
In the text

Fig. 3 Histograms of angular alignment deviations according to Eq. (5) for filamentary structures. Top left: Planck 857 GHz compared with Planck 353 GHz and angular orientations according Eq. (4), |b| > 20°. Top right: Planck 857 GHz compared with H I column densities for |vLSR| < 50 km s−1, |b| > 20°. Center left: Planck 857 GHz compared with best fit single channel H I filaments, all sky. Center right: Planck 857 GHz compared with best fit single channel H I filaments, |b| > 20°. Bottom left: Planck 857 GHz compared with best fit H I filaments at a channel width of 16 km s−1, all sky. Bottom right: Planck 857 GHz compared with best fit H I filaments at a channel width of 16 km s−1, |b| > 20°.

In the text

	Fig. 4 Distribution of H I velocity dispersions $\mathcal{V}$ along FIR filaments at channel widths of 1, 8, and 16 km s−1.
In the text

	Fig. 5 H I velocity dispersions $\mathcal{V}$ along FIR filaments at a channel width of 1 km s−1 depending on simulatednoise enhancements Nf = 2, 4, and 8 of thetypical observed noise level.
In the text

Fig. 6 Velocity field in the case of H I filaments derived for data with a 1 km s−1 channel width at an arbitrarily selected central position, l = 160°, b = 30° in gnomonic projection. The field size is 27. °7; only velocities in the range − 30 < vLSR < 30 km s−1 are displayed. Top: velocities from original telescope data. Bottom: results derived after simulating a noise amplification by a factor of Nf = 4.

In the text

	Fig. 7 Orientation angles according to Eq. (3). Top: from the FIR distribution at 857 GHz.Middle: from the best fit H I distribution. Bottom: distribution of orientation angles according to Eq. (4) calculated from Stokes parameters U and Q of the 353 GHz FIR distribution.
In the text

	Fig. 8 Distribution of harmonic mean H I Doppler temperatures in coherent FIR and H I filaments.
In the text

	Fig. 9 Distribution of H I Doppler temperatures of coherent FIR and H I filaments compared with a Gaussian three-component fit. Top: high latitude data, |b| > 20°. Bottom: all sky.
In the text

	Fig. 10 Distribution of H I Doppler temperatures of local enhancements (blobs) along filaments and their approximations by Gaussian distributions (green). To distinguish the distributions from H I phases, we separate the CNM (blue) from the warmer phase, LNM and WNM (red).
In the text

	Fig. 11 Histograms of the velocity distributions for filamentary structures and their approximations by Gaussian distributions (green). To distinguish the distributions from H I phases, we separate the CNM (blue) from the warmer phase, LNMand WNM (red).
In the text

	Fig. 12 Distribution of the H I polarization pHI according to Eq. (16), all sky and for high latitudes at |b| > 20°.
In the text

Fig. 13 Distribution of the polarization angle dispersion $\mathcal{S}$ along FIR filaments at a spatial resolution of 18′. The peak at $\mathcal{S}={52}^{°}$, annotated with a vertical line, indicates that many positions have significant deviations of polarization angles along the filaments. The black and red lines mark filaments with eigenvalues λ− < −1.5 K deg−2, and blue andorange lines represent filaments limited to λ− < −3 K deg−2.

In the text

	Fig. 14 Two-dimensional density distribution functions for the polarization angle dispersion $\mathcal{S}$ along FIR filaments at a spatial resolution of 18′ as a function of the H I polarization fraction. Top: high Galactic latitudes. Bottom: all sky.
In the text

	Fig. 15 Two-dimensional density distribution functions for the polarization angle dispersion $\mathcal{S}$ along FIR filaments at a spatial resolution of 18′ as a function of the polarization fraction derived from Stokes parameters at 353 GHz. Top: high Galactic latitudes. Bottom: all sky.
In the text

	Fig. 16 Two-dimensional density distribution for the polarization angle dispersion $\mathcal{S}$ along FIR filaments in molecular clouds at a spatial resolution of 18′ as function of polarization fractions in H I filaments (top) and at 353 GHz (bottom). We used the same molecular cloud sample as considered by Planck Collaboration Int. XX (2015).
In the text

	Fig. 17 Two-dimensional density distribution showing the relation between polarization fractions of filaments in H I and at 353 GHz.
In the text

	Fig. 18 All-sky Mollweide displays of derived H I parameters. Top: velocity distribution for coherent FIR and H I filaments in the velocity range − 30 < vLSR < 30 km s−1. Middle: distribution of the velocity dispersions along filaments for $0<\mathcal{V}<25$ km s−1. Bottom: spatial distribution of harmonic mean H I Doppler temperatures for 2 < log(TD) < 4.
In the text

	Fig. 19 Distribution of local filament curvatures C from H I filaments and FIR at 857 GHz compared with a power tail ∝ C−13∕7 according to Schekochihin et al. (2002).
In the text

Fig. 20 Distribution of observed local filament curvatures |C| at a resolution of 18′ in comparison to curvatures |K| from Eq. (18) (Eq. (9) of Schekochihin et al. 2002). Top: curvatures for FIR filaments at 857 GHz and in H I derived from Eq. (3). Bottom: curvatures from 353 GHz Stokes parameters Q and U, using Eq. (4). As indicated by the vertical line, the fit was only applied to |C| > 1 deg−1.

In the text

	Fig. 21 Two-dimensional density distribution of observed local filament curvatures C at 857 GHz in comparison to eigenvalues λ− from the Hessian analysis as a measure of the filament strength.
In the text

Fig. 22 Two-dimensional density distribution of observed local filament curvatures C in comparison to H I polarization fraction. Top: high latitude sky. Middle: selected fields with strong H I anisotropies and associated radio-polarimetric filaments (Kalberla & Kerp 2016; Kalberla et al. 2017). Bottom: two-dimensional density distribution of spurious H I polarization derived by modeling a random distribution of orientation angles along the line of sight.

In the text

	Fig. 23 Two-dimensional histograms of filament curvatures C in central parts of the filaments versus changes in curvature radius, characterized by the ratio RC for channels offset from vfil by 1, 4, and 16 km s−1 (top to bottom). Displayed are the all sky data (left) and at the high latitude data (right).
In the text

	Fig. 24 Two-dimensional histograms of filament curvatures C in central parts of the filaments versus changes in curvature radius, characterized by the ratio RC for curvatures at the closest positions off the filaments. We integrate RC for velocitiesat channels |vLSR − vfil| by 1, 4, and 16 km s−1 (top to bottom). Displayed are the all sky data (left) and the high latitude data (right).
In the text

Fig. 25 Two-dimensional histograms of filament curvatures C in central parts of the filaments versus changes in curvature radius, characterized by the ratio RC for curvatures at the closest positions off the filaments. To the left we display ratios RC derived fromH I column densities, and to the right RC ratios are given for FIR intensities at 857 GHz. Top panels display distributions restricted to high Galactic latitudes. Bottom panels: from all sky data.

In the text