# Thermal and turbulent properties of the Warm Neutral Medium in the solar neighborhood

Antoine Marchal, Marc-Antoine Miville-Deschenes

► To cite this version:

Antoine Marchal, Marc-Antoine Miville-Deschenes. Thermal and turbulent properties of the Warm Neutral Medium in the solar neighborhood. The Astrophysical Journal, 2021, 908 (2), 10.3847/1538-4357/abd108. [hal-03079687](https://hal.science/hal-03079687v1)

HAL Id: hal-03079687

https://hal.science/hal-03079687v1

Submitted on 1 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Thermal and Turbulent Properties of the Warm Neutral Medium in the Solar Neighborhood

Antoine Marchal<fnref n="1" /> and Marc-Antoine Miville-Deschênes<fnref n="2" />

Received 2020 November 10; revised 2020 December 3; accepted 2020 December 5; published 2021 February 23

## Abstract

The transition from the diffuse warm neutral medium (WNM) to the dense cold neutral medium (CNM) is what sets the initial conditions to the formation of molecular clouds. The properties of the turbulent cascade in the WNM, essential to describe this radiative condensation process, have remained elusive in part due to the difficulty in mapping out the structure and kinematics of each H I thermal phase. Here we present an analysis of a 21 cm hyper-spectral data cube from the GHIGLS H I survey where the contribution of the WNM is extracted using ROHSA, a Gaussian decomposition tool that includes spatial regularization. The distance and volume of the WNM emission is estimated using 3D dust extinction map information. The thermal and turbulent contributions to the Doppler line width of the WNM were disentangled using two techniques, one based on the statistical properties of the column density and centroid velocity fields, and the other on the relative motions of CNM structures as a probe of turbulent motions. We found that the volume of WNM sampled here (5.2 × 10⁵ pc³), located at the outer edge of the Local Bubble, shows thermal properties in accordance with expected values for heating and cooling processes typical of the solar neighborhood: $P_{th}/k_B = (4.4 \pm 2.6) \times 10^3$ K cm⁻³, $n = 0.74 \pm 0.41$ cm⁻³, and $T_k = (6.0 \pm 1.3) \times 10^3$ K. The WNM has the properties of sub/trans-sonic turbulence, with a turbulent Mach number at the largest scale probed here ($l = 130$ pc) of $M_s = 0.87 \pm 0.15$, a density contrast of $\sigma_{\rho/\rho_0} = 0.6 \pm 0.2$, and velocity and density power spectra compatible with $k^{-11/3}$. The low Mach number of the WNM provides dynamical conditions that allow the condensation mode of thermal instability to grow freely and form CNM structures, as predicted by theory.

*Unified Astronomy Thesaurus concepts: Solar neighborhood (1509); Warm neutral medium (1789); Interstellar atomic gas (833)*

## 1. Introduction

As for the universe as a whole, the hydrogen atom is the most abundant element in galaxies, being in neutral (H I), ionized (H II), or molecular (H₂) form. The neutral phase is of primary interest for understanding star formation in galaxies; most of the mass and the volume of the interstellar medium is dominated by the diffuse phases, and the H I to H₂ transition is a key step that leads to the formation of cold, dense and self-gravitating structures, where stars form. We know that this transition occurs in the cold neutral medium (CNM) ($T_k \sim 50$ K, $n \sim 50$ cm⁻³), but the physical processes that lead to the formation of these structures is still a matter of debate. One key element in understanding the formation of CNM clouds in galaxies seems to be related to the turbulent cascade acting in the precursor environment, the warm neutral medium (WNM) ($T_k \sim 8000$ K, $n \sim 0.3$ cm⁻³).

It is now well established from a theoretical perspective that the H I can be described as a multi-phase turbulent flow (e.g., Hennebelle & Audit 2007; Saury et al. 2014). In this picture, the CNM occupies only a few percent of the volume. It is the result of a condensation process occurring in the volume-filling WNM that behaves more like a classical compressible, isothermal turbulent flow. Numerical studies have shown that the outcome of the condensation process depends on the specific properties of turbulence in the WNM (e.g., Seifried et al. 2011; Saury et al. 2014; Bellomi et al. 2020). In fact, many have suggested that cold and dense interstellar clouds form by moderately supersonic compression of warm atomic gas streams (e.g., Hennebelle et al. 2008). These studies show

that strongly supersonic turbulence leads to dynamical times that are too short to leave enough time for the cooling and for the condensation process to be efficient. On the other hand, too weak turbulence does not provide enough density fluctuations and mixing to enable the formation of thermally stable cold structures. The other essential parameter is the pressure of the WNM; it has to be high enough to bring the gas into the thermally unstable range (Wolfire et al. 1995b, 2003).

Theory and numerical experiments have been of tremendous help in understanding this complex nonlinear process but this knowledge rests, for now, on a relatively small number of observations. In fact, many of the fundamental physical properties of the WNM remain elusive for now, like its exact kinetic temperature, volume-filling fraction, thermal pressure, and turbulent Mach number. Most of these quantities are deduced either by theory or indirectly by comparison with other phases. There is still an unresolved debate about how the WNM is distributed in interstellar space. Is it the extended skin of cold structures (McKee & Ostriker 1977) or a widespread inter-cloud medium? The question of its relation with ionized phases of the interstellar medium (ISM) (warm and hot) is still very much open (Cox 2005).

Early on, the ubiquity of a large and spectrally smooth 21 cm component, seen in emission but not in absorption, was reported. Dickey et al. (1979) mentioned that this warm, “not strongly absorbing,” gas corresponds to probably up to 75% of the column density at high Galactic latitudes. The fact that the CNM is subdominant in the mass budget of the H I was confirmed later: the CNM accounts for about 10%–30% of the mass in the diffuse ISM at high Galactic latitudes

<footnote n="1">¹Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 St. George Street, Toronto, ON M5S 3H8, Canada; amarchal@cita.utoronto.ca</footnote>

<footnote n="2">²AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, F-91191 Gif-sur-Yvette, France</footnote>
(Dickey et al. 1979; Heiles & Troland 2003; Haud & Kalberla 2007; Kalberla & Haud 2018; Murray et al. 2018, 2020) and possibly more in the vicinity of molecular clouds (40% according to Nguyen et al. 2019).

On the other hand Dickey et al. (1979) could not say how much of the WNM is warm ($T_k >> 1000$ K) and lukewarm (a few hundred K), a topic still debated. They had realized that very long integration times would be required to detect this gas in absorption and, to this day, the detection of the WNM in absorption remains a great challenge. So far, the highest temperatures estimated from absorption measurements are in the thermally unstable range (Mebold et al. 1982; Carilli et al. 1998; Dwarakanath et al. 2002; Kanekar et al. 2003; Begum et al. 2010; Murray et al. 2018; Nguyen et al. 2019) even with an optical depth sensitivity of $\sigma_{\tau} \sim 10^{-3}$ per channel. The highest temperature ($7200^{+1800}_{-1200}$ K) was obtained by Murray et al. (2014) who could detect a WNM signature by stacking 19 absorption spectra.

Under these conditions, to date the estimate of the WNM properties has been done using 21 cm emission profiles. Two techniques were used: one is based on absorption measurements against radio sources where the comparison between the absorption and emission profiles in the vicinity of the source provides a reliable description of the narrow CNM components. The contribution from the cold gas can then be removed from the emission profile, providing an estimate of the WNM column density and line width (e.g., Heiles & Troland 2003). The second method is to decompose emission spectra only, using a sum of Gaussians. That was done early on (e.g., Mebold 1972) and more recently using fully sampled data cubes (Haud & Kalberla 2007; Kalberla & Haud 2018). The average Doppler line width of the WNM at high Galactic latitudes found by these studies is 7.6–9.5 km s⁻¹. This method has the advantage of not being limited to lines of sight crossing radio sources and therefore can provide larger statistics of the WNM gas. Its main drawback is the difficulties related to the Gaussian decomposition (see Marchal et al. 2019, hereafter M19, for a discussion of the many pitfalls).

For both methods, the information extracted from the emission profiles is limited to the column density, velocity centroid, and Doppler line width. The difficulty to estimate the length on the line of sight over which the WNM emission is coming from, and the difficulty to separate thermal and turbulent contributions of the line, make the determination of the WNM density, filling factor, and turbulent Mach number currently elusive.

Nevertheless these quantities are essential to understanding the ISM evolution. In fact the process by which the diffuse gas condenses, increasing its density by several order of magnitudes, implying gathering matter over very large volumes, is central in the general evolution of galaxies (Cox 2005). The efficiency with which the CNM forms and the timescale over which these structure exist depend significantly on the dynamical properties of the inter-cloud medium, especially the pressure and turbulent Mach number of the WNM. This study is an attempt to extract more information from 21 cm emission data cubes, and bring new observational constraints on the dynamical properties of the warm gas in the solar neighborhood.

The paper is organized as follows. In Section 2, we present the data used in this work and the Gaussian decomposition performed to model its multiphase structure. In Section 3, we

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Column density map of the local-velocity component and intermediate-velocity component gas observed in NEP.</figcaption></figure>

analyze the thermal phases at local velocities, their spatial distribution along the line of sight, and volume-filling factors. In Section 4.3, we describe the methodology used to disentangle thermal and turbulent contributions from the Doppler line width of the local WNM. Gas properties (thermodynamic and turbulent) are analyzed in Section 5. Section 6 examines the static and dynamic scales of thermal instability (TI). A summary is provided in Section 7.

## 2. H I Spectral Data and Decomposition

### 2.1. Data

The 12° × 12° NEP field (or simply NEP) studied in this paper, located at $(l, b) = (96°:40, 30°:03)$ was part of the GHIGLS<fnref n="3" /> H I survey (Martin et al. 2015) obtained with the Auto-Correlation Spectrometer (ACS) at the Green Bank Telescope. NEP data have a channel spacing $\Delta v = 0.807$ km s⁻¹, an effective velocity resolution of about 1 km s⁻¹ and a spectral coverage $-365 < v$ [km s⁻¹] < 320. The spatial resolution of the ACS data was about 9/4. NEP was chosen from among other fields of the GHIGLS H I survey to avoid the effect of velocity crowding and self-absorption, important at low latitudes. This is also the field with the largest spatial coverage, allowing us to probe the statistical properties of turbulence in the WNM over almost two order of magnitude in spatial scales. Although there is a high-velocity component in NEP, it could be easily isolated in velocity. Only the intermediate-velocity component (IVC) and the local-velocity component (LVC) in the spectral range $-90 \le v$ [km s⁻¹] $\le 90$ were kept for the purpose of this work. Figure 1 shows the total column density map of NEP in this range.

### 2.2. Gaussian Decomposition

We performed a multiphase separation of NEP using the publicly available code **ROHSA**<fnref n="4" /> M19. **ROHSA** is a multi-

<footnote n="3">³ http://www.cita.utoronto.ca/GHIGLS/</footnote>

<footnote n="4">⁴ https://github.com/antoinemarchal/ROHSA</footnote>
<table><caption>**Table 1**
Mean Kinematic Properties (in km s⁻¹) of Gaussians Inferred from NEP Data</caption><thead><tr><th></th><th>G<sub>1</sub></th><th>G<sub>2</sub></th><th>G<sub>3</sub></th><th>G<sub>4</sub></th><th>G<sub>5</sub></th><th>G<sub>6</sub></th><th>G<sub>7</sub></th><th>G<sub>8</sub></th><th>G<sub>9</sub></th><th>G<sub>10</sub></th><th>G<sub>11</sub></th><th>G<sub>12</sub></th></tr></thead><tbody><tr><td>⟨μ<sub>n</sub>⟩</td><td>-74.1</td><td>-53.9</td><td>-44.7</td><td>-35.0</td><td>-22.9</td><td>-12.6</td><td>-4.8</td><td>-1.3</td><td>0.2</td><td>10.9</td><td>40.8</td><td>75.9</td></tr><tr><td>⟨σ<sub>n</sub>⟩</td><td>9.7</td><td>6.2</td><td>3.5</td><td>5.1</td><td>5.3</td><td>4.7</td><td>1.6</td><td>6.3</td><td>1.9</td><td>7.5</td><td>12.5</td><td>9.3</td></tr></tbody></table>

Gaussian-fitting code whose function is to decompose hyper-spectral observations into a sum of spatially coherent components. In addition, ROHSA was developed so that each resulting component has a similar velocity dispersion across the 2D position space, facilitating the association of each component to a given phase of the neutral ISM (WNM, LNM,<fnref n="5" /> or CNM).

The decomposition used here is that presented in Section 4 of M19. The model $\tilde{T}_b(v_z, \theta(r))$ used to fit the measured brightness temperature $T_b(v_z, r)$ at a projected velocity $v_z$ and coordinates $r$ is

$$ \tilde{T}_b(v_z, \theta(r)) = \sum_{n=1}^{N} G(v_z, \theta_n(r)), \quad (1) $$

with $\theta(r) = (\theta_1(r), \dots, \theta_n(r))$ and where

$$ G(v_z, \theta_n(r)) = a_n(r) \exp \left( - \frac{(v_z - \mu_n(r))^2}{2\sigma_n(r)^2} \right) \quad (2) $$

is a Gaussian parameterized by $\theta_n = (a_n, \mu_n, \sigma_n)$ with $a_n \ge 0$ being the amplitude, $\mu_n$ the position, and $\sigma_n$ the standard deviation 2D maps of the nth Gaussian profile across the plane of the sky. The parameters $\hat{\theta}$ are obtained by minimizing the cost function described in M19. It includes, for each parameter map $(a_n, \mu_n, \sigma_n)$, a Laplacian filtering that penalizes the small spatial frequencies whose strength is controlled by a hyper-parameter. An additional term, minimizing the variance of $\sigma_n$ is added to ensure the multiphase separation. ROHSA was initialized with a sum of $N=12$ Gaussians to ensure a complete encoding of the signal with spatially coherent components and each of the hyper-parameters $\lambda_a, \lambda_\mu, \lambda_\sigma$, and $\lambda'_\sigma$ has been set to 1000. We refer the reader to M19 for a detailed discussion about the choice of the five (including $N$) user parameters of ROHSA.

Mean velocity $\langle \mu_n \rangle$ and mean velocity dispersion $\langle \sigma_n \rangle$ of the 12 Gaussian components $G_n$ are tabulated in Table 1. We refer the reader to M19 for a visualization of their column density maps, velocity fields, and dispersion velocity fields.

### 3. H I Components in the NEP Field

#### 3.1. Local and IVCs

The Gaussian components extracted from NEP (see Table 1) can be separated in two groups; the first ($|\langle \mu_n \rangle| > 20$ km s⁻¹) is composed of intermediate velocity clouds (IVCs), the second ($|\langle \mu_n \rangle| < 20$ km s⁻¹) represents the LVC. Figure 2 shows column density maps of both components (top row) and the ratio of IVC column density fraction over total (bottom left). Many structures appear at all scales with significant variations of the column density over the field, from ~1 to $6 \times 10^{20}$ cm⁻². Interestingly the average and variations of the column density

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Top left: column density map of the local gas in NEP. Top right: column density map of intermediate-velocity components (IVCs). Bottom left: IVC mass fraction $f_{\text{IVC}}$ map. Bottom right: dust optical depth map from Planck Collaboration et al. (2014).</figcaption></figure>

are very similar for the LVC and IVC components. These variations are directly reflected in the $f_{\text{IVC}}$ map.

The dust optical depth at 353 GHz (i.e., 850 µm), $\tau_{353}$ (Planck Collaboration et al. 2014), also shown in Figure 2 (bottom right), has a structure very similar to the total (LVC + IVC) integrated H I emission (see Figure 1). The average column density deduced from the dust optical depth at 353 GHz is about $\langle N_H^{\tau_{353}} \rangle \simeq 6.0 \times 10^{20}$ cm⁻², assuming $\tau_{353}/N_H = 6.3 \times 10^{-27}$ of Planck Collaboration et al. (2014). This value is identical to the value deduced from the total 21 cm integrated emission (sum of the IVC and LVC): $\langle N_{\text{H I}}^{\text{tot}} \rangle = 6.0 \times 10^{20}$ cm⁻². This implies that the IVCs, which constitute about 40% of the H I column density in NEP, contain dust with a similar dust-to-gas ratio than the local gas, as found by Planck Collaboration et al. (2011).

#### 3.2. Thermal Phases at Local Velocities

For the purpose of our analysis dedicated to the study of the WNM, we went a step further in the analysis of the result of the Gaussian decomposition and produced a multi-phase view of the LVC component. To do so we combined the Gaussian components at LVC velocities into three groups of similar width: WNM, LNM, and CNM. Specifically, Gaussian $G_8$ and

<footnote n="5">⁵ LNM stands for lukewarm neutral medium.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Top: integrated column density $N_{HI}^{WNM}$, $N_{HI}^{LNM}$ and $N_{HI}^{CNM}$ fields of the three-phase model. Bottom: $f_{WMN}$, $f_{LNM}$ and $f_{CNM}$ column density fraction maps. Letters A, B, and C indicate the three anti-correlated regions of warm neutral medium and lukewarm neutral medium discussed in the text.</figcaption></figure>

$G_{10}$ are added to obtain the WNM, Gaussian $G_6$ forms the LNM, and Gaussian $G_7$ and $G_9$ are combined to form the CNM (see Table 1).

Integrated column density fields of each phases are shown in Figure 3 (top row) as well as the mass fraction with respect to the total LVC column density (bottom row) The column density maps of the WNM and LNM are anti-correlated in some areas, three in particular that we have annotated A, B, and C (see Figure 3, top left). At each of these three positions, the low WNM column density areas seem to correspond to elongated structures in the LNM. These features appear clearly in the mass fraction maps. On the other hand, we note that the column density map of the CNM, composed of clumpy and filamentary structures, does not appear to be reflected in any other phases, except maybe for region C. Interestingly, in this region the CNM structures seem to have a smaller spatial extension than LNM structures.

The one-point probability distribution functions (PDFs) of the column density and mass fraction, for the three phases, are shown in Figures 4 and 5, respectively. Average and standard deviation values of each quantity are tabulated in Table 2. Each

phase shows a rather complex column density distribution that is not particularly well described by a log-normal. The median column density of the whole LVC is $2.25 \times 10^{20} \text{ cm}^{-2}$. It is dominated by the WNM component, which is 2.5 times larger than the LNM and 10 times larger than the CNM. On average, 64% of the mass is in the WNM, 28% is in the LNM, and only 8% is in the CNM. However, as Table 2 shows, the standard deviation of these maps is relatively high (35%, 53%, and 135% of the median value for the WNM, LNM, and CNM, respectively). This reflects the important variations seen in Figure 3 (bottom row) and Figure 5, especially strong for the CNM which is very intermittent spatially, with a significant fraction of the field (55%) with a mass fraction below 1%.

Such a spatial variability of the column density and mass fraction in each phase is to be expected because of the dynamical nature of TI. The values found in NEP are just one instance of the possible distributions. Interestingly, the rather larger variability of the mass fractions observed here compares well with values that are found all over the sky with absorption surveys. As shown in Figures 4 and 5, for each phase, the ranges of $N_{HI}$ and $f$ found in NEP are included in the broader
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Probability distribution function (PDF) of the column density of cold neutral medium (CNM) (blue), lukewarm neutral medium (LNM) (green) and warm neutral medium (WNM) (red). PDFs of WNM, LNM and CNM inferred from the 21-Sponge survey (Murray et al. 2018) are shown in magenta, lime and cyan, respectively.</figcaption></figure>

distributions deduced from the 21-Sponge survey data by Murray et al. (2018). On the other hand, we note that some $N_{HI}$ values reported by Murray et al. are significantly larger than the maximum values found in NEP, but most of them are found at lower Galactic latitudes than NEP, where the lines of sight are longer. Therefore, because of the different line of sight lengths sampled in the two studies, a comparison of the column density statistics has a limited value. On the other hand, the mass fraction of the different phases is something that should be less dependent on the length on the line of sight.<fnref n="6" /> We found that the CNM mass fraction does not exceed 0.35 in NEP while it reaches more than 0.8 in 21-Sponge. We note that 61% of the $f_{CNM}$ 21-Sponge values larger than 0.35 are at lower Galactic latitudes than in NEP but in this case the latitude cannot explain this discrepancy. Local effects like the Local Bubble where the CNM fraction is likely to be smaller than the Galactic average could explain the difference seen here.

A comparison with the results of Murray et al. (2018) highlights the fact that NEP has a low CNM fraction which makes it a particularly well suited field to study the properties of the diffuse inter-cloud medium.

### 3.3. Spatial Distribution of Gas along the Line of Sight

In order to make a quantitative study of the physical properties of the local H I in NEP, one needs to estimate the distance and depth of the emitting gas. The decomposition of the 21 cm data with ROHSA enables a deblending of the LVC and IVC components, as well as a separation of the different thermal phases of the LVC, but for such high Galactic latitude fields, the velocity information provided by the 21 cm data cannot be used to estimate a kinematic distance. In order to put

<figure><img src="image_7.png" /><figcaption>**Figure 5.** As Figure 4 but for the mass fraction.</figcaption></figure>

<table><caption>**Table 2**
Mean Properties of Column Density Map for the Local Velocity Component</caption><thead><tr><th></th><th>⟨N<sub>H I</sub>⟩<br>10<sup>18</sup> cm<sup>-2</sup></th><th>σ<sub>NH I</sub><br>10<sup>18</sup> cm<sup>-2</sup></th><th>σ<sub>NH I</sub>/⟨N<sub>H I</sub>⟩</th><th>⟨f⟩</th><th>σ<sub>f</sub></th></tr></thead><tbody><tr><td>WNM</td><td>146.5</td><td>51.3</td><td>0.35</td><td>0.64</td><td>0.13</td></tr><tr><td>LNM</td><td>56.9</td><td>30.1</td><td>0.53</td><td>0.28</td><td>0.11</td></tr><tr><td>CNM</td><td>13.5</td><td>18.2</td><td>1.35</td><td>0.08</td><td>0.06</td></tr></tbody></table>

constraints on the distance and depth of the gas, we rely on 3D tomography.

For the last two decades, and even more so now that Gaia data are available, significant effort has been put into mapping the ISM in 3D using tomography techniques (Lallement et al. 1995, 2003, 2018, 2019; Vergely et al. 2010; Green et al. 2014, 2018; Rezaei et al. 2017). Of these studies, the only one that mapped the local atomic ISM directly is Vergely et al. (2001) who applied a regularized Bayesian inversion method on interstellar Lyα absorption measurements toward 454 stars. In that way they could produce a 60 pc resolution map of the 3D density distribution of the neutral hydrogen within 250 pc of the Sun. Even with this coarse spatial resolution, their 3D map gives us an idea of how the H I is distributed in the direction of NEP. A large cavity (the Local Bubble) is observed around the Sun with a volume density of about 0.1 cm⁻³. This cavity seems to extend over 200 pc, the limit reached by this study.

No direct 3D mapping of H I has been performed beyond 250 pc. However, it is possible to approach it using a 3D dust extinction map. Using the same formalism as Vergely et al. (2001), Lallement et al. (2019) used Gaia DR2 photometric data combined with 2MASS to derive extinction measurements toward stars in a 6 × 6 × 0.8 kpc³ volume around the Sun. Here we use this product to compute the mean dust extinction per parsec profile, $a_v$, in the direction of NEP (see Figure 6, top). The presence of the Local Bubble seen in Vergely et al. (2001), appears clearly in this profile. From 150 to 250 pc, there is a

<footnote n="6">⁶ $f_{CNM}$ does not vary with latitude in a plane-parallel model with two different scale heights for the CNM and the WNM.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Top: mean extinction profile $a_v$/pc in NEP based on the 3D extinction dust map performed by Lallement et al. (2019). Bottom: mean H I density profile in NEP for an $N_H/E(B-V)$ ratio between 6 and $10 \times 10^{21}$ (gray), and total H I local density profile from Dickey & Lockman (1990) (blue). Vertical dashed lines indicate the size of the plateau, used as the typical depth $L_z$ of the fluid along the line of sight.</figcaption></figure>

smooth increase of the dust extinction per parsec, also seen in the $N_{H\text{I}}$ distribution of Vergely et al. (see their Figure 6, top left). This makes it difficult to estimate a clear frontier for the Local Bubble. Beyond this first wall, $a_v$ remains rather constant for ~130 pc and then starts to decrease smoothly, on a 100 pc scale.

From the $a_v$ profile, one can infer a gas volume density using the $N_H/E(B-v)$ ratio:

$$n_{H}(l) = \frac{a_{v}(l)}{R_{v}} \times \frac{N_{H}}{E(B - V)}. \qquad (3)$$

Recently, the classical value of $N_H/E(B-v) = 5.8 \times 10^{21}$ cm$^{-2}$ mag$^{-1}$ estimated by Bohlin et al. (1978) has been the subject of significant revisions. Using different techniques, Planck Collaboration et al. (2014), Liszt (2013), Lenz et al. (2017), Zhu et al. (2017), and Nguyen et al. (2018) all obtained significantly higher values in the range, $8 \pm 1 \times 10^{21}$ cm$^{-2}$ mag$^{-1}$.

To visualize a realistic range of possibilities, we present in Figure 6 (bottom) the $n_H(l)$ density profile deduced from dust extinction, assuming $R_v = 3.1$ and considering a $N_H/E(B-v)$ ratio between 6 and $10 \times 10^{21}$ cm$^{-2}$ mag$^{-1}$. For the sake of comparison, we over-plot the H I density profile from Dickey & Lockman (1990) which describes the average H I volume density profile, $n(z)$, in the inner Galaxy, assuming cylindrical symmetry. This average H I model and the model computed from the dust extinction curve appear to be roughly consistent in the range $l \sim 200-400$ pc where most of the local H I mass is. This agreement is reassuring but one should not expect a perfect match as the Dickey & Lockman (1990) model does not take into account the particular conditions present in the solar neighborhood, i.e., the Local Bubble.

At this point it is interesting to compare the column density computed by integrating the dust extinction curve

$$N_{H}^{A_{v}} = \int n_{H}(l) dl, \qquad (4)$$

to the column density inferred from the 21 cm line emission. Integrating Equation (4) using the profile $a_v$ shown in Figure 6,

and assuming $N_H/E(B-v) = 8 \times 10^{21}$ cm$^{-2}$ mag$^{-1}$, we obtain $N_H^{A_v} = 2.4 \times 10^{20}$ cm$^{-2}$. This is significantly lower than the average column density estimated from the dust optical depth or the total H I emission: $N_H^{\tau_{353}} = N_{H1}^{\text{tot}} = 6.0 \times 10^{20}$ cm$^{-2}$. Because they are related to dust, both $N_H^{\tau_{353}}$ and $N_H^{A_v}$ trace the total amount of hydrogen atoms along the line of sight (H I, H II, and H$_2$). The fundamental difference between these two estimates is that $N_H^{\tau_{353}}$ traces dust emission to infinity while $N_H^{A_v}$ stops at $l \sim 1000$ pc, limited by the Gaia sample. This difference could be explained if another dusty cloud were present beyond the detection limit of Lallement et al. (2019).

Interestingly it appears that the median integrated column density of the local gas (LVC) inferred from the 21 cm line is $N_{H1}^{\text{LVC}} = 2.3 \times 10^{20}$ cm$^{-2}$, a value in very good agreement with that deduced from the dust profile $a_v$ (relative difference of 7.4 %). Note that even if we consider the lower and upper value of $N_H/E(B-v)$ ($6-10 \times 10^{21}$ cm$^{-2}$ mag$^{-1}$), the relative difference between $N_H^{A_v}$ and $N_H^{\text{LVC}}$ remains within about $\pm 30\%$. This suggests that the IVC in NEP is located beyond 1 kpc from the Sun (i.e., in the Galactic halo, at $z > 575$ pc).

Building on this, in the following we will assume that the $a_v$ profile provides a description of the distribution of the LVC component along the line of sight. Then it becomes possible to estimate the physical sizes $L_x$ (longitude), $L_y$ (latitude), and $L_z$ (depth) of the LVC in the region of NEP. Using the size of the plateau delimited by two vertical dashed lines in Figure 6, we approximate the typical depth of the H I along the line of sight, $L_z \sim 130$ pc. At $l = 300$ pc, halfway down $L_z$, the $12^\circ \times 12^\circ$ region on the plane of the sky translates into $L_x = L_y \sim 63$ pc.

### 3.4. Volume-filling Factors

Based on the previous 3D tomography analysis, we will now assume that the local velocity H I component is coming from a slab 130 pc deep along the line of sight. The question now is: how are the thermal phases distributed in this volume $V$? In other words, what is the volume-filling factor of each phase?

To explore this it is useful to turn to numerical simulations, like that of Saury et al. (2014). Figure 7 shows the two-dimensional PDF of the mass fraction and the volume-filling factor for the WNM (red), LNM (green), and CNM (blue) (see top and right distributions for their respective one-point PDFs) for a $10 \times 10 \times 40$ pc$^3$ cube, the same one used in M19 (see their Section 3). This specific simulation represents one realization of a thermally bistable H I fluid. The mass fraction in each phase depends on the initial conditions and on the properties of the turbulent forcing. It appears that this one reproduces rather well the mass fractions observed in NEP (see Figures 5 and 7, top).

Because the CNM has a density about 100 times larger than the WNM, it occupies a very small part of the volume. In the case of NEP where the CNM mass fraction is low, $f_M^{\text{CNM}} = 0.08$, the CNM volume-filling factor is likely to be $f_V^{\text{CNM}} \sim 0.001$. Therefore, most of the volume is filled by the low-density LNM and WNM thermal components, a property shared with the numerical simulation used here.

From the numerical simulation, we note that the LNM and WNM mass fractions and volume-filling factors are correlated (see Figure 7); on average, a larger mass fraction implies a larger volume-filling factor. We also note that the structures of the LNM and WNM in the $f_M - f_V$ plane are anti-correlated.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Two-dimensional probability distribution function of the volume-filling factor as function of the mass fraction for the WNM (red), LNM (green), and CNM (blue) along each line of light of a 10 × 10 × 40 pc³ region with moderate CNM fraction (same region as in M19 of the numerical simulation of thermally bistable turbulence performed by Saury et al. 2014).</figcaption></figure>

This reflects the fact that these two thermal phases are dynamically linked and that their average densities are similar. This link between the LNM and WNM translates into an anti-correlation in the column density as well (Saury et al. 2014), something that is seen in NEP (see regions A, B, and C in Figure 3). This opens the possibility to use the mass fraction of the LNM and WNM as an indicator for their respective volume-filling factors. From this, and neglecting the volume occupied by the CNM, we conclude that the LNM and WNM fill about 30% and 70% of the volume of H I respectively (see Table 2).

This coarse estimate of $f_y$ was obtained in the idealized context where the totality of the gas is neutral. In fact this is an upper limit as some fraction of the volume in the direction of NEP is likely to be occupied by the diffuse warm ionized medium (WIM). The joint analysis of pulsar dispersion measures (DMs) and diffuse Hα emission by Gaensler et al. (2008) (see also Reynolds 1977; Berkhuijsen et al. 2006) revealed that the WIM is best described by a collection of ionized structures (as opposed to a continuous medium) filling about 10% of the volume (increasing with Galactic altitude) and with a global scale height of $1830_{-250}^{+120}$ pc. In this context it is difficult to evaluate precisely the fraction of the volume occupied by the WIM in the 130 pc slab we are considering here. In addition we note that the correspondence between $N_H^{tot}$ and $N_H^{tot}$ suggests that the column density of ionized gas $N_H$ is low in the direction of NEP, or that the WIM is deficient in dust (the distinction between these two possibilities remains impossible to date). In any case one should keep in mind that 10%–15% of the volume considered here could be filled with ionized gas.

# 4. The WNM: Disentangling Thermal and Turbulent Observable Properties

The decomposition of the 21 cm data and the grouping of the widest components at local velocities has allowed us to build a model of the data cube of the WNM emission originating from the Galactic disk in the direction of NEP. This contains information not only on the column density of the WNM but also on the gas temperature and on its velocity component along the line of sight. In what follows, building on the previous constraints on the depth of local H I we expose how we analyze the modeled cube to extract information about the thermal and turbulent observable properties of the warm H I phase.

## 4.1. Observable

In emission, the observed 21 cm brightness temperature at position $\mathbf{r}$ and velocity $v_z$, $T_b(v_z, \mathbf{r})$, depends on the variations along the line of sight of the gas volume density, $\rho$, the component of the velocity field along the los, $v_z$, and on the kinetic temperature $T_k$ (Miville-Deschênes & Martin 2007). In the optically thin approximation the dependence of $T_b(v_z, \mathbf{r})$ on these physical quantities can be expressed as

$$T_b(v_z, \mathbf{r}) \propto \int_0^{L_z} \rho(\mathbf{r}, z) \exp\left(-\frac{(v_z' - v_z(\mathbf{r}, z))^2}{2\sigma_{\text{th}}^2(\mathbf{r}, z)}\right) dz, \quad (5)$$

where the integral is over the depth $L_z$ of the medium. Note that the thermal velocity dispersion, $\sigma_{\text{th}}$, also varies with position $z$ along the line of sight as it depends on the kinetic temperature field, $T_k(\mathbf{r}, z)$.

From this, the zeroth, first, and second moments of the emission, the column density, $N_{H1}(\mathbf{r})$, the velocity centroid, $C(\mathbf{r})$, and the velocity dispersion, $\sigma_{T_b}(\mathbf{r})$, can be expressed the following way:

$$N_{H1}(\mathbf{r}) = \int_{0}^{L_z} \rho(\mathbf{r}, z) dz, \quad (6)$$

$$C(\mathbf{r}) = \frac{1}{N_{H1}(\mathbf{r})} \int_{0}^{L_z} \rho(\mathbf{r}, z) v_z(\mathbf{r}, z) dz, \quad (7)$$

and

$$\sigma_{T_b}^2(\mathbf{r}) = \frac{1}{N_{H1}(\mathbf{r})} \int_0^{L_z} \rho(\mathbf{r}, z) v_z^2(\mathbf{r}, z) dz - C^2(\mathbf{r}) + \frac{k_B \overline{T_k}(\mathbf{r})}{m_H}, \quad (8)$$

where $k_B$ is the Boltzmann constant, $m_H$ is the hydrogen atom mass, and $\overline{T_k}(\mathbf{r})$ is the average kinetic temperature along the line of sight, weighted by density:

$$\overline{T}_k(\mathbf{r}) = \frac{1}{N_{H1}(\mathbf{r})} \int_0^{L_z} \rho(\mathbf{r}, z) T_k(\mathbf{r}, z) dz. \quad (9)$$

## 4.2. Density Contrast

First, we analyze the density contrast of the WNM in NEP. To do so we used the formalism developed by Brunt et al. (2010) to reconstruct the variance of a 3D physical field from 2D observation assuming only statistical isotropy of the underlying 3D physical field and knowing the physical scales
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Left: density contrast as a function of volume-filling factor $f_V^{WMN}$ (or effective depth $L_z' = f_V^{WMN}$ $L_z$ on the top axis). Green and orange areas show the volume occupied by the LNM and warm ionized medium, respectively. Right: amplitude of turbulent motions as a function of $f_V^{WMN}$. The horizontal blue line shows the turbulent velocity dispersion obtained from the cloud–cloud velocity dispersion of CNM structures. The red dashed line and red area show our final estimate of $f_V^{WMN} = 0.5 \pm 0.2$.</figcaption></figure>

(in 3D) of the fluid. This formalism is general enough to be applied to any kind of field, including density and velocity. First, we use this formalism to infer the density contrast $\sigma_{\rho/\rho_0}$ of the 3D density field from the column density contrast of its projection along the line of sight $\sigma_{N/N_0} = 0.346$. Using Parseval’s theorem, Brunt et al. (2010) showed that $\sigma_{N/N_0}$ and $\sigma_{\rho/\rho_0}$ are linked by the ratio

$$R = \left( \frac{\sigma_{N/N_0}}{\sigma_{\rho/\rho_0}} \right)^2 \quad (10)$$

$$= \frac{\left( \sum_{k_x = -L_x/2+1}^{L_x/2} \sum_{k_y = -L_y/2+1}^{L_y/2} P_{\rho}^{3D}(k) - P_{\rho}^{3D}(0) \right)}{\left( \sum_{k_x = -L_x/2+1}^{L_x/2} \sum_{k_y = -L_y/2+1}^{L_y/2} \sum_{k_z = -L_z/2+1}^{L_z/2} P_{\rho}^{3D}(k) - P_{\rho}^{3D}(0) \right)} \quad (11)$$

where $P_{\rho}^{3D}(k)$ is the power spectrum of $\rho$. For a given field of size $L_x \times L_y$, two parameters control the ratio $R$: the slope of $P_{\rho}^{3D}(k)$, and the depth of the cube, $L_z$, over which velocity fluctuations are averaged.

An important point has to be made here. Because of the fact that the WNM does not fill the volume fully, the scale over which the density fluctuations are averaged is not $L_z = 130$ pc but it is over an effective depth $L_z' = f_V^{WMN} \times L_z$. It is natural to expect, also due to $f_V^{WMN}$, that the power spectrum of the column density (projected 2D field), usually taken as a direct measure of the 3D density power spectrum, will in fact differ from $P_{\rho}^{3D}(k) \propto k^{-11/3}$ (see Section 5.2.1 and Appendix B). To account for this effect, we assumed $P_{\rho}^{3D}(k) \propto k^{-11/3}$, representative of a sub/trans-sonic turbulence.

Assuming $P_{\rho}^{3D}(k) \propto k^{-11/3}$, $L_x = L_y = 63$ pc and $L_z = 130$ pc, the density contrast depends on the WNM volume-filling factor. Figure 8 (left) shows the recovered value of $\sigma_{\rho/\rho_0}$ for the range $0 < f_V^{WMN} < 1$. Based on this analysis, the maximum possible value is $\sigma_{\rho/\rho_0} = 1.14$, but given our previous estimate of the

WNM volume-filling factor, $f_V^{WMN} \lesssim 0.6$, the contrast is likely to be $\sigma_{\rho/\rho_0} \lesssim 0.72$. Note that this estimate will be refined using further constrains (see Section 4.3.2) and $\sigma_{\rho/\rho_0} = 0.6 \pm 0.2$ will be adopted in the following.

### 4.3. Thermal and Turbulent Velocity Dispersions

The width of the emission line is given by Equation (8). In this equation, the first two terms on the right do not depend on $T_k$; they represent the second-order moment of the z-component of the velocity field, weighted by density. The third term depends only on $T_k$; it represents the thermal broadening of the line. Equation (8) is often written in a more compact form:

$$\sigma_{T_b}^2 = \sigma_{v_z}^2 + \sigma_{\text{th}}^2 \quad (12)$$

The Doppler line width $\sigma_{T_b}$ represents the width of the 21 cm emission line. Equations (5)–(8) show that it is a quadratic sum of the thermal and turbulent motions of the gas along the line of sight, both weighted by density. This observational mixture makes their separation impossible when considering the emission spectrum of a single line of sight. Therefore, a significant challenge in trying to extract the physical properties of turbulence of the WNM, for example its turbulent Mach number, is to estimate the thermal and turbulent contributions to the Doppler line width. In order to do that, we investigated two independent methods. The first uses the fluctuations of the centroid velocity field of the WNM $C(r)$ and the second uses the CNM components extracted with ROHSA as tracer particles of the WNM velocity field. For each method, we make the assumption that the statistical properties of the velocity field are isotropic (i.e., $v_z$ is statistically representative of the three directions).

#### 4.3.1. Use of the Centroid Velocity Field of the WNM

The $C(r)$ map of the WNM in NEP is shown in Figure 9. To our knowledge it is the first time that the centroid velocity of
<figure><img src="image_11.png" /><figcaption>Figure 9. Centroid velocity field of the WNM at local velocities in NEP.</figcaption></figure>

the WNM has been visualized. This map shows shallow fluctuations ($\sigma_{v_z} = 2.06$ km s⁻¹). As shown previously (see Equation (7)) the velocity centroid, $C(r)$, involves a complex combination of density and velocity fluctuations. Because of this mixture, using $C(r)$ to infer the statistical properties of $v_z(r, z)$ can only be done under very specific conditions. As shown by several studies (e.g., Miville-Deschênes et al. 2003; Levrier 2004) the statistics of the 3D velocity field, $v_z(r, z)$, can be recovered from $C(r)$ when density fluctuations are small compared to the mean density of the fluid ($\sigma_{\rho/\rho_0} \lesssim 0.5$; see Ossenkopf et al. 2006). In Section 4.2 we found $\sigma_{\rho/\rho_0} \lesssim 0.7$. Therefore, we think that, in the very specific case of the WNM velocity field under study here, it is a fair approximation to use $C(r)$ as probe of $v_z(r, z)$.

As stated in Section 4.2, Equation (11) can also be applied on the centroid velocity field to infer the velocity dispersion of the 3D velocity field $\sigma_{v_z}$ from the velocity dispersion of its projection along the line of sight, $\sigma_{\langle v_z \rangle} = 2.06$ km s⁻¹. Figure 8 (right) shows the recovered value of $\sigma_{v_z}$ for the range $0 < f_V^{\text{WNM}} < 1$. The maximum possible value is $\sigma_{v_z} = 6.8$ km s⁻¹, but given our previous estimate of the WNM volume-filling factor, $f_V^{\text{WNM}} \lesssim 0.6$, the turbulent contribution to the line width is likely to be $\sigma_{v_z} \lesssim 4.3$ km s⁻¹.

### 4.3.2. Use of CNM Structures as Tracer Particles

To constrain further the amplitude of turbulent motions in the WNM we propose to use the cloud–cloud velocity dispersions $\sigma_{v_z}$ of CNM structures as being representative of the turbulent velocity field $\sigma_{v_z}$ of the inter-cloud medium. In this picture, the volume occupied by each dense CNM structure is very small compared to the volume in which they move. By analogy to fluid mechanics laboratory experiments, these CNM structures can be seen as tracer particles added to the fluid (WNM) to highlight the statistical properties of these turbulent motions. This method makes the assumption that the CNM structures observed in NEP are homogeneously distributed in the $V$.

To compute the velocity dispersion of CNM structures, we use the ensemble of Gaussians from components $G_7$ and $G_9$

(see Table 1). We remind the reader that each Gaussian $i$ has three parameters, namely its amplitude, $a_i$, its projected velocity along the line of sight, $\mu_i$, and its velocity dispersion, $\sigma_{v_z}$. The velocity dispersion of CNM structures, $\sigma_{v_z}$, is defined as the square of the weighted variance of the velocity parameters of the Gaussian sample:

$$ \sigma_{v_z, c} = \sqrt{\frac{\sum_{i=0}^{N} w_i (\mu_i - \bar{\mu}^*)^2}{\sum_{i=0}^{N} w_i}}, \quad (13) $$

where

$$ \bar{\mu}^* = \frac{\sum_{i=0}^{N} w_i \mu_i}{\sum_{i=0}^{N} w_i}, \quad (14) $$

and

$$ w_i = a_i \sigma_i. \quad (15) $$

We find a velocity dispersion $\sigma_{v_z} = 3.5$ km s⁻¹. Going back to the method presented in the previous section, if this value is representative of the WNM turbulent motion along the line of sight, it implies $f_V^{\text{WNM}} = 0.43$ (see the blue dotted line in Figure 8 (right)).

Combining constraints from the centroid velocity field of the WNM and the cloud–cloud velocity dispersions of CNM structures, we adopted values that are the average of the two methods, and uncertainties based on the limits set by them. Therefore we conclude that $\sigma_{v_z} = 3.8 \pm 0.4$ km s⁻¹, and $f_V^{\text{WNM}} = 0.5 \pm 0.2$. Using this better estimate of $f_V^{\text{WNM}}$, we refine our estimate of the density contrast $\sigma_{\rho/\rho_0} = 0.6 \pm 0.2$ (see Figure 8 (left)).

## 5. Gas Properties of the WNM in NEP

A number of physical, thermodynamic, and turbulent properties of the WNM in NEP, as calculated in separate subsections below, are summarized in Table 3.

### 5.1. Thermodynamic Properties

The kinetic temperature, volume density, and thermal pressure are

$$ T_k = \frac{\sigma_{\text{th}}^2 m_H}{k_B}, \quad (16) $$

$$ n = \frac{\langle N_H \rangle}{L_z'}, \quad (17) $$

and

$$ P_{\text{th}}/k_B = nT_k, \quad (18) $$

respectively. Subtracting $\sigma_{v_z} = 3.8$ km s⁻¹ from the observed mean velocity dispersion (see Equation (12)), the mean thermal velocity dispersion is $\sigma_{\text{th}} = 7.0 \pm 0.7$ km s⁻¹. Using Equations (16)–(18), we find that $T_k = (6.0 \pm 1.3) \times 10^3$ K, $n = 0.74 \pm 0.41$ cm⁻³, and $P_{\text{th}}/k_B = (4.4 \pm 2.6) \times 10^3$ K cm⁻³.

Figure 10 shows a $P_{\text{th}}/k_B-n$ diagram where the pressure and density inferred in this work are shown in red. For the sake of comparison, this figure also shows the range of pressure measured in the CNM at high Galactic latitudes by Jenkins & Tripp (2011) (blue line), as well as the standard equilibrium curves of Wolfire et al. (2003) computed for $N_c = 1 \times 10^{20}$ cm⁻²
<table><caption>**Table 3**
Mean Properties of the Local Emission of the WNM in NEP</caption><thead><tr><th>Quantity</th><th>Symbol</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td colspan="4"><strong>Observable</strong></td></tr><tr><td>Column density</td><td>⟨N<sub>H</sub>⟩</td><td>1.47 × 10<sup>20</sup></td><td>cm<sup>-2</sup></td></tr><tr><td>Largest scale observed</td><td>L<sub>z</sub></td><td>130 ± 20</td><td>pc</td></tr><tr><td>Doppler velocity dispersion at L<sub>z</sub></td><td>σ<sub>Tb</sub></td><td>8.0 ± 0.6</td><td>km s<sup>-1</sup></td></tr><tr><td>Volume-filling factor</td><td>f<sub>V</sub><sup>WNM</sup></td><td>0.5 ± 0.2</td><td></td></tr><tr><td colspan="4"><strong>Thermodynamic</strong></td></tr><tr><td>Thermal velocity dispersion</td><td>σ<sub>th</sub></td><td>7.0 ± 0.7</td><td>km s<sup>-1</sup></td></tr><tr><td>Kinetic temperature</td><td>T<sub>k</sub></td><td>(6.0 ± 1.3) × 10<sup>3</sup></td><td>K</td></tr><tr><td>Volume density</td><td>n</td><td>0.74 ± 0.41</td><td>cm<sup>-3</sup></td></tr><tr><td>Thermal pressure</td><td>P<sub>th</sub>/k<sub>B</sub></td><td>(4.4 ± 2.6) × 10<sup>3</sup></td><td>K cm<sup>-3</sup></td></tr><tr><td colspan="4"><strong>Turbulent cascade</strong></td></tr><tr><td>Density contrast at L<sub>z</sub></td><td>σ<sub>p/ρ<sub>0</sub></sub>(L<sub>z</sub>)</td><td>0.6 ± 0.2</td><td></td></tr><tr><td>Turbulent velocity dispersion at L<sub>z</sub></td><td>σ<sub>vz</sub>(L<sub>z</sub>)</td><td>3.8 ± 0.4</td><td>km s<sup>-1</sup></td></tr><tr><td>Turbulent sonic Mach number at L<sub>z</sub></td><td>M<sub>s</sub>(L<sub>z</sub>)</td><td>0.87 ± 0.15</td><td></td></tr><tr><td>Sound speed</td><td>C<sub>s</sub></td><td>7.67 ± 0.82</td><td>km s<sup>-1</sup></td></tr><tr><td>Turbulent velocity dispersion at 1 pc</td><td>σ<sub>vz</sub>(1)</td><td>0.76 ± 0.1</td><td>km s<sup>-1</sup></td></tr><tr><td>Turbulent sonic Mach number at 1 pc</td><td>M<sub>s</sub>(1)</td><td>0.17 ± 0.03</td><td></td></tr><tr><td>Mean free path</td><td>λ</td><td>(4.4 ± 2.4) × 10<sup>-4</sup></td><td>pc</td></tr><tr><td>Kinematic molecular viscosity</td><td>ν</td><td>(5.5 ± 3.1) × 10<sup>20</sup></td><td>cm<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Knudsen number at L<sub>z</sub></td><td>K<sub>n</sub></td><td>(3.4 ± 1.8) × 10<sup>-6</sup></td><td></td></tr><tr><td>Reynolds number at L<sub>z</sub></td><td>Re</td><td>(3.3 ± 1.8) × 10<sup>5</sup></td><td></td></tr><tr><td>Sonic scale</td><td>λ<sub>s</sub></td><td>(1.9 ± 1.1) × 10<sup>2</sup></td><td>pc</td></tr><tr><td>Relative strength of solenoidal to compressive modes</td><td>ζ</td><td>0.5 ± 0.4</td><td></td></tr><tr><td>Dissipation scale</td><td>η</td><td>(1.0 ± 0.4) × 10<sup>-2</sup></td><td>pc</td></tr><tr><td>Dissipation time</td><td>t<sub>η</sub></td><td>50 ± 18</td><td>kyr</td></tr><tr><td>Energy transfer rate</td><td>ε</td><td>(1.2 ± 0.4) × 10<sup>-4</sup></td><td>L<sub>☉</sub> M<sub>☉</sub><sup>-1</sup></td></tr><tr><td>Ambipolar diffusion scale</td><td>l<sub>AD</sub></td><td>(3.2 ± 0.9) × 10<sup>-2</sup></td><td>pc</td></tr><tr><td colspan="4"><strong>Thermal instability</strong></td></tr><tr><td>Cooling time</td><td>t<sub>cool</sub></td><td>≈2.3</td><td>Myr</td></tr><tr><td>Dynamical time at sonic scale</td><td>t<sub>dyn</sub>(λ<sub>s</sub>)</td><td>25 ± 12</td><td>Myr</td></tr><tr><td>Cooling length</td><td>λ<sub>cool</sub></td><td>≈17.5</td><td>pc</td></tr><tr><td>Condensation criterion at sonic scale</td><td>I(λ<sub>s</sub>)</td><td>0.09 ± 0.04</td><td></td></tr><tr><td>Field length</td><td>λ<sub>F</sub></td><td>0.10 ± 0.06</td><td>pc</td></tr></tbody></table>

and three different values of FUV interstellar radiation field strength $G_0$ = (0.84, 1.7, 3.4) (in units of the Draine 1978 field strength), representative of expected variations in the diffuse ISM.<fnref n="7" />

We note that the WNM pressure estimated here is in complete agreement with the average value found in the local CNM by Jenkins & Tripp (2011) ((3.8 ± 1.5) × 10³ K cm⁻³). In addition, considering the realistic range of $G_0$ shown in Figure 10, our results show that the values of pressure and density are compatible with a WNM in thermal equilibrium calculated with standard, solar neighborhood properties for the cooling and heating mechanisms.

<figure><img src="image_12.png" /><figcaption>**Figure 10.** $P_{th}/k_B-n$ diagram: the red cross shows the values obtained in the WNM of NEP. The standard model ($N_c = 1 \times 10^{20} \text{ cm}^{-2}$) of Wolfire et al. (2003) is over-plotted for three different values of the FUV interstellar radiation field strength $G_0$ = (0.84, 1.7, 3.4). The blue line shows the mean pressure of the cold neutral medium obtained by Jenkins & Tripp (2011).</figcaption></figure>

## 5.2. Statistical Properties of the Turbulent Energy Cascade

The separation of the thermal and turbulent contribution to the line width, as well as the estimate of the depth of the line of sight, allow us to quantify the statistical properties of the turbulent energy cascade in the WNM of NEP.

### 5.2.1. Inertial Range

The spatial power spectra of $N_H(r)$ and $C(r)$ shown in Figure 11 (see Appendix A for the procedure) are both well described by single power laws. This confirms the existence of an inertial range in NEP, over the projected spatial range of the observation from a few parsecs up to $L_x = L_y$. Note that the depth along the line of sight, $L_z = 130$ pc, where the turbulent line width $\sigma_{vz}(L_z)$ is measured, is significantly larger than the largest-scale measured in projection on the sky $L_x = L_y = 63$ pc. From Figure 11 alone it is impossible to ensure that an inertial range exists up to $l \sim L_z$. The back of the WNM layer studied here is at a distance of ~400 pc (see Figure 6) which, at the latitude of NEP, corresponds to a height of about 200 pc. The WNM volume under study here is thus within the WNM layer (Dickey & Lockman 1990; HWHM ~265 pc). As it is likely that the outer scale of interstellar turbulence is of the order of, or larger than, the disk thickness (Wolfire et al. 2003), in the following we will assume that the inertial range of turbulence goes up to $l \sim L_z$.

The slopes of the spatial power spectra of $N_H(r)$ and $C(r)$ are -3.1 ± 0.1 and -3.3 ± 0.1, respectively. Martin et al. (2015) performed a similar analysis using a multiphase decomposition of the emission of the local gas in NEP defined between -20.5 < $v_z$(km s⁻¹) < 47.9. They found an exponent -2.7 ± 0.1 for the column density, slightly shallower than that obtained here. Note however, that the WNM modeled by Martin et al. (2015) is a combination of our WNM and LNM phases, which could explain this difference.

The two slopes obtained here for $N_H(r)$ and $C(r)$ are flatter than the expected -11/3 value for compressible sub/trans-sonic turbulence (Kim & Ryu 2005). As noted in Section 4.2, the effect of a partial volume-filling factor of the WNM could

<footnote n="7">⁷ These variations can be due to the distribution of stars in V which locally increase or decrease $G_0$, and/or variations due to the density field topology (Parravano et al. 2003). The latter induce a natural shielding of the radiation field, causing a decrease proportional to the column density of the gas.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Spatial power spectrum model of the column density map (left) and centroid velocity field (right) of the WNM in NEP. $P(k)$ are given by the black dots. Gray points represent the noise component $N(k)$. Red points show the corrected power spectrum $(P(k) - N(k))/B(k)$. The dotted line shows the model fits between the vertical dashed lines, chosen to select $k$ ranges avoiding noise or systematics (high- and low-$k$ limits, respectively). The orange line shows the result of the global modeling of Equation (A1). Exponents $\gamma_N = -3.04 \pm 0.08$ and $\gamma_{(v_z)_z} = -3.25 \pm 0.06$.</figcaption></figure>

explain that difference. To investigate this, we used fractional Brownian motion simulations to which we applied partial volume-filling factors (see Appendix B). We found that whatever the statistical properties of the masking field that simulate $f_v^{WM}$ < 1, the power spectrum slope of its projection is flatter than the true 3D value. We note, though, that the effect is less dramatic for the velocity centroid than for the column density field.

Therefore, from the power spectrum of the centroid velocity (or column density) itself, it is impossible to conclude on the nature of the turbulence, especially on its turbulent Mach number. Nevertheless we note that our results are compatible with $\gamma_v = \gamma_N = -11/3$, the expected value for compressible sub/trans-sonic turbulence.

### 5.2.2. Turbulent Sonic Mach Number and Reynolds Number

We quantify the strength of turbulence using the turbulent sonic Mach number and the turbulent Reynolds number. Note that these two are functions of scale. The turbulent sonic Mach number is

$$ \mathcal{M}_s(L_z) = \frac{\sqrt{3} \sigma_{v_z}(L_z)}{C_s} \qquad (19) $$

where the adiabatic sound speed is

$$ C_s = \sqrt{\frac{\gamma k_B T_k}{\mu m_H}}, \qquad (20) $$

where $\mu = 1.4$ is the molecular weight for the atomic Galactic composition, and $\gamma = 5/3$ the adiabatic index of monoatomic gas. Using Equation (19), we find $\mathcal{M}_s(L_z) = 0.87 \pm 0.15$, showing that the WNM in NEP is in a subsonic regime at scale $L_z$. This result is reminiscent of the low Mach number found for the WIM by Gaensler et al. (2011).

An important property of the turbulent velocity field is that it follows a scaling law dictated by the nature of the turbulent

cascade and controlled by an exponent $q$. In a sub/trans-sonic regime, we expect $q = 1/3$ (equivalent to $-11/3$ for the power spectrum scaling law). The so-called $\sigma_v - l$ relation is

$$ \sigma_{v_z}(l) = \sigma_{v_z}(1) l_{pc}^{1/3}, \qquad (21) $$

where $\sigma_{v_z}(1) = \sigma_{v_z}(L_z)L_z^{-1/3}$ is the velocity dispersion of the turbulent velocity field along the line of sight at 1 pc, and $l_{pc}$ is the scale in parsecs. Here we find that $\sigma_{v_z}(1) = 0.76 \pm 0.1$ km s$^{-1}$. Interestingly, this value is equivalent to the average found in much denser regions of the ISM; Miville-Deschênes et al. (2017a) found $\sigma_{v_z}(1) \sim 0.8$ km s$^{-1}$ for molecular clouds in the solar neighborhood. This echoes the insight of Hennebelle & Falgarone (2012) who mentioned that molecular clouds and the H I are likely to be part of the same turbulent cascade. Similarly, combining Equations (19) and (21), it follows that

$$ \mathcal{M}_s(l) = \mathcal{M}_s(1) l_{pc}^{1/3}, \qquad (22) $$

which gives $\mathcal{M}_s(1) = 0.17 \pm 0.03$.

In the following, we assume a mean free path $\lambda = 1/\sigma n$ with $\sigma = 1 \times 10^{-15}$ cm$^2$ for the hydrogen cross-section (Lequeux 2012) and a thermal velocity $v_{th} = \sqrt{8/\pi} C_s$ to find the kinematic viscosity

$$ \nu = \frac{1}{3} \lambda v_{th}. \qquad (23) $$

Values of these properties are given in Table 3.

For an isothermal gas, the Knudsen number $Kn(L_z) = \lambda/L_z$ and $\mathcal{M}_s$ are linked directly to the Reynolds number

$$ Re(L_z) = \sqrt{\frac{\pi}{2} \frac{\mathcal{M}_s(L_z)}{Kn(L_z)}}. \qquad (24) $$

This quantifies the ratio of advection and diffusion in a turbulent fluid. We find $Re(L_z) = (3.3 \pm 1.8) \times 10^5$ (see also Table 3 for
the corresponding Knudsen numbers), corresponding to a regime of fully developed turbulence at scale $L_z$ (Reynolds 1883).

### 5.2.3. The Sonic Scale

Using Equation (22), the sonic scale is

$$ \lambda_s = L_z M_s^{-3} (L_z), \qquad (25) $$

which satisfies $M_s(\lambda_s) = 1$. We find $\lambda_s = (1.9 \pm 1.1) \times 10^2$ pc. Above this scale, the turbulent cascade starts to be supersonic and is expected to be described by $q = 0.5$ (equivalent $P(k) \propto k^{-4}$). A better sampling of high spatial scales could confirm a break in the scaling law of the projected velocity field and allow a parallel and complementary determination of $\lambda_s$.

### 5.2.4. Relative Strength of Solenoidal to Compressive Modes

Following Federrath et al. (2008), we estimate the relative strength of solenoidal to compressive modes $\zeta$ assuming

$$ \sigma_{\rho/\rho_0}^2 = b^2 M_s^2 (L_z) \qquad (26) $$

where $b = 1 - (2/3)\zeta$ for a 3D geometry. We find $\zeta = 0.5 \pm 0.4$. This value is compatible with a mixture of solenoidal and compressive modes, close to the natural state for which the energy is naturally distributed between the two ($\zeta = 0.5$).

A method to estimate $\zeta$ from hyper-spectral observations was proposed by Brunt & Federrath (2014) and recently applied in the Orion B molecular cloud by Orkisz et al. (2017). A similar application on NEP is beyond the scope of this paper but would provide a direct comparison to the value obtained here.

### 5.2.5. Energy Transfer and Dissipation

Here we consider two mechanisms of energy transfer and dissipation along the turbulent cascade: the kinematic molecular diffusion, and the ambipolar diffusion.

The dissipation scale $\eta$, on which the smallest eddies dissipate the turbulent energy into heat through viscosity, is

$$ \eta = L_z Re(L_z)^{-3/4}. \qquad (27) $$

We find $\eta = (1.0 \pm 0.4) \times 10^{-2}$ pc which is an order of magnitude higher than the mean free path of hydrogen atoms $\lambda$. The dissipation timescale $t_\eta$ is

$$ t_\eta = \frac{\eta^2}{\nu}. \qquad (28) $$

As Table 3 shows, we find $t_\eta = 50 \pm 18$ kyr. Finally, combining the dissipation scale and the dissipation time, the energy transfer rate $\epsilon$ is

$$ \epsilon = \frac{\eta^2}{t_\eta^3}. \qquad (29) $$

We find $\epsilon = (1.2 \pm 0.4) \times 10^{-4} L_\odot M_\odot^{-1}$. This is comparable to the value obtained by Hennebelle & Falgarone (2012) in the CNM ($\epsilon_{CNM} \sim 10^{-3} L_\odot M_\odot^{-1}$). As noted by these authors, $\epsilon_{CNM}$ is considerably lower than the energy transferred to the ISM via UV and visible radiation from stars. It is therefore also the case for the WNM.

We now consider the case where ambipolar diffusion due to ion-neutral friction is the dominant mechanism of energy dissipation. Following Miville-Deschênes et al. (2017b), the

ambipolar diffusion typical scale is

$$ l_{AD} = \sqrt{\frac{\pi}{\mu m_H}} \frac{B}{2X \langle \sigma v \rangle n^{3/2}} \qquad (30) $$

where $X = n_e/(n_e + n)$ is the ionization ratio with an electron density $n_e = 0.0213 \text{ cm}^{-3}$ (Berkhuijsen et al. 2006), $\langle \sigma v \rangle = 2 \times 10^{-9} \text{ cm}^3 \text{ s}^{-1}$ is the collision rate between ions and neutrals, assumed to be the Langevin rate, and $B = 6 \mu\text{G}$ is the typical value of the magnetic field strength in the WIM/WNM of the solar neighborhood (Beck 2001). We find $l_{AD} = (3.2 \pm 0.9) \times 10^{-2}$ pc. This is typically a factor 3–4 higher than the dissipation scale inferred from the molecular viscosity. Given the possible variations of ionization fraction or magnetic field strength, these two scales could be comparable. Therefore, it seems difficult to conclude which mechanism dominates the turbulent energy dissipation in the WNM of NEP.

## 6. Thermal Instability

It is believed that the condensation mode of TI is responsible for the thermal condensation of the warm neutral phase of the ISM that leads to the formation of CNM structures in galaxies. Building on Section 5, the thermal and turbulent properties of the WNM in NEP allow us to evaluate if perturbations around the mean thermodynamic state of the gas allow the condensation mode of TI to develop and then to grow freely. Values calculated in this section are tabulated in Table 3.

### 6.1. Development of the Condensation Mode

In an idealized non-viscous static fluid in thermal equilibrium, the isobaric criterion for development of the condensation mode of thermal instability can be expressed as

$$ \left(\frac{\partial P}{\partial n}\right)_{L=0} < 0 \qquad (31) $$

(Field 1965; Wolfire et al. 1995a). One might expect the WNM to be in a thermal pressure and volume density range which, if slightly compressed, can fall into the range where Equation (31) is satisfied. As Figure 10 shows, the thermodynamical state of the WNM in NEP is compatible with the development of the condensation mode of TI.

### 6.2. Can the Condensation Mode Grow Freely? Cooling and Dynamical Time Scales

When TI develops, condensation is possible if the cooling time of the fluid element is shorter than its dynamical time<fnref n="8" /> (Hennebelle & Pérault 1999):

$$ \mathcal{I}(l) = \frac{t_{\text{cool}}}{t_{\text{dyn}}(l)} < 1, \qquad (32) $$

where the cooling time is

$$ t_{\text{cool}} = \frac{k_B T_k}{(\gamma - 1) n \Lambda(T_k)}. \qquad (33) $$

<footnote n="8">⁸ Otherwise, the energy lost by radiation is small when compared to the increase of internal energy and the process can be considered as adiabatic (i.e., no transfer of heat between the fluid element and its surrounding medium).</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Dynamical time vs. scale (red). The top axis gives the corresponding turbulent sonic Mach number. The blue line and shaded area show the sonic scale and its error. Vertical and horizontal black dashed lines show the cooling length and cooling time, respectively.</figcaption></figure>

Here $\Lambda(T_k)$ is the net cooling function. The dynamical time of the compression is

$$t_{\text{dyn}}(l) = \frac{l_{\text{pc}}}{C_s}, \qquad (34)$$

where $l_{\text{pc}}$ is the linear size of the homogeneous compression.

In a turbulent fluid, compression events arise on a continuous range of scales characterized by the energy cascade. Figure 12 shows the dynamical time as function of scale (see the top axis for the corresponding turbulent sonic Mach number) which depends linearly on the sound speed of the WNM. The horizontal black dashed line shows the cooling time, $t_{\text{cool}} \sim 2.3$ Myr, based on the standard model ($N_c = 1 \times 10^{20}$ cm$^{-2}$ and $G_0 = 1.7$) from Wolfire et al. (2003). Above this line, $\mathcal{I}$ is satisfied, ensuring that TI can grow freely. The intersection between $t_{\text{dyn}}(l)$ and $t_{\text{cool}}$ corresponds to the linear scale that satisfies $\mathcal{I} = 1$ and is called the cooling length, $\lambda_{\text{cool}} \simeq 18$ pc.$^9$ In other words, compression arising on linear scales $l_{\text{pc}} > \lambda_{\text{cool}}$ also satisfies Equation (32). By construction, we find that $M_s(\lambda_{\text{cool}}) = 0.45$, showing that a too weak turbulent flow prevents TI growing freely to enable the formation of CNM structures, in accordance with numerical experiment of bistable turbulent flows (Saury et al. 2014).

For comparison, the blue line and shaded area show the sonic scale and its error $\lambda_s = (1.9 \pm 1.1) \times 10^2$ pc, which satisfies $M_s(\lambda_s) = 1$. At this scale, the information travels at sound speed $C_s$ and $t_{\text{dyn}}(\lambda_s) = 25 \pm 12$ Myr. We find $\mathcal{I}(\lambda_s) = 0.09 \pm 0.04$.

### 6.3. Field Length

During this condensation phase, along with heating and cooling, the fluid is subjected to thermal diffusion. The length at which radiative heating and cooling processes become comparable to thermal diffusion is the Field length$^{10}$

$$\lambda_{\text{F}} = \sqrt{\frac{\kappa(T_k) T_k}{n^2 \Lambda(T_k)}}, \qquad (35)$$

(Field 1965; Begelman & McKee 1990) where $\kappa = 2.5 \times 10^3 T_k^{1/2}$ erg cm$^{-1}$ K$^{-1}$ s$^{-1}$ is the thermal conductivity for hydrogen atoms (Parker 1953). We find $\lambda_F = 0.10 \pm 0.06$ pc. $\lambda_F$ is considerably lower than $\lambda_{\text{cool}}$, ensuring that compression of linear size, $l_{\text{pc}} > \lambda_{\text{cool}}$, will be more affected by cooling than thermal diffusion.

## 7. Summary

We have presented a study of the thermal and turbulent properties of the multiphase neutral ISM in the solar neighborhood based on the analysis of 21 cm emission data of a high Galactic latitude field (NEP) from the GHIGLS H I survey (Martin et al. 2015). In particular we have presented an original method to extract information from a 21 cm emission data cube, allowing us to constrain the physical properties of the WNM.

The specificity of this study relies on the following aspects.

1. The Gaussian decomposition tool ROHSA (Marchal et al. 2019) was used to provide a spatially coherent model of the local phases: WNM, LNM, and CNM. This allowed us to produce maps of the WNM column density and centroid velocity.

2. The distance and physical scales of the H I gas observed in projection were determined using a 3D dust extinction map.

3. The thermal and turbulent contribution to the WNM line widths were disentangled using the observed properties of its centroid velocity field and CNM structures within it.

Being able to isolate the WNM from the total H I emission cube, to locate where it is coming from along the line-of-sight, and to estimate the contribution of turbulent motions to the line width allowed us to estimate the physical properties of the WNM in great detail. The main conclusions are as follows.

1. The WNM average density ($n = 0.74 \pm 0.41$ cm$^{-3}$), kinetic temperature ($T_k = 6.0 \pm 1.3 \times 10^3$ K), and thermal pressure ($P_{\text{th}}/k_B = 4.4 \pm 2.6 \times 10^3$ K cm$^{-3}$) are compatible with expected values for the solar neighborhood conditions assuming typical cooling and heating processes.

2. The mass fraction of each phase shows significant variations over the $12^\circ \times 12^\circ$ field. On average the phase proportions in the local velocity component are CNM 8%, LNM 28%, and WNM 64%. We estimated the volume-filling fraction of the WNM to be $f_V^{\text{WNM}} = 0.6 \pm 0.2$.

3. The WNM has the properties of a trans-sonic turbulent flow, with low density contrast ($\sigma_{\rho/\rho_0} = 0.6 \pm 0.2$), density and velocity power spectra compatible with $k^{-11/3}$, and a low Mach number at the largest scale probed ($M_s(l = 130$ pc) = 0.87 $\pm$ 0.15). The WNM turbulent velocity dispersion at a scale of 1 pc ($\sigma_v(l) = 0.76 \pm 0.1$ km s$^{-1}$) is similar to what is observed in denser phases of the ISM.

4. Our determinations of the dynamic and static scales involved in the thermal condensation of the WNM confirm that the relatively low strength of turbulence coupled to the thermal state of the gas allow the

<footnote>$^9$ The cooling length is also the scale at which the WNM is nonlinearly unstable (Audit & Hennebelle 2005).</footnote>

<footnote>$^{10}$ The Field length is also called the conduction length and can be seen as the typical size of the front between the WNM and the CNM.</footnote>
condensation mode of TI to grow freely to form cold dense structures observed in the CNM.

5. We confirm the theoretical prediction of the cooling length $\lambda_{\text{cool}} \sim 18$ pc (Hennebelle & Pérault 1999) in the WNM of the solar neighborhood. In addition we place an observational constraint on the Field length, also coherent with theoretical expectations.

This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC). Part of this work was supported by Hyperstars, a project funded by the MASTODONS initiative of the CNRS mission for interdisciplinarity. This work took part under the program Milky Way-Gaia of the PSI2 project funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02. We gratefully acknowledge Peter G. Martin and Benjamin Godard for enlightening conversations, Claire Murray for providing us the result of the 21-Sponge survey analysis, and Ludovic Montier for providing us an Healpix version of the 3D dust map used in this work. We also gratefully acknowledge John Dickey and Robert Benjamin for reviewing this work as members of the A.M. thesis jury. We thank the anonymous referee whose comments and suggestions have improved this manuscript.

*Software:* `matplotlib` (Hunter 2007), `NumPy` (van der Walt et al. 2011), and `Astropy`,<fnref n="11" /> a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013, 2018).

## Appendix A
### Power Spectrum Analysis

Following Martin et al. (2015), the power spectrum $P(k)$ of $C(r)$ (eqv. $N_H(r)$) is the azimuthal average of the modulus of its Fourier transform, and is modeled as

$$P(k) = B(k) \times P_0 k^{\gamma} + N(k), \quad (A1)$$

where $P_0$ is the amplitude of the power spectrum, $\gamma$ is the scaling exponent, $B(k)$ is the beam of the instrument assuming a 2D Gaussian of FWHM = 9/4, and $N(k)$ is the noise estimated by taking the power spectrum of empty channels of the PPV cube. Finally, edges of the image are apodized using a cosine function to avoid effects due to the replication done by the Fourier transform algorithm.

## Appendix B
### Impact of Volume-filling Factor on the Statistics of Projected Quantities

The interpretation of power-law slopes of projected quantities, like the column density and the centroid velocity, has been the subject of several theoretical and numerical studies over the past 20 years (e.g., Miville-Deschênes et al. 2003). The general consensus is that the power spectrum of column density, estimated from optically thin lines, is a reliable proxy for the power spectrum of the 3D density field. As mentioned earlier, the centroid velocity is a reliable proxy for the velocity field only in the limit of small 3D density contrasts.

All the previous studies devoted to the understanding of the link between the statistical properties of projected and 3D fields are based on the assumption that the fluid that projects on the sky fills the 3D volume completely. In reality this is rarely

satisfied. For instance, the properties of molecular clouds are often analyzed using line tracers like $^{12}$CO or $^{13}$CO that depend strongly on the local gas volume density. This provides information only about the densest part of the fluid in 3D.

One could think that H I is less affected by this density effect and that the 21 cm emission would be a more reliable tracer of the whole 3D volume. This is most probably the case for the whole 21 cm emission but then the difficulty lies in the estimate of the statistical properties of a 3D multi-phase fluid with density contrasts of the order of 100–1000.

The separation of the thermal phases from the 21 cm emission data with ROHSA allows the extraction of the low-density WNM, opening the possibility to study interstellar turbulence by comparing its multi-scale statistics of density and velocity to those of controlled numerical experiments of isothermal turbulence. In this context, one interesting question is to what extent the fact that WNM does not occupy the full volume on the line of sight affects the relationship between the statistical properties of the 3D fields (density and velocity) and the projected quantities (column density and centroid velocity). To our knowledge, this has never been explored specifically.

The evaluation of the effect of a partial filling of a fluid on the projected quantities is clearly beyond the scope of this paper. In fact it is a rather ill-defined problem as one needs to define the statistics of the physical fields (density, velocity) as well as the 3D shape of the volume occupied by the fluid. That shape of that mask is likely to be defined by a density threshold in the case of CO emission or CNM structures. For the more diffuse inter-cloud H I medium (LNM and WNM), the shape of the volume is probably related to the temperature field, linked to heating and cooling processes. So it could be partly related to the density field but also to the radiation field intensity field. A proper study of this effect should be done with dedicated numerical simulations of the multi-phase ISM. In order to explore the main trends of the expected effects, we present a short study based only on fractional Brownian motion fields.

### B.1. Density Field

Let us first look at the effect on the column density. To do so we have constructed 3D density fields $n_H(r)$ over a grid of 128 × 128 × 256 pixels. We have chosen to produce 3D cubes with a depth ($L_z$) twice the size in [x, y] to mimic an observation like NEP. The statistics of $n_H(r)$ is assumed to follow a log-normal distribution, with $\delta(n_H)/\langle n_H \rangle = 1$, and a 3D power spectrum slope of $\gamma_n = -11/3$, reminiscent of what is seen for compressible, sub- or trans-sonic turbulent flows.

With such a low value of $\delta(n)/\langle n \rangle$, the power spectrum of the column density $N_H = \sum_z n_H dz$ integrated over the full volume has the same slope as the 3D density field: $\gamma_N = \gamma_n$. The question is: what is the power spectrum slope of $N_H$ if a fraction of the volume is removed from $n_H$? We have explored two possibilities: (1) a threshold in density and (2) a mask independent of density. In the first case we have gradually lowered to zero the density values above a given threshold. As the threshold is lowered, the volume-filling factor of the gas gets smaller.

In the second case, $n_H(r)$ is lowered to zero in the region based on a second independent field $m(r)$. At this point it is difficult to evaluate what should be the statistics of $m(r)$. We decided to use a positive, low-contrast ($\sigma(m)/\langle m \rangle = 1/5$) Gaussian field with a power spectrum slope of $-11/3$. As for the first case, the regions of $n_H$ are lowered to zero but here

<footnote n="11">¹¹ http://www.astropy.org</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Power spectrum slope of the column density resulting from the projection of 3D density fields with partial 3D sampling.</figcaption></figure>

we use the criteria on $m$, where $m > m_0$. As $m_0$ is lowered, larger and larger regions of $n_H$ decrease to zero, lowering the volume filling factor.

Figure 13 shows the power spectrum slope of $\gamma_N$ recovered in both cases, as a function of the volume-filling factor $f_v$. The black and red points correspond to case (1) and (2) respectively. As expected, when the fluid occupies the full volume ($f_v = 1$), we recover the result $\gamma_N = \gamma_n$. The main result we obtained is that, as the volume-filling factor is lowered, the power spectrum slope increases (the $P(k)$ is flatter). Interestingly it seems that $\gamma_N$ reaches an almost constant value for $f_v < 0.8$, but that value depends on the way the 3D mask is built.

The main result here is the fact that a volume-filling factor lower than unity introduces a systematic bias in the determination of the 3D density slope from the column density. A value of $f_v < 1$ systematically flattens the power spectrum. The power spectrum slope $\gamma_N$ obtained from observations provides an upper limit of the true value of $\gamma_n$. The effect is important. According to our little experiment it could be as much as $\gamma_n = \gamma_N - 1$. This has important consequences on the use of the column density power spectrum slope to constrain the properties of turbulence, the Mach number in particular (Kim & Ryu 2005).

## B.2. Velocity Field

We developed a similar experiment for the velocity field. We built uncorrelated 3D density and velocity cubes, $n(\mathbf{r})$ and $v_z(\mathbf{r})$, using the same method as previously. The velocity field has a Kolmogorov slope ($\gamma_v = -11/3$). The masking is done on $n$ by removing pixels in 3D where the density is above a given threshold $n > n(0)$. Then the centroid velocity field $C(\mathbf{r})$ is estimated by calculating the average velocity on each line of sight, using only the unmasked pixels. Figure 14 shows the power spectrum slope of $\gamma_C$ recovered in each case, as a function of the volume-filling factor $f_v$. As for the column density field case, the partial filling of the 3D volume introduces a bias in the projected velocity field. We note, though, that the effect is less dramatic for the velocity centroid than for the column density field.

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Power spectrum slope of the centroid velocity field as a function of the volume-filling factor of the fluid. The masked regions are defined using a gradual threshold on a log-normal density field with a power spectrum slope $\gamma_n$. The different curves show the result for a range of $\gamma_n$.</figcaption></figure>

## ORCID iDs

Antoine Marchal https://orcid.org/0000-0002-5501-232X

Marc-Antoine Miville-Deschênes https://orcid.org/0000-0002-7351-6062

## References

* Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

* Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

* Audit, E., & Hennebelle, P. 2005, A&A, 433, 1

* Beck, R. 2001, SSRv, 99, 243

* Begelman, M. C., & McKee, C. F. 1990, ApJ, 358, 375

* Begum, A., Stanimirovic, S., Goss, W. M., et al. 2010, ApJ, 725, 1779

* Bellomi, E., Godard, B., Hennebelle, P., et al. 2020, A&A, 643, A36

* Berkhuijsen, E. M., Mitra, D., & Mueller, P. 2006, AN, 327, 82

* Bohlin, R. C., Savage, B. D., & Drake, J. F. 1978, ApJ, 224, 132

* Brunt, C. M., & Federrath, C. 2014, MNRAS, 442, 1451

* Brunt, C. M., Federrath, C., & Price, D. J. 2010, MNRAS, 405, L56

* Carilli, C. L., Dwarakanath, K. S., & Goss, W. M. 1998, ApJL, 502, L79

* Cox, D. P. 2005, ARA&A, 43, 337

* Dickey, J. M., & Lockman, F. J. 1990, ARA&A, 28, 215

* Dickey, J. M., Salpeter, E. E., & Terzian, Y. 1979, ApJ, 228, 465

* Draine, B. T. 1978, ApJS, 36, 595

* Dwarakanath, K. S., Carilli, C. L., & Goss, W. M. 2002, ApJ, 567, 940

* Federrath, C., Klessen, R. S., & Schmidt, W. 2008, ApJL, 688, L79

* Field, G. B. 1965, ApJ, 142, 531

* Gaensler, B. M., Haverkorn, M., Burkhart, B., et al. 2011, Natur, 478, 214

* Gaensler, B. M., Madsen, G. J., Chatterjee, S., & Mao, S. A. 2008, PASA, 25, 184

* Green, G. M., Schlafly, E. F., Finkbeiner, D., et al. 2018, MNRAS, 478, 651

* Green, G. M., Schlafly, E. F., Finkbeiner, D. P., et al. 2014, ApJ, 783, 114

* Haud, U., & Kalberla, P. M. W. 2007, A&A, 466, 555

* Heiles, C., & Troland, T. H. 2003, ApJ, 586, 1067

* Hennebelle, P., & Audit, E. 2007, A&A, 465, 431

* Hennebelle, P., Banerjee, R., Vázquez-Semadeni, E., Klessen, R. S., & Audit, E. 2008, A&A, 486, L43

* Hennebelle, P., & Falgarone, E. 2012, A&ARv, 20, 55

* Hennebelle, P., & Péault, M. 1999, A&A, 351, 309

* Hunter, J. D. 2007, CSE, 9, 90

* Jenkins, E. B., & Tripp, T. M. 2011, ApJ, 734, 65

* Kalberla, P. M. W., & Haud, U. 2018, A&A, 619, A58

* Kanekar, N., Subrahmanyan, R., Chengalur, J. N., & Safouris, V. 2003, MNRAS, 346, L57

* Kim, J., & Ryu, D. 2005, ApJL, 630, L45
Lallement, R., Babusiaux, C., Vergely, J. L., et al. 2019, A&A, 625, A135

Lallement, R., Capitanio, L., Ruiz-Dern, L., et al. 2018, A&A, 616, A132

Lallement, R., Ferlet, R., Lagrange, A. M., Lemoine, M., & Vidal-Madjar, A. 1995, A&A, 304, 461

Lallement, R., Welsh, B. Y., Vergely, J. L., Crifo, F., & Sfeir, D. 2003, A&A, 411, 447

Lenz, D., Hensley, B. S., & Doré, O. 2017, ApJ, 846, 38

Lequeux, J. 2012, The Interstellar Medium (Berlin: Springer)

Levrier, F. 2004, A&A, 421, 387

Liszt, H. S. 2013, ApJ, 780, 10

Marchal, A., Miville-Deschênes, M.-A., Orieux, F., et al. 2019, A&A, 626, A101

Martin, P. G., Blagrave, K. P. M., Lockman, F. J., et al. 2015, ApJ, 809, 153

McKee, C. F., & Ostriker, J. P. 1977, ApJ, 218, 148

Mebold, U. 1972, A&A, 19, 13

Mebold, U., Winnberg, A., Kalberla, P. M. W., & Goss, W. M. 1982, A&A, 115, 223

Miville-Deschênes, M.-A., Levrier, F., & Falgarone, E. 2003, ApJ, 593, 831

Miville-Deschênes, M.-A., & Martin, P. G. 2007, A&A, 469, 189

Miville-Deschênes, M.-A., Murray, N., & Lee, E. J. 2017a, ApJ, 834, 57

Miville-Deschênes, M.-A., Salomé, Q., Martin, P. G., et al. 2017b, A&A, 599, A109

Murray, C. E., Lindner, R. R., Stanimirovic, S., et al. 2014, ApJL, 781, L41

Murray, C. E., Peek, J. E. G., & Kim, C.-G. 2020, ApJ, 899, 15

Murray, C. E., Stanimirović, S., Goss, W. M., et al. 2018, ApJS, 238, 14

Nguyen, H., Dawson, J. R., Lee, M.-Y., et al. 2019, ApJ, 880, 141

Nguyen, H., Dawson, J. R., Miville-Deschênes, M. A., et al. 2018, ApJ, 862, 49

Orkisz, J. H., Pety, J., Gerin, M., et al. 2017, A&A, 599, A99

Ossenkopf, V., Esquivel, A., Lazarian, A., & Stutzki, J. 2006, A&A, 452, 223

Parker, E. N. 1953, ApJ, 117, 431

Parravano, A., Hollenbach, D. J., & McKee, C. F. 2003, ApJ, 584, 797

Planck Collaboration, Abergel, A., Ade, P. A. R., et al. 2011, A&A, 536, A24

Planck Collaboration, Abergel, A., Ade, P. A. R., et al. 2014, A&A, 571, A11

Reynolds, R. J. 1977, ApJ, 216, 433

Reynolds, O. 1883, RSPTA, 174, 935

Rezaei, Kh., S., Bailer-Jones, C. A. L., Hanson, R. J., & Fouesneau, M. 2017, A&A, 598, A125

Saury, E., Miville-Deschênes, M.-A., Hennebelle, P., Audit, E., & Schmidt, W. 2014, A&A, 567, A16

Seifried, D., Schmidt, W., & Niemeyer, J. C. 2011, A&A, 526, A14

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, CSE, 13, 22

Vergely, J.-L., Ferrero, R. F., Siebert, A., & Valette, B. 2001, A&A, 366, 1016

Vergely, J.-L., Valette, B., Lallement, R., & Raimond, S. 2010, A&A, 518, A31

Wolfire, M. G., Hollenbach, D., McKee, C. F., Tielens, A. G. G. M., & Bakes, E. L. O. 1995a, ApJ, 443, 152

Wolfire, M. G., McKee, C. F., Hollenbach, D., & Tielens, A. G. G. M. 1995b, ApJ, 453, 673

Wolfire, M. G., McKee, C. F., Hollenbach, D., & Tielens, A. G. G. M. 2003, ApJ, 587, 278

Zhu, H., Tian, W., Li, A., & Zhang, M. 2017, MNRAS, 471, 3494