Turbulent power distribution in the local interstellar medium

TL;DR

This study analyzes the angular power distribution of the local interstellar medium using HI survey data, revealing phase-dependent turbulence characteristics and hierarchical structures across different gas phases.

Contribution

It provides the first detailed analysis of the power spectra for different HI phases in the local ISM, highlighting phase-specific turbulence properties and correlations.

Findings

Power spectra follow a power law with indices ≥ -3 for local gas.

Cold filamentary structures contribute excess power at small scales.

Different phases exhibit distinct but correlated distribution properties.

Abstract

Context: The interstellar medium (ISM) on all scales is full of structures that can be used as tracers of processes that feed turbulence. Aims: We used HI survey data to derive global properties of the angular power distribution of the local ISM. Methods: HI4PI observations on an nside = 1024 HEALPix grid and Gaussian components representing three phases, the cold, warm, and unstable lukewarm neutral medium (CNM, WNM, and LNM), were used for velocities $|v_{\mathrm{LSR}}| \leq 25$ \kms. For high latitudes $|b| > 20\deg$ we generated apodized maps. After beam deconvolution we fitted angular power spectra. Results: Power spectra for observed column densities are exceptionally well defined and straight in log-log presentation with 3D power law indices $\gamma \geq -3$ for the local gas. For intermediate velocity clouds (IVCs) we derive $\gamma = -2.6$ and for high velocity clouds (HVCs) $\gamma = -2.0$. Single-phase power distributions for the CNM, LNM, and WNM are highly correlated and shallow with $ \gamma \sim -2.5$ for multipoles $l \leq 100$. Excess power from cold filamentary structures is observed at larger multipoles. The steepest single-channel power spectra for the CNM are found at velocities with large CNM and low WNM phase fractions. Conclusions: The phase space distribution in the local ISM is configured by phase transitions and needs to be described with three distinct different phases, being highly correlated but having distributions with different properties. Phase transitions cause locally hierarchical structures in phase space. The CNM is structured on small scales and is restricted in position-velocity space. The LNM as an interface to the WNM envelops the CNM. It extends to larger scales than the CNM and covers a wider range of velocities. Correlations between the phases are self-similar in velocity.

Equations

$$P(l)=C_{\mathrm{l}}(l)\propto l^{\gamma },$$

$$P_{\mathrm{obs}}(l)=P(l)\cdot B_{\mathrm{data}}^2(l)+N_{\mathrm{oise}}(l).$$

$$P(l)=P_{\mathrm{obs}}(l)/B_{\mathrm{data}}^2(l).$$

$$P_{\mathrm{TOT}}=P_{\mathrm{WNM}}+P_{\mathrm{LNM}}+P_{\mathrm{CNM}}+2[P_{\mathrm{WxC}}+P_{\mathrm{WxL}}+P_{\mathrm{LxC}}].$$

$$f=\mathrm{lo}{\mathrm{g}}_{10}(N_{\mathrm{Hobs}}/10^{20})(0.25\pm 0.02)+(0.87\pm 0.02),$$

$$f_{\mathrm{P}}{(}_{v1}^{v2})=\frac{{\int }_{v1}^{v2}{T}_{\mathrm{bP}}(v_{\mathrm{LSR}})\delta v_{\mathrm{LSR}}}{{\int }_{v1}^{v2}{T}_{\mathrm{b}}(v_{\mathrm{LSR}})\delta v_{\mathrm{LSR}}},$$

$$\mathfrak{R}(l)=P_{\mathrm{ON}}(l)/P_{\mathrm{OFF}}(l)$$

$$\mathfrak{E}(l)=P_{\mathrm{ON}}(l)-P_{\mathrm{OFF}}(l)$$

$$S_l(v_1,v_2)=C_l^2(v_1,v_2)/[C_l(v_1,v_1)\cdot C_l(v_2,v_2)],$$

$$P_l(\mathrm{A},\mathrm{B})=C_l^2(\mathrm{AxB})/[C_l(\mathrm{AxA})\cdot C_l(\mathrm{BxB})]$$

$${\mathfrak{R}}_{\mathrm{VCA}}(l)=P_{\mathrm{broad}}(l)/P_{\mathrm{narrow}}(l),$$

$$\widetilde{I_{\mathrm{obsDump}}}=\widetilde{B_{\mathrm{tel}}}\ast \widetilde{I_{\mathrm{dump}}}+\widetilde{N_{\mathrm{dump}}}.$$

$$\widetilde{I_{\mathrm{obsHEAL}}}=\widetilde{B_{\mathrm{grid}}}\ast (\widetilde{B_{\mathrm{tel}}}\ast \widetilde{I_{\mathrm{dump}}}+\widetilde{N_{\mathrm{dump}}})=\widetilde{B_{\mathrm{data}}}\ast \widetilde{I_{\mathrm{sky}}}+\widetilde{B_{\mathrm{grid}}}\ast \widetilde{N_{\mathrm{dump}}}.$$

$$P_{\mathrm{obs}}=B_{\mathrm{data}}^2{P}_{\mathrm{sky}}+B_{\mathrm{grid}}^2{P}_{\mathrm{Noise}}+2P_{\mathrm{cross}},$$

$$P_{\mathrm{sky}}=P_{\mathrm{obs}}/B_{\mathrm{data}}^2-P_{\mathrm{Noise}}{B}_{\mathrm{grid}}^2/B_{\mathrm{data}}^2-2P_{\mathrm{cross}}/B_{\mathrm{data}}^2$$

$$P=P_{\mathrm{obs}}/B_{\mathrm{data}}^2-P_{\mathrm{Noise}}/B_{\mathrm{tel}}^2-2P_{\mathrm{cross}}/B_{\mathrm{data}}^2$$

$$T_{\mathrm{Noise}}(v_{\mathrm{LSR}})={\sigma }_{\mathrm{av}}[T_{\mathrm{sys}}+T_{\mathrm{B}}(v_{\mathrm{LSR}})]/T_{\mathrm{sys}},$$

$$Var(\widehat{C_l})/[N_l+C_l{]}^2=2/(2l+1).$$

Full text

*Aims. *We used H i survey data to derive global properties of the angular power distribution of the local ISM.

Methods. HI4PI observations on an nside = 1024 HEALPix grid and Gaussian components representing three phases, the cold, warm, and unstable lukewarm neutral medium (CNM, WNM, and LNM), were used for velocities $|v_{\mathrm{LSR}}|\leq 25$ km s-1. For high latitudes |b|>20\hbox{{}^{\circ}} we generated apodized maps. After beam deconvolution we fitted angular power spectra.

*Results. * Power spectra for observed column densities are exceptionally well defined and straight in log-log presentation with 3D power law indices $\gamma\geq-3$ for the local gas. For intermediate velocity clouds (IVCs) we derive $\gamma=-2.6$ and for high velocity clouds (HVCs) $\gamma=-2.0$. Single-phase power distributions for the CNM, LNM, and WNM are highly correlated and shallow with $\gamma\sim-2.5$ for multipoles $l\leq 100$. Excess power from cold filamentary structures is observed at larger multipoles. The steepest single-channel power spectra for the CNM are found at velocities with large CNM and low WNM phase fractions.

*Conclusions. *The phase space distribution in the local ISM is configured by phase transitions and needs to be described with three distinct different phases, being highly correlated but having distributions with different properties. Phase transitions cause locally hierarchical structures in phase space. The CNM is structured on small scales and is restricted in position-velocity space. The LNM as an interface to the WNM envelops the CNM. It extends to larger scales than the CNM and covers a wider range of velocities. Correlations between the phases are self-similar in velocity.

1 Introduction

Turbulence is ubiquitous in the interstellar medium (ISM). Neutral hydrogen is abundant and because of the easily observable H i 21cm emission line it is one of their key diagnostics. However, there are more bearings than just diagnostics. Already in the 1950s, soon after the detection of the 21cm line, Hendrik van de Hulst had “the hope that you could do something by making a turbulence spectrum, a new trick to describe things.”111Interview of Hendrik van de Hulst on 6 September 1978 at Leiden by W. T. Sullivan III, https://www.nrao.edu/archives/Sullivan/sullivan_interviewee_vandehulst.shtml van de Hulst was interested in measuring turbulence of the H i distribution to gain insights into the nature of the ISM. Today we have a better background to formulate the question of how far turbulence is shaping the ISM and in particular the bistable H i distribution. Turbulent motions determine the temporal and spatial structure of the gas pressure. If pressure fluctuations are sufficiently large, they drive some of the gas into the cold phase. Thus, the H i phase composition is affected by turbulence. We should also expect that changes in the H i phase composition have a noticeable imprint on the observable turbulent density and velocity distribution and focused on the H i multiphase structure.

For an introduction to topics addressed in the first paragraph we refer to the reviews about interstellar turbulence by Elmegreen & Scalo (2004) and Scalo & Elmegreen (2004). Our contemporary understanding of the multiphase structure of the ISM is based on seminal papers by McKee & Ostriker (1977) and Wolfire et al. (2003). Heating and cooling processes invoke thermal instabilities that tend to segregate the H i into two distinct stable phases, the cold medium (CNM) and the warm medium (WNM). For at least two decades there has been mounting evidence that a significant fraction of the H i is in an intermediate unstable state, the lukewarm medium (LNM); for the most recent census of the phase fractions we refer to Murray et al. (2018b) and Kalberla & Haud (2018). Turbulence is considered as a driving mechanism that tends to produce strong nonlinear fluctuations in all the thermodynamic variables that are affecting thermal instabilities. Audit & Hennebelle (2005) find in this case large fractions of thermally unstable gas that would not exist without turbulent forcing. These fractions increase with the amplitude of the turbulent forcing and the cold and thermally unstable gas tends to be organized in filamentary structures. As a result, the standard two-phase model may need to be replaced by a phase continuum with at least three phases; we refer to the excellent review by Vázquez-Semadeni (2012) and references therein.

Analyzing Arecibo data, cold filamentary H i structures that are aligned with the magnetic field have been found by Clark et al. (2014) at high Galactic latitudes. Similar large-scale structures all over the sky correlated with the magnetic field orientation implied by *Planck * 353 GHz polarized dust emission were reported by Kalberla et al. (2016). Some of the observed anisotropic H i structures show a clear association with magneto-ionic features (Kalberla & Kerp (2016), (Kalberla et al., 2017), and (Jelić et al., 2018)). Derived power spectra for the H i are in this case consistent with Kolmogorov turbulence, but anisotropies in narrow velocity intervals increase on average with spatial frequency, both as predicted by Goldreich & Sridhar (1995) for incompressible magnetohydrodynamic (MHD) turbulence. These observations strongly support an MHD origin of the turbulence similar to the conception brought forward by Heiles & Troland (2005) and Heiles & Crutcher (2005) who argued for an equipartition between turbulent and magnetic field energy.

H i observations are typically organized as data cubes in position-position-velocity (PPV). Power spectra of emission from narrow velocity channels are affected by velocity fluctuations that cause shallower slopes than spectra derived from broad channels. This is the result of an excess of small features from unrelated structures that blend by velocity crowding on the line of sight. The basic recipes used to analyze these data and disentangle velocity and density fields, called velocity channel analysis (VCA), were given by Lazarian & Pogosyan (2000). Velocity fluctuation can in principle mimic density structures, an effect described as velocity crowding by Burton (1972) and Lazarian & Pogosyan (2000). Recently, a vivid debate has been raised about structures seen in H i channel maps. Lazarian & Yuen (2018) interpret cold filamentary structures observed by Clark et al. (2014) as being caused by velocity caustics. Clark et al. (2019) object and interpret these filaments as genuine H i density structures that are associated with dust. Yuen et al. (2019) reinforce arguments given by Lazarian & Yuen (2018).

In our analysis we make use of a Gaussian decomposition of the HI4PI survey (Kalberla & Haud, 2018). We determine discrete spatial power spectra for CNM, LNM, and WNM at multipoles $l\la 1023$. When separating H i phases we discover limitations on VCA, constraints specified previously by Kolmogorov (1941) as restrictions on locally homogeneous and isotropic structures. The filamentary features debated by Lazarian & Yuen (2018), Clark et al. (2019), and Yuen et al. (2019) are identified as coherent CNM structures, resulting from phase transitions.

The observed cold filamentary H i column density structures are likely fibers with cylindrical geometry (Clark et al., 2014) or projections from edge-on sheets and have velocity gradients perpendicular to the sheets. Kalberla & Kerp (2016) and Kalberla et al. (2017) report in addition a steepening of the spectral indices at velocities that are dominated by cold gas, indicating phase transitions caused by thermal instabilities. This steepening is opposite to predictions from numerical simulations by Saury et al. (2014), but agrees well with more recent investigations by Wareing et al. (2019). The relation between spectral indices and phase transitions is one of our topics. We also consider correlations in position and velocity space, relating the more extended LNM as an envelope to the CNM.

The local H i needs to be distinguished from intermediate and high velocity clouds (IVCs and HVCs, respectively). These features are thought to be located in the Galactic halo, IVCs at kpc distances and probably originating from a Galactic fountain and HVCs at larger distances with an external origin (Wakker & van Woerden (1997) and van Woerden et al. (2004)). Both cloud types have a two-phase structure, but are distinct from Galactic H i because of their morphology and their distribution in center velocities and velocity widths (Haud, 2008). We consider the question of how much the IVC and HVC power spectra differ from the local distribution.

Several turbulence studies are available in the literature (Crovisier & Dickey (1983), Kalberla & Mebold (1983), Green (1993), Stanimirovic et al. (1999), Deshpande et al. (2000) ,Dickey et al. (2001), Elmegreen et al. (2001), Stanimirović & Lazarian (2001), Miville-Deschênes et al. (2003), Khalil et al. (2006), Miville-Deschênes & Martin (2007), Chepurnov et al. (2010b), Roy et al. (2010), Dedes & Kalberla (2012), Roy et al. (2012), Pingel et al. (2013), Martin et al. (2015), Kalberla & Kerp (2016), Kalberla et al. (2017), Blagrave et al. (2017), Pingel et al. (2018), and Choudhuri & Roy (2019)), yet a major problem is that the scatter of derived power indices is appreciable. The published range for 3D power law indices is $-2.2>\gamma>-4$.

We study how closely spectral indices depend on the composition of the H i in different phases, but we must also consider that the large scatter in $\gamma$ may be caused by the derivation of the results from small patches on the sky under very different physical conditions: absorption data from H i with considerable optical depth up to emission from diffuse high latitude gas, distant H i close to the Galactic plane up to local gas. Different distances imply variable linear scales perpendicular to but also along the line of sight. Last but not least, instrumental issues such as beam effects, instrumental noise, and apodization may cause uncertainties for the derived parameters. In some of the publications these instrumental issues are not discussed, others mention these effects and argue that beam smoothing and noise do not affect their analysis. Only a few publications consider instrumental biases in detail when fitting spectral indices. Kalberla & Kerp (2016, Fig. 22) demonstrated that a 3D spectral index of $\gamma=-2.97$ can steepen to $\gamma=-3.4$ if the beam correction is disregarded. Apodization and and instrumental noise can also lead to biases (Kalberla et al., 2017, Appendix A).

Biases caused by neglecting instrumental issues may even be larger than noted in the previous example. Blagrave et al. (2017) rectify a previously determined value of $\gamma=-3.6\pm 0.2$ (Miville-Deschênes et al., 2003) toward Ursa Major to $\gamma=-2.68\pm 0.14$ and explain biases in their Appendix E. The independent confirmation with $\gamma=-2.68\pm 0.07$ by Kalberla & Kerp (2016) demonstrates that a proper data reduction can deliver reliable results that are telescope independent. Particularly problematic is the slope of $-3.6$ that was considered in the review by Hennebelle & Falgarone (2012, Sect. 4.22, Fig. 10) as characteristic for high latitude H i emission. This incorrect value has in turn propagated to other publications and is still cited in the most recent papers without mentioning the revised value of $\gamma=-2.68$ for the Ursa Major region.

This paper is organized as follows. In Sect. 2 we briefly describe the data reduction. Details are given in Appendix A, and throughout the paper we discuss possible instrumental biases. In Sect. 3 we derive the spatial power distribution for different H i phases and correlations between WNM, LNM, and CNM power spectra. Narrowband (thin slices in velocity) spectral indices are derived in Sect. 4, dependences of spectral indices on the velocity channel width (thick slices) are discussed in Sect. 5. In Sect. 6 we discuss restrictions on the analysis of the turbulent flow in relation to the Kolmogorov (1941) theory for turbulence of locally homogeneous and isotropic structures. Spectral indices for the density and velocity correlation functions are discussed in Sect. 7. In Sect. 8 we derive power spectra for IVCs and HVCs. In Sect. 9 we discuss the power at low multipoles and the outer scale of turbulence. We summarize and discuss our results in Sect. 10.

2 Derivation of power spectra

We use high resolution 21 cm line data from the HI4PI survey (HI4PI Collaboration et al., 2016). To calculate the spatial power spectra we use ANAFAST provided by version 3.40 of the HEALPix software package222https://sourceforge.net/projects/healpix/ (Górski et al., 2005). For multipoles $l$ this routine calculates the power spectrum $P(l)$ of a HEALPix map333For calculations of the CMB, power spectra of the kind $l(l+1)C_{\mathrm{l}}$ are usually considered. Here we do not follow this convention, but use simply $C_{\mathrm{l}}$. This way the power exponent $\gamma$ increases by two.

[TABLE]

where $C_{\mathrm{l}}(l)$ are the correlation coefficients; we also define a power law index $\gamma$. Since our signal is bandwidth limited, we use in general $l<l_{\mathrm{max}}=1023$.

The parameter $P(l)$, however, is not the observed power spectrum. The observations are affected by instrumental issues, data processing, and even by selection effects from windowing particular regions on the sky. First, the intensity distribution on the sky is smoothed by the effective beam function $B_{\mathrm{data}}$ of the telescope. This includes beam smoothing as well as smoothing caused by the gridding process, both causing an artificial steepening of the observed power distribution. In addition, some noise power $N_{\mathrm{oise}}$, depending on instrument and observing method, is added by the receiving system

[TABLE]

The beam-corrected noise contribution $N_{\mathrm{oise}}(l)/B_{\mathrm{data}}^{2}(l)$ can be critical since beam effects suffer from a vanishing beam response $B_{\mathrm{data}}(l)$ at high multipoles.

The observed power distribution is further degraded by the window function, causing a convolution of the observable power distribution with the Fourier transform of the window function. These effects can be mitigated by the use of a proper apodization (Harris, 1978). Finally, fitting the power spectra includes statistical and systematical errors.

All these issues are discussed in detail in Appendix A. There we specify in particular the dependence of the noise term $N_{\mathrm{oise}}(l)$ on gridding and data processing and demonstrate that the noise term is unimportant for our application, using HI4PI data with an excellent signal-to-noise ratio. For $l<1023$ we obtain the clean and transparent relation

[TABLE]

3 Spatial power distribution for different H i phases

For $l\ga l_{\mathrm{crit}}\sim 8$ the beam corrected power law spectra $P(l)$ for the observed column density distributions are in most cases close to a straight line in log-log presentation and can be fit well with a constant spectral index (see, e.g., Fig. 32). The spectral index for $l\la l_{\mathrm{crit}}$ is rather flat and not well defined. The critical multipole $l_{\mathrm{crit}}$ depends on the outer scale of turbulence $L$ and the scale height $H$ of the turbulent medium, $l_{\mathrm{crit}}\sim 2\pi H/L$. We postpone the detailed discussion of the low multipoles to Sect. 9.

Turbulence in a two-phase H i medium was previously considered by Martin et al. (2015). In their study the CNM was defined by Gaussian components with Doppler temperatures $T_{\mathrm{D}}<443$ K. These authors found in the case of the north ecliptic pole that the spectral index of $-2.86\pm 0.04$ for the total H i changes to $-1.9\pm 0.1$ for the CNM and $-2.7\pm 0.1$ for the WNM. In addition to changes in the spectral index they found “an additional uncharacterized noise component in the $N_{{H\textsc{i}}}$ maps near the pixel scale, reflected at high spatial frequencies in the power spectrum.”

Our analysis is based on a three-phase all-sky Gaussian decomposition of the local H i gas (Kalberla & Haud, 2018). The different phases were shown to have different spatial distributions. The CNM is clumpy and embedded in the spatially more extended LNM which covers in addition a larger velocity spread around the CNM with narrow lines. Observed column densities are anti-correlated. The WNM with broad lines is embedding both LNM and CNM. These phase space relations should have correspondences in the associated power distributions.

3.1 Power spectra for CNM, LNM, and WNM column densities

The power spectra in Fig. 1 were calculated for the column density distributions shown by Kalberla & Haud (2018) in their Fig. 9 on top. Figure 1 shows the spatial power distributions for column densities from Gaussian components assigned to WNM (top left), LNM (top right), and CNM (bottom left) for a single channel at $v_{\mathrm{LSR}}=0$ km s*-1*.

Using the sum of all phases to calculate the power spectra for all Gaussian components results in the distributions shown in the bottom panel of Fig. 1 on the right-hand side. These are essentially the power spectra shown in Fig. 32 with minor deviations. Brightness temperatures restored from Gaussian components may deviate within the noise from observations, furthermore components that most probably are caused by radio frequency interference or instrumental problems have been rejected. Such a Gaussian based cleaning process has been proven a very efficient tool to eliminate residual instrumental problems (Kalberla & Haud, 2015). The power spectra of the resulting deviations from the observed H i distribution are plotted in the lower right panel of Fig. 1 to demonstrate that these deviations, including the noise floor, are far below the individual power spectra (see Appendix A). The shape of the power distribution from the expected scatter caused by the Gaussian decomposition (Fig. 1, bottom right) does not reflect the derived strong increase in the power for the individual phases. We conclude that the enhanced power at high multipoles $l\ga 100$ must be significant and caused by CNM structures on small scales. These structures were described by Clark et al. (2014) and Kalberla et al. (2016) as filamentary. They are unresolved by HI4PI and even by the narrow Arecibo beam, indicating a linear scale of $\la 0.1$ pc.

It is highly unexpected to find small-scale structures in the WNM power spectra at high multipoles; however, this is just a consequence of the correlations between the H i phases. Phase transitions generate CNM structures on small scales with some unstable LNM around the CNM (Kalberla & Haud, 2018). It is easy to understand that such structures are also present in the more extended LNM, but they are necessarily also reflected in the WNM power distribution with enhancements at high multipoles. The WNM distribution is not smooth at high spatial frequencies.

The power distribution $P_{\mathrm{TOT}}$ from all Gaussian components in the lower right panel of Fig. 1 can also be calculated by summing up all auto power spectra for the individual phases and all associated cross power spectra $P_{\mathrm{WxC}}$, $P_{\mathrm{WxL}}$, and $P_{\mathrm{LxC}}$ between the phases that are needed for a complete description of a three phase medium:

[TABLE]

The cross terms are shown in Figs. 2 and 3. In particular, multipoles $l\ga 100$ are of interest. The auto power enhancements for $P_{\mathrm{WNM}}$, $P_{\mathrm{LNM}}$, and $P_{\mathrm{CNM}}$ at high multipoles are caused by pronounced systematic anti-correlations in the cross power between WNM and LNM ($P_{\mathrm{WxL}}$), and LNM and CNM ($P_{\mathrm{LxC}}$). There is only a weak anti-correlation between WNM and CNM ($P_{\mathrm{WxC}}$), an effect that was noted by Kalberla & Haud (2018, Sect. 4.2). We note that in Figs. 2 and 3, unlike all the other power law plots, we do not use a logarithmic scale. The negative cross power on low multipoles is significant.

The power distributions shown in Fig. 1 for column densities in individual phases deviate significantly from the power spectra of the total H i, combining all phases. All single-phase slopes, derived at $10<l<100$, are shallow, with $\gamma\sim-2.5$ in comparison to $\gamma\sim-3$ for $P_{\mathrm{TOT}}$. In all cases we find enhanced power for $l\ga 100$, strongly increasing to high multipoles. These issues can only be explained by the fact that the column density distributions in different phases are highly correlated. The cross terms in Eq. 4 are important and describe the interplay between different phases.

It is important to make this point clear. As an example we consider a layered structure, consisting of several completely independent and uncorrelated sheets of different phases along the line of sight. In this case the cross terms in Eq. 4 would vanish, resulting in $P_{\mathrm{TOT}}=P_{\mathrm{WNM}}+P_{\mathrm{LNM}}+P_{\mathrm{CNM}}$. This, obviously, is not observed. Another example for vanishing cross terms is the case of disjunct spatial distributions. Separate power spectra from different hemispheres, as shown in Fig 30, may be added since the cross terms for disjunct regions are zero. The established assumption, in the case of a two-phase medium, that “the density fluctuations are spatially separated in two media and therefore their correlation is likely to be negligible” (Lazarian & Pogosyan, 2000, Sect. 4.3) is not justified444In Eq. 51 of Lazarian & Pogosyan (2000) the cross term between CNM and WNM in the case of a two-component medium is defined without the factor two in our Eq. 4.. The H i with its different phases needs to be considered a composite. It is not possible to describe the H i with several independent phases.

3.2 Optical depth effects

The H i survey data we analyzed may suffer from optical depth effects. For the CNM the optical depth is expected to increase with decreasing spin temperature, though a general physically significant relationship cannot be established (Heiles & Troland, 2003). An increase in the optical depth for the clumpy CNM can lead to a systematical underestimation of the observed column densities and this implies that the power at high multipoles is underestimated. In turn, observed power spectra may be steepened artificially.

We use an empirical correction derived by Lee et al. (2015) from Arecibo observations in direction to the Perseus molecular cloud which is consistent with data from Heiles & Troland (2003) and which was confirmed later by Murray et al. (2018a) and Murray et al. (2018b). Observed column densities $N_{\mathrm{Hobs}}$ are accordingly biased by a factor

[TABLE]

and we apply a correction $N_{\mathrm{Hcorr}}=f\cdot N_{\mathrm{Hobs}}$ for $f>1$, this concerns CNM with $N_{\mathrm{Hobs}}\ga 0.6\leavevmode\nobreak\ 10^{20}\ {\rm cm^{-2}}$. Using this correction we recalculate the power spectra for $v_{\mathrm{LSR}}=0$ km s*-1*, as shown in Fig. 1. We find in all cases that power spectra for the CNM and the total H i get flatter. However, the resulting changes of the spectral indices are typically less than half of the rather low formal uncertainties of the fit. The spectral indices are within the errors unaffected by optical depth effects; only in one case do we find a slight change. In the case of the total all-sky H i column densities the index changes from $\gamma=-2.776\pm 0.004$ to $\gamma=-2.768\pm 0.003$. We conclude that optical depth effects cannot explain the observed steepening of the single-channel power spectra close to zero velocities. Optical depth effects are noticeable at high column densities in the Galactic plane. High latitude power spectra are essentially unaffected and we draw our conclusions predominantly from these data.

3.3 Power spectra for CNM, LNM, and WNM phase fractions

The power spectra discussed in Sect. 3.1 are caused by column density fluctuations of individual H i phases. Hence there are two effects competing with each other, fluctuations in column density and in phase composition. We want to determine power spectra merely related to individual H i phases, and therefore consider average phase fractions along the line of sight.

Phase fractions $f_{\mathrm{P}}(_{v1}^{v2})$ depend in general on the velocity range and, following Kalberla & Haud (2018), for $v1<v_{\mathrm{LSR}}<v2$ are defined as

[TABLE]

where $T_{\mathrm{b}}$ is the observed brightness temperature, while $T_{\mathrm{bP}}$ stands for the brightness temperature contribution from phase P, P = CNM, LNM, or WNM. Since the dependences on column densities cancel in this expression, the weights for the decomposition into different phases are identical for all lines of sight. Derived multiphase power spectra are essentially deconvolved for column density effects. For a visualization of the spatial column density distributions for different phases in comparison to associated phase fractions, we refer to Kalberla & Haud (2018, Figs. 9 and 10).

In Fig. 4 we show power spectra from single-channel phase fractions at the velocity $v_{\mathrm{LSR}}=0$ km s*-1*. Comparing the power spectra from Fig. 4 with those from Fig. 1 shows that now the spectra for all-sky and high latitudes are in far better agreement than before, but spectra for |b|>20\hbox{{}^{\circ}} are slightly steeper than the all-sky spectra. We observe significantly increased power at multipoles $l\la l_{\mathrm{crit}}$. In the case of phase fractions the cosecant latitude effect of the Galactic disk is removed and possible biases (see Sect. 9) for low multipoles are minimized. The spectral indices for the phase fractions are $-1.8>\gamma>-2.2$ and are significantly flatter than $-2.14>\gamma>-2.48$ for the column density distributions discussed in Sect. 3.1. This flattening is most pronounced for the LNM and can be explained by the absence of column density fluctuations.

Phase fractions need to sum up to one, the bottom panel on the right side of Fig. 4 shows in the all-sky case the completely negligible power spectrum for deviations caused by the Gaussian decomposition (black). We note that this power spectrum now reflects the increased uncertainties at low multipoles from baseline uncertainties at the field boundaries discussed in Sect. A.2.3. The red curve shows apodization errors at high latitudes. They cause a considerable ringing, but fortunately such unavoidable errors are also not significant since they drop off very fast at high multipoles.

Phase fractions, defined by Eq. 6 are useful to study the gas composition independent from its total column density. These data do not represent any reasonable model for the multiphase H i distribution, but are useful to consider how strictly power spectra depend on the phase composition. Phase transitions cause shallow power spectra for $l\la 100$ and enhanced power at higher multipoles.

4 Velocity dependent narrowband spectral indices

We consider now velocity dependences in the spectral indices. Figure 5 shows the average distribution of phase fractions for the local ISM at high latitudes in the velocity range $-16<v_{\mathrm{LSR}}<16$ km s*-1* (Kalberla & Haud, 2018, Sect. 3). There are considerable fluctuations in the phase fractions and an imprint on the power distribution for individual phases at these velocities may be possible.

We consider first the model distributions of phase fractions that are unaffected by column density fluctuations. The velocity dependences of narrowband spectral indices for multipoles $10<l<100$ are shown in Fig. 6 for all-sky and in Fig. 7 for |b|>20\hbox{{}^{\circ}}. We observe a general steepening of the spectral indices close to $v_{\mathrm{LSR}}\sim 1$ km s*-1*. The strongest effect is found for the CNM. The enhancement of the average CNM phase fraction at this velocity (Fig. 5) clearly leads to a pronounced steepening of the CNM power spectra. The steepening of the spectral indices of the other phases is noteworthy, but may be a consequence of the coupling of the power spectra between different phases according to Eq. 4.

Next we consider the power spectra derived from the Gaussian components of different H i phases. Figures 8 and 9 show still local minima for the CNM spectral indices, however the spectral indices for the other phases are now mostly anti-correlated; the WNM has a particularly smooth distribution. Spectral indices derived from power spectra of the total observed H i also show smooth changes. For high latitudes (Fig. 9, red) we find a steepening in anti-correlation to the WNM, but with local minima that are affected by the CNM and LNM. In the all-sky case (Fig. 8, black) we observe a smooth trend, but influences from the CNM and LNM are weak.

Velocity dependences for the spectral indices have been observed before (Kalberla & Kerp (2016) and Kalberla et al. (2017)). Recently, Choudhuri & Roy (2019) have presented data that show a steepening of the spectral indices at velocities that are considered to be representative for the CNM. Spectral indices for 13CO 2–1 and 12CO 3–2 channel maps in the Perseus cloud were observed by Sun et al. (2006) with a steepening of the indices that can be described as an average structure of the index spectrum similar to the line profile.

Numerical studies of the condensation of WNM into CNM structures caused by turbulence and thermal instabilities were conducted by Saury et al. (2014). They found that turbulence plays a key role in the structure of the cold medium; it can induce the formation of CNM when the WNM is pressurized and put it in a thermally unstable state. They found evidence for subsonic turbulence with a shallow power index $\gamma\sim-2.4$. Recent high resolution simulations by Wareing et al. (2019) have shown that the thermal instability dynamically can form sheets and filaments on typical scales of 0.1 to 0.3 pc, explaining that as the structure grows in the simulation, the density power spectrum rapidly rises and steepens. Wareing et al. (2019, Fig. 10) derive $\gamma\sim-3,$ but the steepening is most prominent on scales below 1 pc.

The preferred width of the filaments in the range 0.1 to 0.3 pc is consistent with unresolved cold small-scale structures observed by Clark et al. (2014) and Kalberla et al. (2016). Furthermore, these cold filaments have characteristic radial velocities close to zero km s*-1*. The single-phase power spectra show a strong increase at the corresponding multipoles $l\sim 1023$.

A steepening of thin velocity slice H i power spectra in a narrow velocity range, associated with a decrease in the WNM fraction and the coexistence of cold anisotropic CNM filamentary structures was reported before for three fields at intermediate latitudes by Kalberla et al. (2017, Sect. 6). The interpretation was that phase transitions, indicated by changes in H i phase fractions, modify the narrowband power distribution. For a determination of dependences between power spectra and phase transitions these authors proposed isolating particular cold H i components from the rest of the H i distribution in order to differentiate between ON velocities at the line centers and appropriate OFF velocities in the line wings. Accordingly we define here

[TABLE]

as the power ratio and

[TABLE]

as the power excess, relating single-channel ON and OFF power distributions; $P_{\rm OFF}(l)$ is defined as the average single-channel power at two OFF positions.

The choice of the ON velocity is obvious, we observe at $v_{\mathrm{LSR}}=3$ km s*-1* the highest CNM and at the same time the lowest WNM ratio (Fig. 5). The high latitude data also show at this velocity a well-defined minimum for the CNM spectral index (Fig. 9). From this plot the definition for the OFF velocities at $v_{\mathrm{LSR}}=3\pm 5$ km s*-1* is also well defined, the CNM power law steepening is limited to the range $\Delta v_{\mathrm{LSR}}\sim 10$ km s*-1*; we refer to Sects. 5.3 to 5.6 for further discussion of the velocity spread.

Figure 10 shows the power ratios $\mathfrak{R}(l)$. We observe $\mathfrak{R}(l)\ga 1$ for all multipoles except $l\sim 30$. For the WNM $\mathfrak{R}(l)$ is flat but LNM and CNM show enhancements at high multipoles. Strong fluctuations are also observable at multipoles $10\la l\la 50$.

Figure 11 shows the power excess $\mathfrak{E}(l)$. In comparison to Fig. 1 ($\gamma=-2.776\pm 0.004$ for the total H i and $\gamma=-2.37\pm 0.03$ for the CNM) we observe steep power spectra with $\gamma=-2.87\pm 0.02$ for the total H i for $l\ga 10$, and in particular $\gamma=-2.8\pm 0.1$ for the CNM for $10\la l\la 100$. This is by far the steepest power law derived by us for the CNM, and we interpret the steep excess power distribution $\mathfrak{E}(l)$ for the CNM as induced by phase transitions. Thermal instabilities obviously also affect the power spectra for WNM and LNM; they are no longer straight, but bent up systematically with increasing power for high multipoles. The different phases are correlated according to Eq. 4, but the power spectrum for the total H i remains straight.

In agreement with Kalberla et al. (2017) we find strong evidence that the power distribution in the local ISM is significantly steepened by phase transitions. The CNM is filamentary and anisotropic on small scales (Kalberla & Kerp (2016) and Kalberla et al. (2017)), but even the large CNM power ratios $\mathfrak{R}(l)$ at low multipoles appear to be linked to large-scale filaments, oriented along the magnetic field as observed by Clark et al. (2014) and Kalberla et al. (2016).

5 Spectral indices and velocity channel widths

In this section we determine spectral indices for individual H i phases after integrating the column density distributions over several channels with a total full width at half maximum (FWHM) of $\Delta v_{\mathrm{LSR}}$. Figures 12 and 13 show the results. From the all-sky data we note that the spectral indices for WNM and LNM steepen continuously. The WNM has broad lines, resulting in a slow change in spectral index. The characteristic LNM line widths are narrower, corresponding to faster changes of the slope with $\Delta v_{\mathrm{LSR}}$. For the CNM the situation is more complex. In the case of the all-sky data we note a prominent minimum at $\Delta v_{\mathrm{LSR}}=5$ km s*-1* and a rise afterward. At high latitudes the situation is similar, but the CNM spectral index now has a broad minimum around $\Delta v_{\mathrm{LSR}}\sim 10$ km s*-1*.

From the observed H i column densities or the sum of all phases we derive the spectral indices shown in Fig. 14. Spectral indices decrease in a similar way for the all-sky data and for |b|>20\hbox{{}^{\circ}} until $\Delta v_{\mathrm{LSR}}=16$ km s*-1*. Afterward the indices diverge; they fall for all-sky data but rise for |b|>20\hbox{{}^{\circ}}. This result is unexpected, and an explanation is needed.

5.1 Velocity channel analysis

Power spectra extracted from narrowband data (thin velocity slices) are affected by turbulence in two ways. Density fluctuations have an imprint on the observed column densities, but velocity fluctuations also need to be taken into account. To separate effects caused by fluctuations in density or velocity, the basic idea of VCA is to increase the width $\Delta v_{\mathrm{LSR}}$ of the velocity slice. Velocity fluctuations should average out, and when the velocity window $\Delta v_{\mathrm{LSR}}$ is broad enough that all internal velocities are covered (thick velocity slice) the observed emissivity should eventually be dominated by density fluctuations (Lazarian & Pogosyan, 2000, Eq 25). The expected minimum velocity width is $\Delta v_{\mathrm{LSR}}\sim 17$ km s*-1* (Lazarian & Pogosyan, 2000, Sect. 4.3 ), the typical FWHM thermal line width for the WNM. Starting with a thin velocity slice and successively increasing the velocity width, a gradual steepening of the spectral index should be observed until convergence. For our application to H i column densities we observe a steepening (Fig. 14), but in the case of high latitude data only until $\Delta v_{\mathrm{LSR}}\sim 16$ km s*-1*, which is the maximum linewidth for consistent results. Considering the CNM in Figs. 12 and 13 we find a steepening, but subsequently the spectral index flattens again. The question arises whether VCA is applicable to our data. We take up the discussion later in Sect. 7.

Spectral indices close to $\gamma\sim-2.95$ and within the uncertainties of $\Delta\gamma\la 0.1$ independent of velocity widths $0.82<\Delta v_{\mathrm{LSR}}<21.4$ km s*-1* have been reported by Khalil et al. (2006, Tables 3 and 4, and Fig. 29). A constant spectral index is inconsistent with VCA; these authors discuss in their Sect. 6 some other observations that do not show the steepening predicted by VCA. Yuen et al. (2019) note that the analysis of the data agrees well with VCA predictions and revise the Khalil et al. (2006) indices in their Table 3 to $\gamma_{\mathrm{thin}}=-2.6$ and $\gamma_{\mathrm{thick}}=-3.4$. No details are given how this result was obtained, but the authors refer to a forthcoming more detailed paper.

5.2 Velocity channel cross-correlations

Trying to better understand the dependences between turbulent power distribution and the velocity window width $\Delta v_{\mathrm{LSR}}$, we calculate cross-correlations for single-channels in the wings of the line signal. Hence we cross-correlate channels separated by a distance $\Delta v_{\mathrm{LSR}}$ in velocity. This is a differential approach. The cross power spectra that we use to study internal velocity dependences of individual H i phases are, according to the Wiener-Khintchine theorem, the Fourier transforms of the cross-correlation function at particular characteristic lags in velocity. In the case of a signal from a genuine turbulent medium, the cross power spectra should inform us on how well correlated the individual H i phases remain throughout their velocity distribution.

We consider high latitudes first; all-sky data is deferred until Sect. 5.8. Figure 15 shows us the cross-correlations between two channels at velocities $v_{\mathrm{LSR}}=-5$ km s*-1* and $v_{\mathrm{LSR}}=+5$ km s*-1*. For the WNM we find a well-defined power spectrum with a slope of $\gamma=2.41\pm 0.06$ for multipoles $10<l<100$ with enhanced power at high multipoles. For the LNM no power law can be found for $10<l<100$ (only positive values of $C_{l}$ can be plotted), but for this phase we also obtain enhanced power at high multipoles; the velocity cross power distributions for LNM and WNM are almost identical there. In the case of the CNM only a noisy cross-correlation signal at high multipoles can be seen. In comparison to the expected contribution of the uncertainties from the Gaussian decomposition shown in Fig. 1 this signal ($\sim 10^{30}$ cm*-4*) is barely significant. The local minimum of the spectral index for the CNM phase in Fig. 15 at $\Delta v_{\mathrm{LSR}}=10$ km s*-1* indicates that the correlation of the CNM is lost. This implies that on average any CNM structure seen in channels separated by $\Delta v_{\mathrm{LSR}}\ga 10$ km s*-1* is unrelated. Homogeneous CNM clouds have an average FWHM line width of 3.6 km s*-1* (Kalberla & Haud, 2018) and according to Sect. 5.5 only 3.2 km s*-1* for velocities $|v_{\mathrm{LSR}}|\la 5$ km s*-1*. The internal velocity distribution of individual CNM clouds is too narrow to cause a limitation of this kind. The implication must be that CNM clumps in a turbulent flow decouple from each other if their center velocities differ by $\sqrt{10^{2}-3.6^{2}}\ga 9$ km s*-1*. In the following we investigate this remarkable result further; in Sect. 6 an explanation is offered from constraints considered by Kolmogorov (1941). Both the WNM and LNM maintain their correlation at a lag of 10 km s*-1*. Their cross power spectra are identical for $l\ga 200$; in Fig. 13 the LNM data overplot the WNM data. Increasing the velocity separation for two channels further, we also observe for the LNM a decorrelation of the cross power for $\Delta v_{\mathrm{LSR}}\ga 26$ km s*-1*. This is discussed later in Sect. 5.4 and demonstrated in Fig. 19.

5.3 Constraints on spectral coherence

A different way to explore correlations in the velocity domain is by defining the spectral coherence $S_{l}(v_{1},v_{2})$ between two single channels at velocities $v_{1}$ and $v_{2}$

[TABLE]

where $C_{l}(v_{1},v_{2})$ is the cross-correlation between individual PPV channel maps at velocities $v_{1}$ and $v_{2}$ with corresponding autocorrelations $C_{l}(v_{1},v_{1})$ and $C_{l}(v_{2},v_{2})$. This relation has the advantage that it is independent of beam smoothing as long as the H i emission is little affected by noise. It is clear that $S_{l}=1$ for $v_{1}=v_{2}$ and the signal decorrelates for increasing channel separation, but we require for a turbulent flow that some correlation should be maintained if we want to define eddies as coherent patterns of flow velocity, vorticity, and pressure.

We use high latitude data to calculate $S_{l}(v_{1},v_{2})$ for the WNM, LNM, CNM, and the total H i column density distributions at velocities $v_{1}=-5$ km s*-1* and $v_{2}=+5$ km s*-1*. The results are shown in Fig. 16. The WNM distribution is coherent at all multipoles $l$. This is expected since the velocity spread $v_{2}-v_{1}=10$ km s*-1* is small compared to the average FWHM width of 23.3 km s*-1* of the WNM emission (Kalberla & Haud, 2018). In the case of the LNM with an average FWHM 9.6 km s*-1* some spectral coherence remains at high multipoles. This behavior is expected for an unstable phase that is embedded in the WNM, but surrounds a clumpy CNM distribution (see Sect. 5.7). For the CNM with an average FWHM width of 3.6 km s*-1* the fit $\gamma=-1.07\pm 0.09$ at a channel separation of $\Delta v_{\mathrm{LSR}}=10$ km s*-1* is within the uncertainties consistent with random Gaussian noise ($\gamma=-1$). For comparison, the spectral coherence is also shown for the total multiphase H i column density (Fig. 16, red). This is expected close to one, but interestingly the scatter in $S_{l}$ is large compared to the WNM. Clearly the WNM derived from Gaussian components provides a more consistent presentation of the large-scale H i distribution.

We repeat the calculations and determine $S_{l}$ for $v_{2}-v_{1}<10$ km s*-1*; coherence for the CNM is recovered. As an example we use $S_{l}(v_{1},v_{2})$ for $v_{1}=-2$ km s*-1* and $v_{2}=+2$ km s*-1*. This velocity lag corresponds roughly to the average FWHM width of 3.6 km s*-1* for the CNM. The result is shown in Fig. 17. For the WNM the situation is similar to Fig. 16, but for the CNM we find coherence for $l>100$. This situation can be described as follows: power in the CNM exists spatially predominantly on small scales, and structures in velocity space can only be detected for small velocity lags. Increasing the velocity lag leads to a gradual decorrelation; as shown in Fig. 16, the spectral coherence is lost completely for $\Delta v_{\mathrm{LSR}}=10$ km s*-1*; the CNM in the two separate channels decouples if this value is exceeded.

For the LNM we find correlation on larger scales (up to smaller multipoles) and at the same time for larger velocity lags. The LNM surrounds the CNM, occupying a larger volume and at the same time a larger velocity spread. Some evidence for this hierarchical structure was given before by Kalberla & Haud (2018). Increasing the lag $v_{2}-v_{1}$ we find a gradual decorrelation for the LNM, the LNM decouples for $\Delta v_{\mathrm{LSR}}\ga 23$ km s*-1*.

5.4 Self-similarities in velocity space

The missing correlation between CNM clouds with velocity lags $\delta v_{\mathrm{LSR}}\ga 10$ km s*-1*, and the good spectral coherence of the CNM at lags $\delta v_{\mathrm{LSR}}\sim 4$ km s*-1* shown in Fig. 17 do not imply that the CNM distribution is restricted in velocity space. The CNM can exist independently at other bulk velocities $(v_{1}+v_{2})/2$ and other domains embedded in LNM.

To demonstrate this we repeat the calculations and determine $S_{l}(v_{1},v_{2})$ for $v_{1}=6$ km s*-1* and $v_{2}=10$ km s*-1*. The result is shown in Fig. 18 and resembles the CNM in Fig. 17. The local H i is self-similar concerning the properties of the phase distribution in velocity space even though the column densities are significantly different at $v_{\mathrm{LSR}}\sim 0$ and $v_{\mathrm{LSR}}\sim 8$ km s*-1*. Differences in phase fractions (Fig. 5) and spectral indices (Fig. 7 and 9) also exist, but they have no significant imprint on the spectral coherence of the CNM. For the LNM we find in comparison to Fig. 17 an increase in spectral coherence at multipoles $l\la 100$, caused by the reduced contribution of CNM to the H i at velocities $\delta v_{\mathrm{LSR}}\sim 8$ km s*-1*.

Self-similarities are also recognizable for the LNM though the velocity limits are in this case less sharply defined than in the case of the CNM. Figure 19 shows examples of lags between velocities $v=$ -15, 0, and 15 km s*-1*. The spectral coherence in the LNM is lost at large velocity lags, but can be recovered by gradually decreasing the width of the velocity lag, independent of the bulk velocity $(v_{2}+v_{1})/2$. The example in Fig. 19 also shows that the spectral coherence values for the LNM at negative and positive bulk velocities are comparable, while slight differences exist for the WNM.

5.5 Spectral coherence at high multipoles

The dependence of $S_{l}(v_{1},v_{2})$ (Eq. 9) on the lag $\Delta v=v_{2}-v_{1}$ can be best studied at high multipoles. We consider $v_{\mathrm{LSR}}=(v_{2}+v_{1})/2=0$ km s*-1* and use high latitude data at |b|>20\hbox{{}^{\circ}}. To suppress random fluctuations we average $S_{l}(v_{1},v_{2})$ for $l\ga 1000$. The results are shown in Fig. 20 with symbols. We try to model these $S_{l}$ distributions by assuming that they are caused by a sample of H i clumps with characteristic Doppler temperatures and corresponding Gaussian line shapes. This works well for the CNM, but fails in the case of the the LNM. For the CNM we fit a FWHM width of $3.20\pm 0.01$ km s*-1*; this Gaussian is shown in Fig. 20. For comparison we plot for the LNM a Gaussian with FWHM width of 10 km s*-1*. It is obvious that the LNM cannot be approximated by a Gaussian distribution with a single characteristic line width. This implies that the CNM and the LNM must have very different properties. The CNM can be described as a phase that occupies a well-defined range in velocity width. The LNM as an unstable phase may be “best understood as a range of density and temperature values” (Vázquez-Semadeni, 2012). A large scatter in Doppler temperatures (or line widths) observed by Kalberla & Haud (2018, Fig. 7) supports this proposal. The missing stability of the LNM implies that H i gas in this phase has to develop either to CNM or to WNM; the line width distribution gets stretched out in both directions.

The fitted CNM line width of $3.20\pm 0.01$ km s*-1* corresponds to a Doppler temperature of $T_{\mathrm{D}}=223.7\pm 2$ K. This temperature was reported previously by Kalberla et al. (2016) for cold filamentary H i gas. The agreement is at first glance surprising since filamentary structures, derived from unsharp masking (USM), have typically low column densities. Kalberla & Haud (2018) determined an average CNM line width of $3.6\pm 0.1$ km s*-1*, corresponding to a Doppler temperature of 283 K, from profiles with less than eight well-defined Gaussian components. They assumed that the Gaussian analysis would recover only a part of the weak components. The cold filamentary USM structures, however, turn out to be characteristic of the small-scale structure of the CNM if velocities are restricted to the velocity range $-5<v_{\mathrm{LSR}}<5$ km s*-1*.

Unsharp masking aims to extract H i structures at the highest spatial frequencies that can be observed by large single-dish radio telescopes. The Doppler temperature was determined by Kalberla et al. (2016) as the median and also by fitting a log-normal probability density distribution555For a turbulent medium the characteristic PDFs are log-normal (Vazquez-Semadeni, 1994) as a result of the central limit theorem applied to self-similar random multiplicative perturbations. The geometric mean of a log-normal distribution is equal to its median. . Here we select coherent structures at high multipoles and derive the column density weighted geometrical mean Doppler temperature from the correlation functions. The methods for both investigations are quite different, in each case with uncertain assumptions, but the resulting Doppler temperatures agree very well.

A characteristic temperature for the CNM can also be derived directly from the sample of Gaussian components. From narrowband data at high Galactic latitudes with $\delta v_{\mathrm{LSR}}=1$ km s*-1* we derive the lowest median $T_{\mathrm{D}}=208$ K at $v_{\mathrm{LSR}}=0$ km s*-1*. For the velocity range $-5<v_{\mathrm{LSR}}<5$ km s*-1* we obtain a median $T_{\mathrm{D}}=217$ K and a mean $T_{\mathrm{D}}=229$ K. These $T_{\mathrm{D}}$ values are deconvolved for the instrumental line broadening caused by the spectrometers and agree well with the value $T_{\mathrm{D}}=223.7\pm 2$ K, derived by fitting $S_{l}$. No deconvolution is necessary for $S_{l}$ determined from Eq. 9 because bandwidth and beam effects cancel in this case. Thus, the median Doppler temperature $T_{\mathrm{D}}=223$ K for small-scale CNM structures is well defined and confirmed by our current investigations.

Extracting prominent cold filaments (Kalberla et al., 2016) is equivalent to selecting H i components with steep spectral indices (see also Sect. 4). A completely independent determination of H i small-scale structures with the Arecibo telescope by Clark et al. (2014) leads, within the errors, to the same result; these authors found a Doppler temperature of 220 K. Clark et al. (2019) have proven recently that anisotropic magnetized H i small-scale structures and narrow linewidths are dust-bearing density structures. They found that anisotropic small-scale H i channel map structures are correlated in position and position angle with far infrared structures at 857 GHz, implying that this emission originates from a colder, denser phase of the ISM than the surrounding material.

Coherent CNM structures, described by $S_{l}$ according to Eq. 9, can be interpreted as eddies (see Clark et al. 2014, Figs. 3 and 4, or Kalberla & Kerp 2016, Figs. 28 and 29). These structures are caused by a well-defined population of CNM clouds with characteristic log-normal distributions in column density and Doppler temperature (Kalberla & Kerp, 2016, Figs. 12 and 13). The properties of these CNM structures are well defined. They are cold, dusty, magnetized, and aligned with the local magnetic field (Clark et al. (2014) and Clark et al. (2019)). Typical temperatures of the associated dust are $T_{\mathrm{d}ust}\sim 18.5$ K. These structures are embedded in LNM with $T_{\mathrm{d}ust}\sim 19$ K (Kalberla & Haud, 2018). Anisotropies associated with radio-polarimetric filaments, explored by Kalberla & Kerp (2016) and Kalberla et al. (2017), suggest that these structures are shaped by MHD turbulence in the presence of a magnetic field (Goldreich & Sridhar, 1995). Filamentary features from these observations are probably mostly sheets with systematic velocity gradients perpendicular to the field direction. According to Heiles & Crutcher (2005), “edge-on sheets should be edge-on shocks in which the field is parallel to the sheet.” Not all of the data can be explained this way. An alternative interpretation of filaments as fibers, having an approximately cylindrical geometry, was given by Clark et al. (2014).

A characteristic Doppler temperature of 223 K implies for an average spin temperature of 50 K (Heiles & Troland, 2003) a mean Mach number of $M\sim 3.8$ for the CNM, a value that agrees well with the estimate of $M\sim 3.7$ by Heiles & Troland (2005). Using a median spin temperature of 61 to 79 K (Murray et al., 2018b, Table 4) would result in $2.8\la M\la 3.3$.

According to Lazarian & Pogosyan (2000) small-scale structures in the H i do not reflect masses of real clumps, but are caustics produced by the turbulent velocity field via projection; this process is also called velocity mapping (Lazarian & Pogosyan, 2000, Sect. 6). Observational parameterized structures with properties as summarized in the last two paragraphs were newly rejected as being real entities. Lazarian & Yuen (2018) conclude that “the filaments observed by Clark et al. (2014) in thin channel maps can be identified with caustics caused by velocity crowding.” Clark et al. (2019) object and conjecture instead that “small-scale structures in narrow H i channel maps are preferentially cold neutral medium that is anisotropically distributed and aligned with the local magnetic field.” We also interpret these structures as density enhancements, CNM eddies from a turbulent flow, caused by phase transitions and observable at the resolution limit of large single-dish telescopes. Phase transitions may be partly provoked by turbulence (Saury et al., 2014), but cause at the same time back-reactions on the power distribution (see also Sects. 4 and 7.4).

5.6 Linewidth regimes for different phases

The FWHM velocity window width for the steepest CNM power index (Fig. 13) is identical with the lag where spectral coherence for the CNM is lost. This velocity window is also identical with the mean LNM line width of 9.6 km s*-1* determined by Kalberla & Haud (2018). The model assumption of the LNM as a transition between the stable clumpy CNM and diffuse WNM implies that the velocity space covered by the LNM must be at least as large as the velocity spread where we find coherence in the CNM. The LNM is unstable, and there is a high probability of finding it associated with CNM, which causes a coupling between both phases in velocity space.

By generalizing such a hierarchical scheme in the velocity space, we assume that the spectral power for the LNM distribution, which is also embedded in the WNM, should be correlated with the mean FWHM velocity width of 23.3 km s*-1* for the WNM (Kalberla & Haud, 2018). Figure 13 shows that this is the velocity width of the steepest LNM spectral index. For broader velocity windows the LNM spectral index is flattening since the correlation between channels separated by $\Delta v_{\mathrm{LSR}}\ga 23$ km s*-1* is lost.

We conclude that the mean FWHM velocity widths of the different phases are not independent from each other; instead, they reflect correlations between the H i phases. Saury et al. (2014) report from their simulations “that the turbulent motions inside clumps and the relative velocities between them are related to the motions of the WNM from which they were formed,” which means that “a non-negligible fraction of the measured velocity dispersion is caused by the relative motions of the clumps along the line of sight, suggesting that the observed line broadening is likely to be due to the relative clumps motions rather than supersonic turbulence.”

5.7 Spatial phase coherences

To study the spatial relationship between different phases we define, in analogy to Eq. 9, the spatial phase coherence

[TABLE]

between phases A and B. Here $C_{l}(\mathrm{AxB)}$ is the cross-correlation between PPV channel maps for phases A and B, and $C_{l}(\mathrm{AxA)}$ denotes the auto correlation for phase A, and $C_{l}(\mathrm{BxB)}$ for phase B. A and B stand for CNM, LNM, and WNM; $P_{l}$ depends in general not only on the H i phases, but also on the velocity range $v_{1},v_{2}$.

Figure 21 shows $P_{l}$ for single channels at $v_{\mathrm{LSR}}=0$ km s*-1*. Here $P_{l}(\mathrm{CNM,LNM})$ indicates that the CNM is dominating high multipoles $l\ga 100$; $P_{l}(\mathrm{LNM,WNM})$ shows that the LNM is somewhat more extended with $l\ga 50$; $P_{l}(\mathrm{CNM,WNM})$ indicates that CNM and WNM are spatially almost uncorrelated since the LNM is enveloping the CNM, see also Fig. 3. Within the noise a very weak correlation $P_{l}(\mathrm{CNM,WNM})$ remains; this signal is noisy because of the low volume filling factor of the CNM. To find a relation between CNM and WNM it is necessary to constrain the statistics by selecting positions with a significant phase fraction for the CNM (see Kalberla & Haud 2018).

The spatial phase coherence $P_{l}$ is self-similar for changes in bulk velocity $(v_{1}+v_{2})/2$ and velocity spread $v_{2}-v_{1}$ as long as the spectral coherence for the individual H i phases is not violated.

5.8 Galactic plane data – differential Galactic rotation

Next we consider the all-sky data shown in Fig. 14. Due to the strong emission in the Galactic plane, the all-sky cross-correlations for channels with a separation of $\Delta v_{\mathrm{LSR}}=26$ km s*-1* (see Fig. 22) are significant for WNM, LNM, and the sum of all phases. However, what we observe is not a signal that we expect for a genuine turbulent flow. The observed H i column density distribution in the Galactic plane is affected by differential Galactic rotation, causing at Galactic longitudes $Glon$ a $\mathrm{sin}(2\cdot Glon)$ modulation (e.g., Mebold, 1972). In particular, velocities close to zero are affected by velocity crowding (Burton, 1972), leading to a systematic degradation of the signal for $l\la 100$; the assumption of homogeneity and isotropy is no longer valid. High column densities from the Galactic plane dominate the statistics. Checking cross-correlation spectra with smaller velocity separation we can trace this problem to channel separations as low as $\Delta v_{\mathrm{LSR}}=16$ km s*-1*. HI4PI all-sky data, in short, are affected by confusion from the Galactic plane, and are therefore only of limited use to us.

6 Kolmogorov’s local constraint

Kolmogorov (1941) constrained his paper on “turbulence in incompressible viscous fluids for very large Reynolds numbers” to locally homogeneous and isotropic structures. The term “locally” was (except for stationarity in time) specified by Kolmogorov as “restrictions are imposed only on the distribution laws of differences of velocities and not on the velocities themselves.”

A remarkable result of our investigations in Sects. 4 and 5 is that the properties of the turbulent flow in the multiphase ISM are limited concerning velocity differences (lags $(v_{2}-v_{1})$ used by us). For the two phases, CNM and LNM, we find a decorrelation of the flow in velocity (Fig. 15) and describe this phenomenon as a loss in spectral coherence (Fig. 16). The decorrelation depends on spatial properties of the flow. The CNM dominates high multipoles $l\ga 100$ (Fig. 17), and decorrelation happens at lags $v_{2}-v_{1}\sim 10$ km s*-1*. The LNM exists at multipoles $l\ga 50$ (Fig. 21) and decorrelates at lags $v_{2}-v_{1}\sim 23$ km s*-1*. Thus, for both phases we have well-defined domains, and in agreement with Kolmogorov (1941) “the hypothesis of local isotropy is realized with good approximation in sufficiently small domains G of the four-dimensional space.” We find self-similarities, i.e., the properties of the CNM and LNM do not change significantly when shifting the velocities (Figs. 18 and 19), but decorrelation at large velocity differences remains, independent of the “velocities themselves” that are investigated. Turbulence in the local ISM is locally constrained but otherwise is homogeneous and isotropic (for a comparison between the power spectra from two hemispheres, see the lower right panel of Fig. 30).

Restricting turbulence to be locally constrained, Kolmogorov (1941) probably also had limitations from the experimental setup in mind; the domain under investigation should not be “lying near the boundary of the flow or its other singularities.” Except for the layered structure of the H i, affecting according to Sect. 9 multipoles $l\la 9$, and apodization to avoid confusion from the Galactic plane, our investigations are not affected by boundaries. Our global panoramic view is based on a coordinate system that is comoving with the LSR in the center of a locally turbulent flow and should not generate limitations concerning “velocity differences.”

7 VCA revisited

With the results from Sect. 5 we now apply a modified velocity channel analysis (Lazarian & Pogosyan, 2000) to our data. From the theoretical point of view VCA is quite appealing; however, from the channel cross power analysis in the previous section we see that some care is needed to avoid biases caused by a decorrelation of the H i database due to confusion that is, unfortunately, not easily recognizable from the data. From the multiphase composition of the ISM we get additional constraints that are also not recognizable without a decomposition in different phases.

7.1 Characteristic broadband spectral indices

We conclude from Fig. 14 that the best possible estimates for multiphase broadband power spectra in our case should be derived for $\Delta v_{\mathrm{LSR}}=16$ km s*-1*. This is rather narrow compared to the typical FWHM width of 23.3 km s*-1* for the WNM, but interestingly this width is close to the FWHM width of $\Delta v_{\mathrm{LSR}}=16.8$ km s*-1* for the velocity distribution of filamentary H i structures. It has been suggested that these structures indicate processes that feed kinetic energy to the local ISM (Kalberla et al., 2016, Sect. 5.13 and 5.14, Fig. 22). We show the broadband spectral distributions resulting for $\Delta v_{\mathrm{LSR}}=16$ km s*-1* in Fig. 23, and obtain remarkably well-defined results: $\gamma=-2.943\pm 0.003$ for all-sky and $\gamma=-2.944\pm 0.005$ for high latitudes. The derived uncertainties are formal one-sigma errors, but systematic fluctuations on the order of $\Delta\gamma\sim 0.05$ may be more characteristic.

Decomposing broadband distributions for individual phases leads, according to Fig. 23 for multipoles $10<l<100$, to all-sky spectral indices $-2.15>\gamma>-2.5$, comparable to the single-channel results at $v_{\mathrm{LSR}}=0$ km s*-1*, shown in Fig. 1. These results are questionable since they may be affected by spurious effects from differential Galactic rotation. At high latitudes we obtain $\gamma=-2.40$ for the WNM and $\gamma=-2.48$ for the other phases. These values differ only slightly from the single-channel result $\gamma\sim-2.36$ at $v_{\mathrm{LSR}}=0$ km s*-1*. The cross power spectral index for the WNM at high latitudes with a channel separation of $v_{\mathrm{LSR}}=10$ km s*-1* (Fig. 15) is, within the errors, identical to the auto power WNM spectral index $\gamma=-2.40$ for $\Delta v_{\mathrm{LSR}}=16$ km s*-1* (Fig. 23).

At high Galactic latitudes the spectral indices for broadband power spectra may also be determined from Fig. 13, assuming that for single H i phases the minima of the derived distributions are more characteristic. We obtain in this case for CNM and LNM $\gamma\sim-2.5$; the value for the WNM also appears to be consistent with this result, though a minimum is not yet reached at $\Delta v_{\mathrm{LSR}}=50$ km s*-1*. All-sky results according to Fig. 12 cannot be used this way; these data are obviously affected by confusion.

7.2 Modified VCA spectral index for the density field

The observational determined 3D spectral index $\gamma$, according to Lazarian & Pogosyan (2000) or Lazarian (2009, Table 3), can be used to estimate the spectral index of the density correlation function $\Gamma_{\rho}$ and for the velocity correlation function $\Gamma_{v}$, respectively. For data that are averaged over velocity the observed fluctuations in column density must be due to density fluctuations. In the case of a shallow 3D density distribution with $\Gamma_{\rho}>0$ the observed thick slice power index corresponds to $\gamma=-3+\Gamma_{\rho}>-3$; this transformation is essentially a reduction in dimensionality.

Not foreseen in the framework of VCA are limitations due to a decorrelation of the velocity field. We use for CNM and LNM the minima of the spectral index distribution from Fig. 13, both close to $\gamma=-2.5$. This value also appears reasonable for the WNM. We obtain $\Gamma_{\rho}\sim 0.5$. Using indices from Fig. 23 would result in $\Gamma_{\rho}\sim 0.6$ for the WNM and $\Gamma_{\rho}\sim 0.52$ for LNM and CNM. We disregard all-sky data since they are most probably biased by confusion from the Galactic plane. Using observed H i column densities with a 3D spectral index $\gamma\sim-2.94$ (Fig. 23 bottom right) would imply $\Gamma_{\rho}\sim 0.06$. Given all the constraints discussed previously we consider this value (or $\gamma\sim-2.94$ for the 3D density field) as characteristic of the total H i gas phase, but unrelated to the multiphase composition.

7.3 Modified VCA spectral index for the velocity field

For thin velocity slices in the case of a shallow 3D density distribution with $\gamma>-3$, the observed power index corresponds to $\gamma=-3+\Gamma_{\rho}+\Gamma_{v}/2$ (Lazarian, 2009, Table 3). Here $\Gamma_{v}$ is the spectral index of the velocity correlation function. For this relation it is assumed that intensity fluctuations in thin velocity slices are significantly affected by velocity effects. Line emission that is shifted in velocity can mimic density fluctuations, so-called caustics. On Galactic scales this effect is known as velocity crowding (Burton, 1972). For shallow density spectra (i.e., $\gamma>-3$) and independent random velocity and density distributions, the pure velocity effect should dominate the density fluctuations on large scales (Lazarian & Pogosyan 2006), where fluctuations are supposed to be large compared to the mean density. Lazarian & Yuen (2018) conclude that on scales larger than 3 pc, corresponding to multipoles $l\la 100$ at an assumed distance of 100 pc, fluctuations in channel maps are dominated by velocity fluctuations. Thus, velocity caustics should mimic real physical entities such as filaments. Clark et al. (2019) reject this interpretation and find that the H i intensity features are real density structures in a multiphase medium and not velocity caustics. Cold CNM filamentary structures exist predominantly at small scales, or $l\ga 100$, and should (at least for shallow H i power spectra) be largely unaffected by velocity fluctuations.

Thus, according to Lazarian & Pogosyan (2000) and Lazarian (2009) $\Gamma_{v}/2$ can be determined from changes between thin and thick velocity slice power distributions. We use results at high latitudes shown in Fig. 13 and find changes in the power law index in the range 0.12 to 0.16. Accordingly, we estimate $\Gamma_{v}\sim 0.28$ for the individual H i phases. In the case of total H i column densities we obtain from Fig. 14, $\Gamma_{v}\sim 0.3$. These results are low in comparison to the Kolmogorov (1941) index $\Gamma_{v}=2/3$, but intermittency and supersonic motions may modify this index (Kolmogorov, 1962). For a compressible ISM the spectral indices may depend on the sonic Mach number (e.g., Burkhart et al., 2010). A 3D spectral index of $\gamma\sim-2.5$ for $10\la l\la 100$ implies accordingly Mach numbers near eight, though from the Gaussian decomposition lower values are expected: $M_{\mathrm{CNM}}=4.4$, $M_{\mathrm{LNM}}\la 3.6$, and $M_{\mathrm{WNM}}\sim 1.4$ (Kalberla & Haud, 2018).

Our estimates of $\Gamma_{v}$ according to a modified VCA are lower limits only. The HI4PI survey data are limited in velocity resolution to $\delta_{v}=1$ km s*-1* for GASS and $\delta_{v}=1.5$ km s*-1* for EBHIS, the asymptotical thin slice limit is not reached. At $\Delta v_{\mathrm{LSR}}=0$ km s*-1* we should see in Figs. 12 and 13 a horizontally tangential approach. Our resolution is insufficient to derive a significant thin slice velocity limit. The demanding condition $\delta_{v}\ll 1$ km s*-1* for a narrow line width is quite a surprise since the estimate by Lazarian & Pogosyan (2000, Sect. 4.3) was that a thin slice of $\delta_{v}\sim 2.6$ km s*-1* should be a sufficient velocity slice width for the H i. Very thin slices in velocity may cause serious observational problems; the instrumental noise would be significantly increased and the noise term $N_{l}$ then needs to be explicitly taken into account.

7.4 Phase transitions versus velocity caustics

In the previous subsection we derived VCA predicted spectral index changes with constraints on spectral coherence from Sect. 5. Here we explore how the VCA solution depends on multipole $l$. In particular we want to derive how much the CNM power distribution $P_{\mathrm{CNM}}(l)$ is modified by bandwidth changes $\Delta v_{\mathrm{LSR}}$ in velocity. We define the ratio

[TABLE]

for the narrowband (single-channel) power $P_{\rm narrow}(l)$ and broadband power $P_{\rm broad}(l)$, integrated over the CNM coherence width $\Delta v_{\mathrm{LSR}}=10$ km s*-1*, both at identical center velocities.

We chose two velocity settings. In the first case, for $v_{\mathrm{LSR}}=0$ km s*-1*, there are significant fluctuations in the CNM spectral index. The second velocity, $v_{\mathrm{LSR}}=-10$, was chosen such that the CNM spectral index has little fluctuations (see Fig. 9). We show $\mathfrak{R}_{\rm VCA}(l)$ for $v_{\mathrm{LSR}}=0$ km s*-1* in Fig. 24 and for $v_{\mathrm{LSR}}=-10$ in Fig. 25. In each case we observe more power than expected from the increased bandwidth, $\mathfrak{R}_{\rm VCA}(l)>10$. There are in addition two multipole ranges with different slopes in log-log presentation.

For $l\la 100$ we observe a steepening of the spectral index with $\delta\gamma=-0.12\pm 0.04$ at $v_{\mathrm{LSR}}=0$ km s*-1* (Fig. 24) and $\delta\gamma=-0.21\pm 0.04$ at $v_{\mathrm{LSR}}=-10$ km s*-1* (Fig. 25). This steepening is expected from the velocity correlation function (Lazarian & Pogosyan (2000) and Lazarian (2009)), and we derive accordingly $\Gamma_{v}\sim 0.24$ at $v_{\mathrm{LSR}}=0$ km s*-1* and $\Gamma_{v}\sim 0.4$ at $v_{\mathrm{LSR}}=-10$ km s*-1*, consistent with the results from the previous subsection.

For $l\ga 100$ we observe the opposite trend, i.e., the power increases for high multipoles. The sign for $\delta\gamma$ is positive, but the distribution cannot be approximated with a constant spectral index. Could this multipole range be biased? According to Lazarian & Pogosyan (2000) the velocity slice thickness must be increased to steepen the power law distribution; however, this approach leads to contradictions. For an increase in slice thickness beyond $\Delta v_{\mathrm{LSR}}\sim 10$ km s*-1* we observe for low mutilpoles a flattening of the spectral index (see Sect. 5, Fig. 13). The power at high multipoles increases further. The power for $l\ga 100$ cannot be explained as being caused by velocity caustics. It must originate from phase transitions, causing extra power at high multipoles as discussed in Sect. 5.5.

The question arises whether velocity caustics or phase transitions are dominating spectral index changes. The results from Sect. 4 indicate that spectral index changes from phase transitions can amount to $\delta\gamma\sim-0.4$, twice as large as changes determined in this section and by a VCA analysis in velocity width in the previous section. Such changes are usually attributed to velocity caustics. CNM structures at $l\sim 1000$, discussed in Sect. 5.5, need to be interpreted as magnetized dust-bearing density structures (Clark et al., 2019). These entities are organized in larger filamentary structures, up to scales of tens of degrees (Clark et al. (2014) and Kalberla et al. (2016)); the most prominent structure is Loop I. Contemporary observations of such objects with large single-dish telescopes (Arecibo, GBT, Effelsberg, Parkes), also with the DRAO interferometer, interpret these structures consistently as clouds or filamentary cloud complexes, real density structures with a well-defined range of physical parameters, (e.g., Clark et al. (2014), Martin et al. (2015), Kalberla et al. (2016), Blagrave et al. (2017), and Clark et al. (2019)). An interpretation that most of the structures should be due to velocity caustics (Lazarian & Yuen, 2018) is not supported by observations.

8 Intermediate and high velocity clouds

Power spectra for thick velocity slices, discussed in the previous section, may be affected by emission from IVCs for large $\Delta v_{\mathrm{LSR}}$. We have chosen the most prominent IVC emission in the velocity range $-70<v_{\mathrm{LSR}}<-30$ km s*-1* and calculated the power distribution. The IVC emission is dominant on the northern hemisphere. To avoid any possible instrumental biases from a telescope mix we extracted EBHIS data, using the apodization schemes discussed in Sect. A.2.2. Our result is shown in Fig. 26. The all-sky data are obviously affected by spurious effects from differential Galactic rotation and do not represent genuine turbulent features. High latitude data show a straight power spectrum with an index $\gamma=-2.620\pm 0.004$, within possible systematical uncertainties $\Delta\gamma\sim 0.05$ in good agreement with $\gamma=-2.68\pm 0.04$ obtained by Martin et al. (2015), and $-2.60\pm 0.04\la\gamma\la-2.48\pm 0.06$ byBlagrave et al. (2017). In comparison to a thick slice index of $\gamma=-2.94$ for the local H i, this is significantly flatter. No decomposition of the ICV emission in components from different H i phases is applied since such a separation for sources with unknown distances and beam smoothing effects appears to be ambiguous.

For the sake of completeness we also calculated power spectra for HVCs in the velocity range $-150<v_{\mathrm{LSR}}<-100$ km s*-1*. In this case as well we only used data from the northern hemisphere. The resulting power spectrum in Fig. 27 is useless in the all-sky case, but at high latitudes it is straight for $2<l<256$ with $\gamma=-2.0$. This is even flatter than the index $\gamma=-2.59\pm 0.07$ observed by Martin et al. (2015) and $\gamma=-2.85\pm 0.07$ by Blagrave et al. (2017). Our power spectrum shows for multipoles $l>350$ an S-shaped feature that cannot be explained by instrumental biases. For $350\la l\la 700$ there is a steepening to $\gamma\sim-3.2$ with an abrupt flattening at high multipoles. This S-shaped feature comes from high latitude HVCs and is even more pronounced if we apodize data for |b|>30\hbox{{}^{\circ}}. An interpretation without further detailed investigations, which are beyond the scope of the current paper, is difficult. Perhaps this S-feature is caused by fragmentation and dissipation of HVCs on scales near 1∘. The observed differences in spectral indices between our data and those of Martin et al. (2015) and Blagrave et al. (2017) may imply changes between different HVC complexes.

9 Low multipoles – outer scale

It is broadly believed that the outer scale $L$, first mentioned in Sect. 3 and most likely defined by energy injection from old supernova remnant shock waves, must be close to $L\sim 100$ pc (Haverkorn et al., 2008). According to Chepurnov (1998), Cho & Lazarian (2002), and Mertsch & Sarkar (2013) the critical multipole $l_{\mathrm{crit}}$ for a broken power law depends on $L$ and the scale height $H$ of the turbulent medium, $l_{\mathrm{crit}}\sim 2\pi H/L$. The half width at half maximum scale height for the H i layer was determined to $115\la H\la 140$ pc (Dickey & Lockman 1990, Fig. 10 and Kalberla et al. 2007, Fig. 14), resulting in $7\la l_{\mathrm{crit}}\la 9$. This is in good agreement with our finding that the power spectra tend, within the noise, to be shallow at $l\la 9$. Strong multipole disparities are also mostly restricted to $l\la 9$. (e.g., Figs. 1 and 33).

Low multipoles $l\leq 8$ were in all cases fit separately from the rest of the power distribution, using in this case constant weights. Figure 28 shows the derived velocity dependent power law indices. As an average over all fits for $-16<v_{\mathrm{LSR}}<16$ km s*-1* we determine $\gamma=-0.66\pm 0.13$ for the all-sky data and $\gamma=-0.65\pm 0.09$ for |b|>20\hbox{{}^{\circ}}.

Our estimate $l_{\mathrm{crit}}\sim 8$ is consistent with Regis (2011) who determined $l_{\mathrm{crit}}=10\pm 3$ for the all-sky Galactic continuum emission from five radio maps and $l_{\mathrm{crit}}=5^{+10}_{-4}$ at high latitudes. Our power law index $\gamma=-0.65\pm 0.09$ for $l\la l_{\mathrm{crit}}$ appears to be in conflict with the proposed value $\gamma=-1$. However, we need to consider that the index, derived for $l_{\mathrm{crit}}\sim 8$, is probably influenced by the planar distribution of the Galactic H i. On large scales the column density distribution is dominated by cosecant latitude effects. Removing this effect leads to a change in the power distribution at low multipoles (see Sect. 3.3, and Fig. 4 in comparison to Fig. 1).

We use all-sky data, but our main conclusions are drawn from high latitudes, |b|>20\hbox{{}^{\circ}}. This region covers 67% of the sky with typical H i column densities below $10^{21}\mathrm{cm}^{-2}$. We compare this region with that used by Planck intermediate results XXX. (2016). To study the angular power spectrum of polarized dust emission at intermediate and high Galactic latitudes, sky fractions between 30% and 80% were investigated by these authors to derive best-fit power law exponents. For multiploes $l>40$ there are no significant changes in the best-fit power law exponents, in particular no changes for sky fractions of 70% and 80%. From these results, since gas and dust at high Galactic latitudes are well correlated with each other (Planck early results. XXIV., 2011), we expect no significant cosecant latitude effects on our analysis of H i column densities for multipoles $l\ga 10$.

10 Summary and discussion

We use high resolution 21 cm line data from the HI4PI survey (HI4PI Collaboration et al., 2016) to determine spatial power spectra for the H i in the local ISM. These panoramic observations enable us to derive spatial power spectra for multipoles $9\la l\la 1023$. The power distribution for multipoles $l<9$ is consistent with the assumption that turbulence has an outer scale close to 100 pc. Most of the observed multiphase H i power spectra are exceptionally straight in log-log presentation and can be fit well with constant spectral indices, $\gamma\sim-2.94$ for the local gas, $\gamma\sim-2.6$ for IVCs, and $\gamma\sim-2$ for HVCs. HVCs deviate from this rule and show some steepening at high multipoles. HI4PI power spectra for the emission from the complete Galactic plane can be contaminated seriously by confusion from differential Galactic rotation. H i data for several low latitude subfields within a range 98\hbox{{}^{\circ}}\la l\la 145\hbox{{}^{\circ}} and $-20.4<v_{\mathrm{LSR}}<0.2$ km s*-1* were analyzed previously by Khalil et al. (2006). They derived $\gamma=-2.95\pm 0.05$, in good agreement with our result for the local gas at high latitudes.

The extraordinarily straight power distribution for the total H i at high latitudes with its well-defined broadband index of $\gamma\sim-2.944\pm 0.005$ (with possible systematic fluctuations on the order of $\Delta\gamma\sim 0.05$; see Sect. A.2.2) is in excellent agreement with dust spectral index $\gamma_{\mathrm{dust}}=-2.9\pm 0.1$ determined by Miville-Deschênes et al. (2016) who also conclude from the straight power distribution that there is no sign of energy dissipation at the 0.01 pc scale. The power spectrum for the dust has to be interpreted as a spectrum of turbulent density fluctuations. Gas and dust are well mixed. The value $\gamma\sim-2.94$ for the total (phase independent) H i distribution must therefore be interpreted as characteristic of the total H i density distribution.

Phase transitions lead to modified distribution laws for individual phases (Saury et al., 2014). The observed filamentary CNM structures are well aligned with the magnetic field (Clark et al. (2014) and Kalberla et al. (2016)), resulting in asymmetries (Kalberla & Kerp (2016) and Kalberla et al. (2017)). The 3D space is no longer populated isotropically (Goldreich & Sridhar, 1995). This probably implies a reduction in dimensionality, suggesting a flattening of the 3D power spectra for those parts of the H i that share the magnetic field orientation, the CNM, and the LNM.

We derive single-phase power spectra for CNM, LNM, and WNM from a Gaussian decomposition with characteristic linewidth regimes. These individual phases show power spectra that deviate significantly from the multiphase power spectra. For all of the phases we get straight power spectra for $10<l<100$ with indices $\gamma\sim-2.5$. In addition for $l\ga 100$ a strong increase in the power up to high multipoles is observed. The single-phase power spectra are, according to Eq. 4, highly correlated; the auto power relations for all phases must be considered together with all cross power terms between these phases. This correlation leads to counterintuitive consequences and contradictions if we try to describe the H i with several independent phases.

We find distinct correlations in velocity space. Coherence in the CNM is limited to a velocity spread of $\sim 10$ km s*-1*. Coherent features at high multipoles are identified as cold clumps with median Doppler temperatures of 223 K. Our analysis recovers structures that were previously described by Clark et al. (2014) and Kalberla et al. (2016). Trying to interpret these CNM structures as velocity caustics leads to contradictions. In particular, according to VCA the spectral index should steepen; however, we observe it to be undefined for velocity widths beyond 10 km s*-1*. These CNM features must be real entities in accordance with the findings of Clark et al. (2019). H i filaments on large scales are made up of aligned small-scale CNM structures. These CNM features show velocity fluctuations on larger scales, but these large-scale structures cannot be explained as being caused by velocity caustics.

Phase transition must be considered as local events with limitations in the cloud-cloud velocity dispersion for the resulting CNM clumps. For different bulk velocities we obtain coherence for different samples of CNM objects. The spectral coherence of the CNM is self-similar and statistically invariant under translation in velocity. The correlation between CNM and LNM implies that the LNM exists in connection with the CNM as a more extended phase, enveloping the CNM spatially but also covering a wider velocity spread of $\sim 23$ km s*-1* as a transition to the more extended WNM. Self-similarity in velocity also applies to the LNM. Hence, associated with the hierarchical wave number scaling for phases with different linewidth domains we obtain a remarkable homogeneity in velocity space. Comparing power spectra for different hemispheres we find a high degree of isotropy.

The remarkable loss of correlation at large velocity lags is not in contradiction to the theoretical setup for the analysis of turbulent flows as considered by Kolmogorov (1941). The treatment of a turbulent flow needs to be restricted to a sufficiently small domain in four-dimensional space. Homogeneity and isotropy was explicitly defined by Kolmogorov as “local” only. These constraints, as described by Kolmogorov, correspond to the coherence conditions derived by us.

Restrictions in the velocity domain imply constraints for a velocity analysis, including the use of velocity centroids and explain discrepancies discovered previously by Kalberla et al. (2017, Appendix C). Our multiphase 3D spectral index $\gamma=-2.94$, derived with velocity constraints, is in conflict with steeper power law indices derived by Dickey et al. (2001) and Stanimirović & Lazarian (2001); however, these authors considered very different conditions. Dickey et al. (2001) got thick slice 3D spectral indices up to $-4.0$ in the Galactic plane for velocities $-20<v_{\mathrm{LSR}}<-80$ km s*-1*, covering distances $1.5\la d\la 13$ kpc in the inner Galaxy. Stanimirović & Lazarian (2001) found a thick slice 3D spectral index of $-3.4$ for the Small Magellanic Cloud, analyzing structures between 30 pc and 4 kpc perpendicular to the line of sight, integrating in velocity over $\Delta v_{\mathrm{LSR}}=90$ km s*-1*.

We use VCA (Lazarian & Pogosyan, 2000) modified with restrictions in velocity space to disentangle density and velocity effects on observed power spectra. Taking velocity decorrelations into account, we derive single-phase estimates of $\Gamma_{\rho}\sim 0.5$ for the spectral index of the density correlation function and $\Gamma_{v}\sim 0.3$ for the spectral index of the velocity correlation function for $10\la l\la 100$. For the total observed H i column density distribution we determine $\Gamma_{\rho}\sim 0.06$ and $\Gamma_{v}\sim 0.3$ for all multipoles $l\ga 9$. The significance of the derived $\Gamma_{v}\sim 0.3$ must be questioned. A steepening of the power spectra as expected for velocity caustics can also be caused by phase transitions. For $v_{\mathrm{LSR}}=3$ km s*-1* we determine a CNM spectral index change $\delta\gamma\sim-0.4$ from phase transitions, larger than the change $\delta\gamma\sim-0.14$, indicated by VCA.

Spectral indices for single-channel power spectra are steepest at velocities close to zero km s*-1* where CNM phase fractions are high and WNM phase fractions are low. We interpret this as an indication of local phase transitions. Phase transitions do not change the total column densities significantly, but decrease locally the line widths. Thus, thermal instabilities produce enhanced power in the line centers at the expense of the power in the wings of the line. At the same time phase transitions cause fragmentation, small-scale structure is generated, giving rise to enhanced power at high multipoles ($l\ga 100$). This process is necessarily accompanied by a decrease in power on scales that are characteristic of the regions that were affected by phase transitions and fragmentation. Power is reduced on intermediate scales and the power spectrum steepens for scales that are little affected by fragmentation ($l\la 100$).

Recent high resolution simulations of thermal instabilities and collapsing cold clumps by Wareing et al. (2019) show that the density power spectrum can rapidly rise and steepen as the structures grow in the simulations. They find sheets and filaments on typical scales of 0.1 to 0.3 pc, consistent with high power for the CNM at large multipoles and also consistent with small-scale structures observed by Clark et al. (2014) and Kalberla et al. (2016). Previously Audit & Hennebelle (2005) and Federrath (2016) have derived a preferred scale of 0.1 pc from simulations. Our results appears to be in conflict with Saury et al. (2014) who find for subsonic turbulence that phase transitions lead to shallow power spectra with $\gamma\sim-2.4$. The discrepancy would be resolved if their analysis were sensitive to the single-phase CNM for which we observe $\gamma\sim-2.37$ at high latitudes.

Steep broadband 3D power law indices $\gamma\sim-3.6$ that were previously considered to be characteristic for H i in emission deviate according to Hennebelle & Falgarone (2012, Fig 10) significantly from power law indices for molecular lines and other data. The single-phase power indices $\gamma\sim-2.5$, derived by us, are in much better agreement with the other indices. This also applies to the comparison by Ghosh et al. (2017, Fig. 11). Cold filamentary H i structures were found to be correlated with *Planck * 353 GHz polarized dust emission (Kalberla et al., 2016). Our broadband CNM power index of $\gamma=-2.48\pm 0.04$ fits well to the mean *Planck * intermediate latitude sky exponent of $\gamma=-2.42$ at 353 GHz for the EE and BB polarized dust power spectra at multipoles $40<l<370$ (Planck intermediate results XXX., 2016). Ghosh et al. (2017) use in their dust model for the south Galactic pole region the CNM as a tracer of the dust polarization angle, and find for EE, BB, and TE modes spectral indices close to $-2.4$, a value that they consider to be characteristic of the turbulent magnetic field as well. This is consistent with the results from Vansyngel et al. (2017) and more recently with Planck 2018 results. XI. (2018), who obtain an index of $-2.42$ for the EE, $-2.54$ for the BB, and $-2.49$ for the TE mode. Adak et al. (2019) obtain for their dust model for the north Galactic pole region an index of $-2.4\pm 0.1$ for EE, BB, and TE modes and the turbulent magnetic field, in excellent agreement with our value of $\gamma=-2.48$ for CNM and LNM at high latitudes.

Acknowledgements.

We acknowledge the second anonymous referee for the constructive criticism. P.K. thanks Jürgen Kerp for the inspiring discussions and Bärbel Koribalski for the help with the LVHIS data. U.H. acknowledges the support from the Estonian Research Council grant IUT26-2, and from the European Regional Development Fund (TK133). This research has made use of NASA’s Astrophysics Data System. EBHIS is based on observations with the 100m telescope of the Max-Planck-Institut für Radioastronomie (MPIfR) at Effelsberg. The Parkes Radio Telescope is part of the Australia Telescope, which is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO. Some of the results in this paper have been derived using the HEALPix package.

Appendix A Observations and data reduction

A.1 HI4PI 21-cm line survey data

We use HI4PI data (HI4PI Collaboration et al., 2016), combining 21 cm line data from Effelsbeg-Bonn H i Survey (EBHIS), observed with the 100m Effelsberg radio telescope (Winkel et al. (2016) and Winkel, Lenz & Flöer (2016)) and from the Galactic All Sky Survey (GASS), observed with the 64m Parkes telescope (McClure-Griffiths et al. (2009), Kalberla et al. (2010), and Kalberla & Haud (2015)). EBHIS data from the first data release (Winkel, Lenz & Flöer, 2016) and GASS data from the final data release (Kalberla & Haud, 2015) were calibrated to a common intensity scale.

The multi-beam receivers scanned the sky, dumping spectra with short integration times (0.5 sec for EBHIS and 5 sec for GASS). Each of the individual dumps was corrected for stray radiation. The telescope data were gridded using Gaussian kernels of 5$.\mkern-4.0mu^{\prime}$4 for EBHIS and 7$.\mkern-4.0mu^{\prime}$5 for GASS, resulting in an effective FWHM beam size of 10$.\mkern-4.0mu^{\prime}$8 for EBHIS, and 14$.\mkern-4.0mu^{\prime}$5 for GASS. An nside = 1024 HEALPix database was chosen, appropriate for an angular resolution of \Theta_{\mathrm{pix}}=3\hbox{.\mkern-4.0mu^{\prime}}44 (Górski et al., 2005). This implies that our database is oversampled; neighboring positions are not independent from each other. Bandwidth limitations need to be taken into account, but a deconvolution (discussed below) for nside = 1024 is rather simple and well defined. For both telescopes a large number of dumps from the multi-beam receiver were averaged to generate profiles at the HEALPix grid positions, leading in each case to well-defined Gaussian beam shapes.

Combining data from different hemispheres we use a border line at \delta=-2\hbox{{}^{\circ}}. To avoid discontinuities caused by different beam sizes, we take a linear interpolation between the two surveys for -4\hbox{{}^{\circ}}<\delta<0\hbox{{}^{\circ}}.

The 21 cm line surveys were observed with different spectrometers. For EBHIS the spectral resolution is $\delta v_{\mathrm{LSR}}=1.29$ km s*-1*; for GASS it is $\delta v_{\mathrm{LSR}}=0.82$ km s*-1*. The average noise for a single channel outside the H i line is 90 mK for the EBHIS and 55 mK for the GASS. We chose to interpolate the data to a common intermediate velocity grid with $\delta v_{\mathrm{LSR}}=1.03$ km s*-1*, which is the velocity resolution of the LAB (Kalberla et al., 2005); the same database has been used previously (Kalberla et al., 2016). This interpolation causes for the EBHIS a slight increase in the average noise to 100 mK.

A.2 Instrumental effects

The observed intensities $\widetilde{I_{\mathrm{obs}}}$ can for individual telescope dumps be described by the convolution between the true distribution $\widetilde{I_{\mathrm{sky}}}$ on the sky and the beam function $\widetilde{B_{\mathrm{tel}}}$ of the telescope. In addition, we need to take into account the instrumental noise $\widetilde{N_{\mathrm{dump}}}$, which affects individual dumps:

[TABLE]

To get the intensity distribution on a HEALPix grid, the dump database is interpolated using a Gaussian weighting with an smoothing kernel $\widetilde{B_{\mathrm{grid}}}$. This implies an additional convolution and we obtain

[TABLE]

For convenience, we define here the smoothing function $\widetilde{B_{\mathrm{data}}}=\widetilde{B_{\mathrm{grid}}}\ast\widetilde{B_{\mathrm{tel}}}$ that was applied to the data.

To calculate the power distribution we need to autocorrelate and perform a Fourier transform on Eq. 13, leading to

[TABLE]

with the cross power $P_{\mathrm{cross}}$ from $\widetilde{B_{\mathrm{data}}}\ast\widetilde{I_{\mathrm{sky}}}$ and $\widetilde{B_{\mathrm{grid}}}\ast\widetilde{N_{\mathrm{dump}}}$

[TABLE]

or with $B_{\mathrm{data}}=B_{\mathrm{grid}}\ B_{\mathrm{tel}}$ simplified to

[TABLE]

This result makes it necessary to correct the beam smoothing with different kernel functions for observed data and the noise term, and also to consider the question of whether the observations are affected by a cross term between noise and data. Smoothing and noise effects are discussed below in separate subsections. We also need to consider that EBHIS and GASS have different beam functions and different noise contributions. In data processing we have to process the two different hemispheres separately.

A.2.1 Apodization

If we want to analyze only a part of the sphere, in our case high latitudes with |b|>20\hbox{{}^{\circ}}, we need to take the window function into account that masks the observations. Since this window is applied to the brightness temperature distribution on the plane of the sky, it implies a smoothing to the observed power distribution $P_{\mathrm{obs}}(l)$. Sharp edges of the window function may cause aliasing or power leakage between different multipoles $l$. To avoid such aliasing it is necessary to smooth the sharp edges of the window function by apodizing (Planck 2015 results, XI, 2016).

Figure 29 shows the raw power spectrum for the observed column densities in a single channel at $v_{\mathrm{LSR}}=0$ km s*-1*. We use this case here to demonstrate our data processing. The upper black profile presents the all-sky power spectrum; below we show two versions of the power spectrum at high latitudes, |b|>20\hbox{{}^{\circ}}. The cyan spectrum in Fig. 29 shows $P_{\mathrm{obs}}(l)$ after setting in ANAFAST the parameter theta_cut_deg (or gal_cut in healpy) to 20∘. We observe strong aliasing caused by an observational window with sharp cutoff at 20∘; an apodization is needed to overcome spurious effects.

Similar to Kalberla & Kerp (2016) we use a cosine taper (a Tukey window, Harris, 1978) for apodization in latitude $b$. After a few tests we chose a taper width with a half-period of 15∘, weighting smoothly from one at |b|=20\hbox{{}^{\circ}} to zero at |b|=5\hbox{{}^{\circ}}. The resulting power spectrum is shown in Fig. 29 in red, overplotted on the unapodized case. Testing a narrower cosine taper, for example with a half-period of 10∘, we found still significant unwanted aliasing effects. A broader cosine taper for |b|=20\hbox{{}^{\circ}} is also not appropriate since it leads to a leakage of Galactic plane emission from |b|<5\hbox{{}^{\circ}}. Any apodization window leads to unwanted side effects (Harris, 1978), but this apodization-ringing decreases as the width of the window increases. Remaining problems (Fig. 29, red curve) are acceptable for our analysis with minor errors, as demonstrated in Fig. 4. The optimal width of the Tukey window that we applied is slightly larger then the \sigma=5\hbox{{}^{\circ}} Gaussian taper used for the apodization map generated by the cosmological parameters team (Planck 2015 results, XI, 2016).

A.2.2 Beam smoothing

Observed brightness temperatures are smoothed by the beam function. In addition we need to take into account that for each cell of the HEALPix grid the data are interpolated. The resulting effective beam, has a FWHM beam width of 10$.\mkern-4.0mu^{\prime}$8 for EBHIS (Winkel et al., 2016). For the GASS we needed to redetermine the FWHM beamwidth since we use the third data release (Kalberla & Haud, 2015) with some improvements over the second release Kalberla et al. (2010). We searched the Local Volume H i Survey (Koribalski et al., 2018) for appropriate point-like narrowband H i sources and used the barred Magellanic irregular galaxy IC 4662 (van Eymeren et al., 2010) which has at $v_{\mathrm{LSR}}=324$ km s*-1* a bright compact 2′ core. In addition, the beam width was determined from the ultra-compact high velocity cloud HVC289+33+251 (Putman et al. (2002) and Brüns & Westmeier (2004)) with a FWHM size of 4$.\mkern-4.0mu^{\prime}$4. Fitting the sizes of these sources we derive a FWHM beamwidth of 14\hbox{.\mkern-4.0mu^{\prime}}5\pm 0\hbox{.\mkern-4.0mu^{\prime}}2 which compares well with the beamwidth of 14$.\mkern-4.0mu^{\prime}$4 determined by Calabretta et al. (2014). This 14\hbox{.\mkern-4.0mu^{\prime}}5 beam is a considerable improvement over the second GASS data release with an effective beam width of 15$.\mkern-4.0mu^{\prime}$9 in the presence of residual instrumental errors.

We consider once more a single channel at $v_{\mathrm{LSR}}=0$ km s*-1* (see Fig. 29). The upper all-sky power spectrum (black) and the lower one (red and apodized for |b|=20\hbox{{}^{\circ}}), both suffer for $l\ga 256$ from beam smoothing. We observe a steepening of the spectral index at high multipoles, which without beam correction may be mistaken as a real effect, for example as steepening caused by a change in slice thickness (Lazarian & Pogosyan, 2000). For the raw data shown in Fig. 29 essentially only a limited range, $l<256$, as indicated by a vertical bar, can be used to fit the spectral index $\gamma$. It is essential to correct for beam smoothing if we want to avoid instrumental biases.

To evaluate beam smoothing effects we separate the emission observed in different hemispheres. We use a 15∘ cosine apodization scheme, extracting data for EBHIS at declinations \delta>13\hbox{{}^{\circ}} and for GASS at \delta<-17\hbox{{}^{\circ}}. Figure 30 (top) shows the correlation coefficients for EBHIS (left panel) and GASS (right panel). The black lines show significant differences in beam smoothing caused by different beam widths. We correct in both cases the observed power spectra by applying corrections based on the effective FWHM beam widths, $P=P_{\mathrm{obs}}/B^{2}_{\mathrm{data}}$ according to the first term in Eq. 15. $B_{\mathrm{data}}(l)$ is telescope dependent and includes beam smoothing and smoothing caused by gridding. The power spectra, including their fits for $8<l<1023$ are shown in magenta. A comparison between observed and beam corrected power spectra shows that differences can be seen for $l\ga 256$, corresponding to angular scales of \theta\sim 180\hbox{{}^{\circ}}/l\sim 42\hbox{{}^{\prime}} or about three to four FWHM beam widths. In the case of noisy data, the beam smoothing effects may not be recognizable until scales corresponding to one beam width. This does not imply that beam corrections are not necessary.

Merging disjunct data from EBHIS and GASS makes it necessary to correct at once for beam effects from both telescopes. For an isotropic power distribution this is fortunately very simple; at any multipole $l$ the power from the two telescopes is summed; we can correct for beam effects by using a weighted average beam function, taking the borderline at \delta=-2\hbox{{}^{\circ}} into account. The lower left plot in Fig. 30 can serve as a proof. Finally, we compare the corrected power spectra in the lower panel to the right side of Fig. 30. The spectra from EBHIS and GASS agree well; the sum of the two power distributions is shown in red. This sum contains apodized regions; therefore, we also plot the unapodized all-sky power spectrum. This comparison shows a nearly isotropic power distribution between the two disjunct hemispheres, but there are systematic changes in the spectral index. We obtain $-2.717>\gamma>-2.810$ with an average of $\gamma=-2.768$, hence a deviation of up to $\Delta\gamma\sim 0.05$ from the all-sky index. This is more than a factor of ten larger than the formal one-sigma uncertainty for an individual fit in the range for $8<l<1023$. The change of the spectral index reflects changes in the source distribution and is predominantly caused by differences at low multipoles.

A.2.3 System noise

It has been demonstrated by Kalberla & Kerp (2016) that a beam correction according to Eq. 2 leads in general to an amplification of instrumental uncertainties at high multipoles. However, for $l\la 1023$ we do not observe an obvious increase in Fig. 30.

To check whether our analysis might be degraded by unrecognized system noise, we generate all-sky random noise maps as expected from the system performance of the two telescopes. The system noise $T_{\mathrm{sys}}$ contains several contributions. Most important is the thermal noise from the receiver system and the elevation dependent ground radiation including spill-over. Next there is the position dependent continuum sky background. All these components are variable, but we use here an average thermal contribution $T_{\mathrm{sys}}=30$ K. The line signal $T_{\mathrm{B}}(v_{\mathrm{LSR}})$ adds to the frequency independent part and the noise contribution $T_{\mathrm{Noise}}(v_{\mathrm{LSR}})$ to the H i data can then be approximated as (Haud, 2000)

[TABLE]

where $\sigma_{\mathrm{av}}$ is the average noise level, determined at velocities without H i line emission.

Figure 31 shows a noise power spectrum for such a simulation of the observational conditions for $T_{\mathrm{B}}(v_{\mathrm{LSR}})\sim 0$ K after deconvolution for the telescope beams at velocities without H i line emission. We overlay this estimate with the power distribution $P_{\mathrm{obs}}(l)/B^{2}_{\mathrm{data}}(l)$ derived for a single channel observed at the velocity $v_{\mathrm{LSR}}=400$ km s*-1*, without significant H i emission. The cross term $2P_{\mathrm{cross}}$ in Eq. 16 is therefore zero. Comparing the observed power distribution with the expected distribution, we find agreement for $l\ga 300$ but significant deviations at low multipoles, typically for single-dish telescopes (e.g., Martin et al. 2015, Figs. 6 and 9). Assuming a constant $\sigma_{\mathrm{av}}$ is only a poor approximation of the observed noise power (see Kalberla et al. 2010, Fig. 9 and Winkel et al. 2016, Fig. 14). The observational setup causes unavoidable deviations from a random noise, in particular blocky structures, which in turn lead to enhanced power at low multipoles. To improve the noise template we need to replace the averages $\sigma_{\mathrm{av}}$, $T_{\mathrm{sys}}$, with position dependent values in Eq. 17. For each of the telescope dumps and HEALPix positions the rms fluctuations in the baseline regions are known and can also be used. This more sophisticated scheme was used to control the Gaussian decomposition (Haud (2000), Kalberla & Haud (2015), and Kalberla & Haud (2018), but here we found that it was not necessary to elaborate Eq. 17 further.

We apply our simple noise model to the single channel H i data at $v_{\mathrm{LSR}}=0$ km s*-1* shown in Fig. 29. In comparison to Fig. 31 the line emission term $T_{\mathrm{B}}(v_{\mathrm{LSR}})$ in Eq. 17 is significant, leading to an enhancement of the noise power. In addition, $T_{\mathrm{B}}(v_{\mathrm{LSR}})\neq 0$ K implies that the noise depends on the intensity distribution on the sky, and $2P_{\mathrm{cross}}/B^{2}_{\mathrm{data}}$ from Eq. 16 needs to be taken into account. In Fig. 32 we compare the estimated noise contribution and the cross power with the beam corrected power spectrum. The noise contribution $P_{\mathrm{Noise}}/B^{2}_{\mathrm{tel}}$ is critical at high multipoles, but due to the good signal-to-noise ratio of the HI4PI data this noise power is, in the worst case, at $l=1023$ more than an order of magnitude below the H i signal $P_{\mathrm{obs}}(l)/B^{2}_{\mathrm{data}}(l)$. The noise cross power $2P_{\mathrm{cross}}/B^{2}_{\mathrm{data}}$ scales with $P_{\mathrm{obs}}(l)/B^{2}_{\mathrm{data}}(l)$, but is about three orders of magnitude below the relevant H i power signal. Altogether, this is a very comfortable result, reflecting that the HI4PI survey has a good signal-to-noise ratio.

For H i in the profile wings considered by us, for example at $|v_{\mathrm{LSR}}|=16$ km s*-1*, the noise power is still more than a factor of five below the line signal. Our noise template according to Eq. 17 can be applied to different velocity widths and we obtain in this case a simple scaling relation. Averaging $n$ channels decreases the noise power contribution by a factor of $n$, but at the same time the H i power increases. For a total bandwidth of $\Delta v_{\mathrm{LSR}}=16$ km s*-1* the expected noise contribution is, in the worst case, more than four orders of magnitude below the observed signal. Noise limitations, including other remaining instrumental uncertainties, are not important for HI4PI data as used by us. For most of our analysis we can safely disregard the noise term $N_{\mathrm{oise}}(l)$ in Eq. 2.

A.3 Intrinsic uncertainties in parameter fitting

We learn from the noise simulations that the scatter of $C_{\mathrm{l}}(l)$ is for white noise increasing at low multipoles $l$ (see Fig. 31). This is a general property of angular power spectra with finite counts (Campbell, 2015, Eq. 5) and we expect for the estimator $\hat{C_{l}}$

[TABLE]

$Var(\hat{C_{l}})$ is usually called cosmic variance. In Fig. 33 we show an example of single-channel data at $v_{\mathrm{LSR}}=0$ km s*-1*. When fitting power spectra we take dispersions $\propto\sqrt{C_{l}^{2}\cdot 2/(2l+1)}$ into account. Fitting is done iteratively with constant weights for an initial estimate of $\hat{C_{l}}$, next using derived weights to improve the error estimate for the final fit with $\chi^{2}/N_{\mathrm{DOF}}\la 1$. We quote for the power spectral indices formal one-sigma uncertainties from the fit. Systematic fluctuations between different hemispheres amount to $\Delta\gamma\sim 0.05$ (Sect. A.2.2) and may be more characteristic for the overall uncertainties. Differences between the initial estimate and the final fit are typically less than or about one sigma. In general we find a large scatter for $C_{l}$ at low multipoles $l\la 8$ with tendencies for even-odd oscillations. Large-scale symmetries of the Galactic H i emission cause an enhanced disparity between even and odd modes (Mertsch & Sarkar, 2013). These disparities disappear mostly for higher multipoles, and therefore we exclude in general multipoles $l\leq 8$ when fitting spectral indices for the H i distribution.