ROHSA: Regularized Optimization for Hyper-Spectral Analysis - Application to phase separation of 21 cm data

TL;DR

ROHSA is a novel Gaussian decomposition algorithm that leverages regularized optimization to analyze hyper-spectral 21 cm data, effectively revealing the multiphase structure of the interstellar medium and overcoming velocity blending issues.

Contribution

This paper introduces ROHSA, a new multi-resolution, regularized optimization algorithm for hyper-spectral analysis, specifically designed to extract multiphase HI structures from 21 cm data.

Findings

ROHSA accurately recovers the multiphase HI structure in synthetic data.

ROHSA reveals detailed multiphase features in real 21 cm observations.

The method effectively maps intermediate-velocity clouds in high Galactic latitude regions.

Abstract

Extracting the multiphase structure of the neutral interstellar medium (ISM) is key to understand the star formation in galaxies. The radiative condensation of the diffuse warm neutral medium producing a thermally unstable lukewarm medium and a dense cold medium is closely related to the initial step which leads the HI-to-H2 transition and the formation of molecular clouds. Up to now the mapping of these phases out of 21 cm emission hyper-spectral cubes has remained elusive mostly due to the velocity blending of individual cold structures present on a given line of sight. To address this problem, we developed a new Gaussian decomposition algorithm named ROHSA (Regularized Optimization for Hyper-Spectral Analysis) based on a multi-resolution process from coarse to fine grid. ROHSA uses a regularized non-linear least-square criterion to take into account simultaneously the spatial coherence of the emission and the multiphase nature of the gas. In order to obtain a solution with spatially smooth parameters, the optimization is performed on the whole data cube at once. The performances of ROHSA were tested on a synthetic observation computed from numerical simulations of thermally bi-stable turbulence. An application on a 21 cm observation of a high Galactic latitude region from the GHIGLS survey is presented. The evaluation of ROHSA on synthetic 21 cm observations shows that it is able to recover the multiphase nature of the HI and physically meaningful information about the underlying three-dimensional fields (density, velocity and temperature). The application on a real 21 cm observation of a high Galactic latitude field produces a picture of the multiphase HI, with isolated, filamentary and narrow structures and wider, diffuse and space filling components. The test-case field used here contains significant intermediate-velocity clouds that were well mapped out by the algorithm.

Figures

Tables

Table 1. Table 1: Mean velocity ⟨ 𝝁 n ⟩ delimited-⟨⟩ subscript 𝝁 𝑛 \langle\bm{\mu}_{n}\rangle and mean velocity dispersion ⟨ 𝝈 n ⟩ delimited-⟨⟩ subscript 𝝈 𝑛 \langle\bm{\sigma}_{n}\rangle of the eight Gaussian components G n subscript 𝐺 𝑛 G_{n} inferred by ROHSA on the 21 cm synthetic observation of the numerical simulation described in Sect.  3.2

	$G_1$	$G_2$	$G_3$	$G_4$	$G_5$	$G_6$	$G_7$	$G_8$
$⟨{𝝁}_n⟩$ [km s-1]	0.2	-1.7	0.5	2.3	-0.1	5.0	-3.7	-2.5
$⟨{𝝈}_n⟩$ [km s-1]	8.2	6.1	0.5	0.6	1.5	1.5	1.8	0.9

Table 2. Table 2: Mean velocity ⟨ 𝝁 n ⟩ delimited-⟨⟩ subscript 𝝁 𝑛 \langle\bm{\mu}_{n}\rangle and mean velocity dispersion ⟨ 𝝈 n ⟩ delimited-⟨⟩ subscript 𝝈 𝑛 \langle\bm{\sigma}_{n}\rangle of the twelve Gaussian components G n subscript 𝐺 𝑛 G_{n} inferred by ROHSA in NEP.

	$G_1$	$G_2$	$G_3$	$G_4$	$G_5$	$G_6$	$G_7$	$G_8$	$G_9$	$G_{10}$	$G_{11}$	$G_{12}$
$⟨{𝝁}_n⟩$ [km s-1]	-74.1	-53.9	-44.7	-35.0	-22.9	-12.6	-4.8	-1.3	0.2	10.9	40.8	75.9
$⟨{𝝈}_n⟩$ [km s-1]	9.7	6.2	3.5	5.1	5.3	4.7	1.6	6.3	1.9	7.5	12.5	9.3

Equations

$$\tilde{T}_{B}\big{(}v_{z},\bm{\theta}(\bm{r})\big{)}=\sum_{n=1}^{N}G\big{(}v_{z},\bm{\theta}_{n}(\bm{r})\big{)}$$

$$G\big{(}v_{z},\bm{\theta}_{n}(\bm{r})\big{)}=\bm{a}_{n}(\bm{r})\exp\left(-\frac{\big{(}v_{z}-\bm{\mu}_{n}(\bm{r})\big{)}^{2}}{2\bm{\sigma}_{n}(\bm{r})^{2}}\right)$$

$$L\big{(}v_{z},\bm{\theta}(\bm{r})\big{)}=\tilde{T}_{B}\big{(}v_{z},\bm{\theta}(\bm{r})\big{)}-T_{B}(v_{z},\bm{r})$$

$$Q(\bm{\theta})=\frac{1}{2}\,\norm{L\big{(}v_{z},\bm{\theta}\big{)}}_{\bm{\Sigma}}^{2}=\frac{1}{2}\,\sum_{v_{z},\bm{r}}\left(\frac{L\big{(}v_{z},\bm{\theta}(\bm{r})\big{)}}{\bm{\Sigma}(\bm{r})}\right)^{2}$$

$$d=\left[\begin{array}{ccc}0& -1& 0\\ -1& 4& -1\\ 0& -1& 0\end{array}\right].$$

$$R(\mathbf{\theta })=\frac{1}{2}\sum _{n=1}^N{\lambda }_{\mathbf{a}}\Vert \mathbf{D}{\mathbf{a}}_n{\Vert }_2^2+{\lambda }_{\mathbf{\mu }}\Vert \mathbf{D}{\mathbf{\mu }}_n{\Vert }_2^2+{\lambda }_{\mathbf{\sigma }}\Vert \mathbf{D}{\mathbf{\sigma }}_n{\Vert }_2^2$$

$$R(\mathbf{\theta },\mathbf{m})=\frac{1}{2}\sum _{n=1}^N$$

$$+{\lambda }_{\mathbf{\sigma }}^{\prime }\Vert {\mathbf{\sigma }}_n-m_n{\Vert }_2^2$$

$$J(\mathbf{\theta },\mathbf{m})=Q(\mathbf{\theta })+R(\mathbf{\theta },\mathbf{m})$$

$$[\hat{\mathbf{\theta }},\hat{\mathbf{m}}]=\underset{\mathbf{\theta },\mathbf{m}}{\mathrm{argmin}}J(\mathbf{\theta },\mathbf{m}),\text{wrt. }{\mathbf{a}}_n\ge 0,\forall n\in [1,\dots ,N]$$

$$\nabla J(\mathbf{\theta },\mathbf{m})=\left[\begin{array}{c}\nabla L(\mathbf{\theta })\times L(\mathbf{\theta })\\ 0\end{array}\right]+\left[\begin{array}{c}{\nabla }_{\mathbf{\theta }}R(\mathbf{\theta },\mathbf{m})\\ {\nabla }_{\mathbf{m}}R(\mathbf{\theta },\mathbf{m})\end{array}\right]$$

$${\mathbf{\theta }}_{(k+1)}={\mathbf{\theta }}_{(k)}-{\alpha }_{(k)}{\mathbf{H}}_{(k)}^{-1}\nabla J\left({\mathbf{\theta }}_{(k)},{\mathbf{m}}_{(k)}\right)$$

$$⟨T_B{⟩}_i=\frac{1}{K_i}\sum _{k\in {\mathcal{V}}_i}{T}_B(v_z,{\mathbf{r}}_k)$$

$${\varphi }_{{\nu }_z}(\text{x})d{v}_z^{\prime }=\frac{1}{\sqrt{2\pi }\Delta (\text{x})}\times \exp \left(-\frac{(v_z^{\prime }-v_z(\text{x}){)}^2}{2{\Delta }^2(\text{x})}\right)d{v}_z^{\prime }$$

$$T_b(v_z,\mathbf{r})=\sum _{z}T(\mathbf{r},z)\left[1-e^{-\tau (v_z,\mathbf{r},z)}\right]e^{-{\sum }_{z^{\prime }<z}\tau (v_z,\mathbf{r},z^{\prime })}$$

$$\tau (v_z,\mathbf{r},z)=\frac{1}{C}\frac{\rho (\mathbf{r},z){\varphi }_{{\nu }_z}(\mathbf{r},z)}{T(\mathbf{r},z)}dz$$

$$T_B^{\mathrm{thin}}(v_z,\mathbf{r})d{v}_z^{\prime }=B(\mathbf{r})$$

$$N_{\mathrm{HI}}^{\mathrm{thin}}(\mathbf{r})=C{\int }_{-\infty }^{+\infty }{T}_B(v_z,\mathbf{r})d{v}_z$$

$$⟨v_z(\mathbf{r})⟩=\frac{{\int }_{-\infty }^{+\infty }{v}_z{T}_B(v_z,\mathbf{r})d{v}_z}{{\int }_{-\infty }^{+\infty }{T}_B(v_z,\mathbf{r})d{v}_z}.$$

$$\displaystyle\nabla L(v_{z},\bm{\theta}(\bm{r}))=\begin{bmatrix}\exp\left(-\frac{\big{(}v_{z}-\mu_{1}(\bm{r})\big{)}^{2}}{2\sigma_{1}(\bm{r})^{2}}\right)\\ \frac{a_{1}(v_{z}-\mu_{1}(\bm{r}))}{\sigma_{1}^{2}}\exp\left(-\frac{\big{(}v_{z}-\mu_{1}(\bm{r})\big{)}^{2}}{2\sigma_{1}(\bm{r})^{2}}\right)\\ \frac{a_{1}(v_{z}-\mu_{1}(\bm{r}))^{2}}{\sigma_{1}^{3}}\exp\left(-\frac{\big{(}v_{z}-\mu_{1}(\bm{r})\big{)}^{2}}{2\sigma_{1}(\bm{r})^{2}}\right)\\ \vdots\\ \frac{a_{N}(v_{z}-\mu_{N}(\bm{r}))^{2}}{\sigma_{N}^{3}}\exp\left(-\frac{\big{(}v_{z}-\mu_{N}(\bm{r})\big{)}^{2}}{2\sigma_{1}(\bm{r})^{2}}\right)\end{bmatrix}$$

$${\nabla }_{\mathbf{\theta }}R(\mathbf{\theta },\mathbf{m})=\left[\begin{array}{c}{\lambda }_{\mathbf{a}}{\mathbf{D}}^t\mathbf{D}{\mathbf{a}}_1\\ {\lambda }_{\mu }{\mathbf{D}}^t\mathbf{D}{\mathbf{\mu }}_1\\ {\lambda }_{\mathbf{\sigma }}{\mathbf{D}}^t\mathbf{D}{\mathbf{\sigma }}_1\\ {\lambda }_{\mathbf{\sigma }}^{\prime }({\mathbf{\sigma }}_1-m_1)\\ ⋮\\ {\lambda }_{\mathbf{\sigma }}^{\prime }({\mathbf{\sigma }}_N-m_N)\end{array}\right]$$

$${\nabla }_{\mathbf{m}}R(\mathbf{\theta },\mathbf{m})=\left[\begin{array}{c}-{\sum }_{\mathbf{r}}{\lambda }_{\mathbf{\sigma }}^{\prime }({\mathbf{\sigma }}_1-m_1)\\ ⋮\\ -{\sum }_{\mathbf{r}}{\lambda }_{\mathbf{\sigma }}^{\prime }({\mathbf{\sigma }}_N-m_N)\end{array}\right]$$

Full text

*Aims. *The goal of this work was to develop a new algorithm to perform separation of diffuse sources in hyper-spectral data. Specifically the algorithm was designed in order to address the velocity blending problem by taking advantage of the spatial coherence of the individual sources. The main scientific driver of this effort was to extract the multiphase structure of the HI from 21 cm line emission only, providing a mean to map each phase separately, but the algorithm developed here should be generic enough to extract diffuse structures in any hyper-spectral cube.

*Methods. *We developed a new Gaussian decomposition algorithm named ROHSA (Regularized Optimization for Hyper-Spectral Analysis) based on a multi-resolution process from coarse to fine grid. ROHSA uses a regularized non-linear least-square criterion to take into account simultaneously the spatial coherence of the emission and the multiphase nature of the gas. In order to obtain a solution with spatially smooth parameters, the optimization is performed on the whole data cube at once. The performances of ROHSA were tested on a synthetic observation computed from numerical simulations of thermally bi-stable turbulence. An application on a 21 cm observation of a high Galactic latitude region from the GHIGLS survey is presented.

Results. The evaluation of ROHSA on synthetic 21 cm observations shows that it is able to recover the multiphase nature of the HI. For each phase, the power spectra of the column density and centroid velocity are well recovered. More generally that test reveals that a Gaussian decomposition of HI emission is able to recover physically meaningful information about the underlying three-dimensional fields (density, velocity and temperature). The application on a real 21 cm observation of a high Galactic latitude field produces a picture of the multiphase HI, with isolated, filamentary and narrow ($\sigma\sim 1-2$ km s-1) structures and wider ($\sigma\sim 4-10$ km s*-1*), diffuse and space filling components. The test-case field used here contains significant intermediate-velocity clouds that were well mapped out by the algorithm. As ROHSA is designed to extract spatially coherent components, it performs well at projecting out the noise.

*Conclusions. *In this paper we are introducing ROHSA a new algorithm performing a separation of diffuse sources in hyper-spectral data on the basis of a Gaussian decomposition. The algorithm makes no assumption about the nature of the sources, except that each one as a similar line-width. The tests we made shows that ROHSA is well suited to decompose complex 21 cm line emission of high Galactic latitude regions, but its design is general enough that it could be applied to any hyper-spectral data type for which a Gaussian model is relevant.

1 Introduction

Star formation in galaxies is strongly linked to the physical processes that govern the evolution of the interstellar medium. Stars form by gravitational collapse of dense ($n>10^{4}$ cm*-3*) and cold ($T\sim 10$ K) structures in molecular clouds but the process that leads to the formation of these over-densities is still unclear. One key element seems to be related to the efficiency of the formation of cold clouds of neutral hydrogen (HI) (Ostriker et al. 2010).

The current vision of the HI comes from an important legacy. Early observations of the 21 cm line showed a significant difference between emission and absorption spectra. On lines of sight crossing radio-sources the HI appears in absorption with very narrow features (a few km s*-1*). In emission the 21 cm line contains these narrow features on top of much boarder spectral structures (10-20 km s*-1*). Clark (1965) was the first to suggest that this might be the signature of a cloud-intercloud medium in pressure equilibrium. Very rapidly, Field (1965); Field et al. (1969) introduced the concept of thermal instability and laid out the theoretical ground of a ”two phase” HI model showing that, at the pressure of the ISM, the heating and cooling processes naturally lead to two thermally stable state : a dense cold neutral medium (CNM - $T\sim 50$ K, $n\sim 50$ cm*-3*) immersed in a diffuse warm neutral medium (WNM - $T\sim 8000$ K, $n\sim 0.3$ cm*-3*). This vision was later complemented by Wolfire et al. (1995, 2003) considering updated heating (dominated by the photo-electric effect on small dust grains) and cooling (dominated by CII - 158 $\mu$m, OI - 63 $\mu$m, L$\alpha$ and electron recombinations onto positive charged grains) processes of the ISM.

This description of the diffuse neutral gas completed a parallel point of view that emerged in the 50’s (e.g., von Weizsäcker 1951; von Hoerner 1951; Chandrasekhar & Münch 1952) and that considered the ISM as a multi-scale turbulent medium. In this case the density and velocity structures are the result of a highly dynamical and out-of-equilibrium medium. In order to reconcile the static/two-phase and the turbulent point of views, several studies have aimed at understanding the production of the CNM in a turbulent and thermally unstable flow using numerical simulations (e.g., Hennebelle & Pérault 1999; Koyama & Inutsuka 2002; Audit & Hennebelle 2005; Hennebelle et al. 2008; Saury et al. 2014). In general, these numerical studies show that the WNM has the properties of a trans-sonic turbulent flow, while the CNM shows a much more contrasted density structure, in accordance with the cloud-intercloud picture. In addition, such studies indicate the presence of a significant fraction of the mass being in the thermally unstable regime, (i.e., with a temperature mid-way between the CNM and WNM stable states). For instance, Saury et al. (2014) showed that 30% of the HI is in the thermally unstable regime. Interestingly they also show that this lukewarm neutral medium (LNM) is spatially located around the cold structures, pointing at the transitional nature of this thermal state.

From the observational standpoint, studies combining 21 cm absorption and emission data have clearly revealed the presence of HI at intermediate/unstable temperatures, typically between 500 and 5000 K (e.g., Heiles & Troland 2003a; Kanekar et al. 2003; Roy et al. 2013a, b; Murray et al. 2015, 2018b). Based on a coherent modeling of emission and absorption spectra Heiles & Troland (2003b); Murray et al. (2015, 2018b) estimated that about 30% of the HI is in the cold CNM phase, 20% in the thermally unstable regime and 50% in the WNM. Nevertheless the fraction of the HI in each phase remains uncertain and large variations are observed: the fraction of the mass in the CNM ranges from $\sim$1% to more than 50% Murray et al. (2018b).

The nature of these variations and how they relate to the dynamical conditions of the gas remains largely unexplored from the observational point of view. One main hurdle in getting access to this information is the fact that our knowledge of the multiphase nature of the HI relies on 21 cm absorption measurements that are limited to lines of sight crossing radio-sources. By nature, this way of observing prevents us from mapping the HI phases. To go further, and really compare with numerical simulation that are, for now, under-constrained by observation, it is mandatory to map the column density structure of each phase and study the spatial variations of their centroid velocity and velocity dispersion. This calls for methods that can extract the information of each HI phase from fully sampled 21 cm emission data only. Huge efforts have been made to map the 21 cm emission of the Galactic HI (recent examples are Taylor et al. 2003; Kalberla et al. 2005; Stil et al. 2006; McClure-Griffiths et al. 2009; Winkel et al. 2016; Peek et al. 2018) and a large amount of data is now available. The information about the multiphase and multi-scale nature of the HI contained in these large hyper-spectral data cubes has remained elusive due to the difficulty to separate the emission from the different phases on each line of sight. In this paper we propose a new method to map out the contribution of each phase to the 21 cm emission. The method is based on a decomposition of the 21 cm emission line with Gaussian profiles and constraints that favors spatially coherent parameters.

This paper is organized as follows. In Sect. 2, we describe the methodology used to develop our Gaussian decomposition algorithm. Evaluation on numerical simulation is presented Sect. 3. In Sect. 4 we present an application on observation. The discussion and summary are presented Sects. 5 and 6..

2 Methodology

2.1 Gaussian decomposition of the 21 cm emission

Very early on after its detection, the 21 cm line was observed to be well described by a sum of a small number of Gaussian components. This was found to be true for the least confused absorption spectra (Muller 1957, 1959; Clark 1965) but also for emission spectra observed away from the Galactic plane (Heeschen 1955). In fact very few spectra at high Galactic latitudes do not comply with that rule whatever the angular resolution of the data. Recently Kalberla & Haud (2018) showed that more than 60% of the spectra over the whole sky can be described by the sum of 7 or less Gaussian components. In the two decades after the detection of the 21 cm line, many studies used the Gaussian decomposition to infer physical parameters from the data (Matthews 1957; Davis 1957; Muller 1957, 1959; Dieter 1964, 1965; Lindblad 1966; Takakubo & van Woerden 1966; Mebold 1972). The fact that a small number of Gaussian components is needed to describe the signal was seen as a convenient way to describe the emission profiles with a small set of parameters (Takakubo & van Woerden 1966). It is also a very strong element in favor of the Gaussian function as a significant descriptor of the underlying physics.

Takakubo (1967) showed that the width of the 21 cm emission line could be grouped in three components ($\sigma\leq 3$ km s*-1*; $3<\sigma<7$ km s*-1*; $\sigma\geq 7$ km s*-1*), a result confirmed later on (Mebold 1972; Haud & Kalberla 2007; Kalberla & Haud 2018). Takakubo (1967) also showed that the narrow 21 cm features are well correlated with Ca+ K line absorption measurements. They concluded that the narrow component, also seen in 21 cm absorption spectra, are likely to be isolated cold clouds (CNM) in the Solar neighborhood. They also showed that the spatial distribution of the centroid velocity, velocity width and column density of the large feature is compatible with a warm (WNM) and diffuse disk that follows Galactic rotation. The second group of Gaussians, with a width ($\sigma$) between 3 and 7 km s*-1*, is generally attributed to gas in the thermally unstable range, but a fraction of them could be caused by blending of narrow features.

The exact mass fraction of gas in each phase (CNM, LNM and WNM) is still a matter of debate. In addition, because this knowledge is based on absorption measurements, there is very little information about the structure on the sky of these phases. Being able to separate the different phases on each line of sight would allow us to study the structure and kinematic of the cold phase and its relationship with the more diffuse gas. In theory, one could expect that Gaussian decomposition of the emission spectra could provide such mapping.

2.2 Limitation of the Gaussian model

There are many pitfalls in the description of the 21 cm emission spectra as a sum of Gaussian components : velocity blending, ambiguities of the number of components, non Gaussian profiles, noise peaks, non-uniqueness of the solution. Another important limitation is the effect of optical depth of the 21 cm line that modifies the shape of the line. More generally, the main opposition to Gaussian decomposition of emission spectra is that any spectrum can be decomposed that way provided that enough Gaussians are used. If that is the case, how can one be sure that the Gaussian representation provides some real physical information about the emitting gas ? For instance, two spatially disconnected cold structures present on the same line of sight could appear at the same projected velocity. In this case their respective emission profiles would be confused. This effect of velocity blending is affecting both the emission and absorption spectra. It is present in the data, more so at lower Galactic latitudes where the length of the line of sight is larger and the number of HI structures increases.

This line of reasoning led Dickey & Lockman (1990) to advise against using Gaussian decomposition to analyze 21 cm spectra. Later on, it continued to be used (Verschuur & Schmelz 1989; Verschuur & Magnani 1994; Poppel et al. 1994; Haud 2000; Verschuur 2004; Begum et al. 2010; Martin et al. 2015; Kalberla & Haud 2018) but overall there was a lost of interest, except for the analysis of absorption spectra. Indeed by comparing nearby absorption and emission spectra one can recover the effect of optical depth on the 21 cm emission. In addition, absorption measurements are only sensitive to cold gas, limiting the velocity blending problem. For these reasons, the Gaussian decomposition continued to be used in this context (Dickey et al. 2003; Kanekar et al. 2003), and especially after the seminal work of Heiles & Troland (2003b) who developed a dedicated formalism which was then used in several other studies since (e.g. Stanimirović & Heiles 2005; Begum et al. 2010; Stanimirović et al. 2014; Lee et al. 2015; Murray et al. 2014, 2015, 2018a). In fact key information about the nature of the HI come from the joint Gaussian decomposition of emission and absorption spectra.

2.3 Development of a new approach

Following what has been done for the comparison of emission and absorption spectra where the Gaussian decomposition is considered valid, we would like to argue that a similar decomposition could be envisaged for emission data only, at least at high Galactic latitudes where the effect of optical depth of the 21 cm line has been shown to be negligible (Murray et al. 2018a).

The fact that physical information could be obtained using a Gaussian decomposition of absorption data reveals the fact that thermal broadening has a significant effect in shaping the line profile or, in other words, that the dynamics of each HI phase is typical of sub or trans-sonic turbulence. When the amplitude of turbulent and thermal motions are commensurate, the line appears smooth and can be represented by a small number of Gaussian components (Miville-Deschênes et al. 2003). If the HI at high Galactic latitude is indeed represented by a two phase medium with small, cold and trans-sonic structures immersed in a relatively low Mach number and warm diffuse phase, the Gaussian representation could bear significant physical information.

The perspective of mapping the phases of the HI is so important that we ventured to explore new ways of decomposing the emission spectra that could be applicable to the high Galactic sky. The main difficulty is the effect of velocity blending for cold structures. As mentioned by Takakubo & van Woerden (1966), there will always be cases where a given spectrum can be fitted with a smaller number of components than its neighbors, if two or more components have similar central velocity and velocity dispersion. One way to avoid this confusion is to look for solutions that have a spatial continuity, or that have a slow spatial variation. Poppel et al. (1994), Haud (2000), Martin et al. (2015) and Miville-Deschênes et al. (2017b) have implemented Gaussian decomposition methods that uses some information about their neighbors in order to favor spatially coherent solutions. Nevertheless, formally these algorithms do not force solutions to be spatially coherent, they simply provide initial guesses to the fit of a single spectrum based on the most likely solutions found in some area around. The optimization is not bound to this initial guess and it can always converge to another solution that would break the spatial smoothness of the parameter space.

The novelty of the algorithm we present here is that it is the first one that imposes the spatial coherence in the determination of the parameters. In order to do that, all the spectra of the emission cube are fitted at the same time. To make sure that the recovered parameters are spatially smooth, specific regularization terms are added to the cost function with non negativity constraints on the amplitude. This algorithm, called ROHSA, is described next.

2.4 ROHSA

ROHSA performs a regression analysis using a regularized nonlinear least-square criterion. We formulate in this section the Gaussian model used as well as the energy terms added to the cost function to take into account simultaneously the spatial coherence of the emission and the multiphase nature of the gas. The quasi-Newton algorithm, L-BFGS-B, used to perform the optimization is then briefly described. Finally, we formulate the algorithm performed by ROHSA based on a multi-resolution process from coarse to fine grid.

2.4.1 Model

The data are the measured brightness temperature $T_{B}(v_{z},\bm{r})$ at a given projected velocity $v_{z}$ across sky coordinates $\bm{r}$. The proposed model \tilde{T}_{B}\big{(}v_{z},\bm{\theta}(\bm{r})\big{)} is a sum of $N$ Gaussian G\big{(}v_{z},\bm{\theta}_{n}(\bm{r})\big{)}

[TABLE]

with \bm{\theta}(\bm{r})=\big{(}\bm{\theta}_{1}(\bm{r}),\dots,\bm{\theta}_{n}(\bm{r})\big{)} and where

[TABLE]

is parametrized by \bm{\theta}_{n}=\big{(}\bm{a}_{n},\bm{\mu}_{n},\bm{\sigma}_{n}\big{)} with $\bm{a}_{n}\geq\bm{0}$ the amplitude, $\bm{\mu}_{n}$ the position and $\bm{\sigma}_{n}$ the standard deviation 2D maps of the $n$-th Gaussian profile across the plan of sky. The residual is

[TABLE]

The estimated parameters $\hat{\bm{\theta}}$ are defined as the minimizer of a cost function that includes the sum of the squares of the residual

[TABLE]

where $\Sigma$ is the standard deviation 2D map of the noise assumed without spatially correlation. In practice this term is estimated using a sequence of empty velocity channels of $T_{B}(v_{z},\bm{r})$.

For each of the $N$ Gaussians, we want to obtain a spatially coherent solution, meaning that for each parameter, the values have to be close for neighboring lines of sight. This can be done by penalizing the small scale spatial fluctuations of each parameter, measured by the energy at high spatial frequencies. The considered high-pass filter is the second-order differences, that is the Laplacian filtering, defined by the 2D convolution kernel

[TABLE]

The following regularization term, containing energy terms is added to the cost function given in Eq. (4)

[TABLE]

where $\bm{D}$ is a matrix performing the 2D convolution using the kernel $d$ and $\lambda_{\bm{a}}$, $\lambda_{\bm{\mu}}$, and $\lambda_{\bm{\sigma}}$ are hyper-parameters than tune the balance between the different terms.

These terms ensure a positive spatial correlation of the model parameters for neighboring pixels. However, each term is free to have large variation across the field at larger scale. Since $\bm{\sigma}_{n}$ contains information about the gas thermodynamics, we design an additional term in the cost function to group Gaussians with similar $\sigma_{n}$. This is implemented in order to favour solution that would produce components ascribable to each phase : WNM, LNM or CNM. In order to do that we add a term, $\lambda^{\prime}_{\bm{\sigma}}\|\bm{\sigma}_{n}-m_{n}\|_{2}^{2}$, that constrains $\bm{\sigma}_{n}$ to be close to an unknown scalar value $m_{n}$. The full regularization term is then

[TABLE]

with $\bm{m}=(m_{1},\dots,m_{N})$ and $\bm{a}_{n}\geq 0$, $\forall\,n\in[1,\dots,N]$. The last two terms in Eq. (2.4.1), representative of a joined constraint imposed on $\sigma_{n}$, allows us to interpret a posteriori and simultaneously the morphology and the thermodynamical state of each component extracted from the data. The full cost function is then

[TABLE]

2.4.2 Optimization algorithm

Unlike $Q(\bm{\theta})$, each energy term proposed in Eq (6) involves linear dependencies on the parameters $\bm{\theta}$. The cost function Eq. (8) is therefore a regularized non-linear least-square criterion. The minimizer

[TABLE]

as no closed form expression and is not directly tractable because of the complexity of the model $\tilde{T}_{B}$ and the size of the unknown and data. The proposed solution relies instead on an iterative optimization algorithm that uses the gradient

[TABLE]

which is tractable since it involves the residual, the Jacobian of the residual $\nabla L(\bm{\theta})$, and 2D convolutions with the kernel $d$ for $\bm{D}$ and $\bm{D}^{t}$. The gradient $\nabla R^{t}(\bm{\theta},\bm{m})=[\nabla_{\bm{\theta}}R^{t}(\bm{\theta},\bm{m}),\nabla_{\bm{m}}R^{t}(\bm{\theta},\bm{m})]$ and $\nabla J(\bm{\theta})$ are detailed in Appendix A.

For the optimization, RHOSA relies on L-BFGS-B (for Limited-memory Broyden–Fletcher–Goldfarb–Shanno with Bounds), a quasi-Newton iterative algorithm described by Zhu et al. (1997) which allows to take into account the positivity constraints of the amplitudes. In this algorithm, after an initialization $\bm{\theta}_{(0)}$, the solution is approached iteratively by

[TABLE]

where $\bm{H}^{-1}\nabla J\left(\bm{\theta},\bm{m}\right)$ is approximated with the L-BFGS formula. The iterations are repeated until one of the two following criteria is met : (1) the total number of evaluations of $J(\bm{\theta},\bm{m})$ and $\nabla J(\bm{\theta},\bm{m})$ exceeds a maximum number of iterations defined by the user; (2) the projected gradient is sufficiently small (i.e, —proj $\nabla J(\bm{\theta},\bm{m})$— / (1+—$J(\bm{\theta},\bm{m})$—) ¡ 10*-10*).

Due to its non-linearity, the least-square criterion $J(\bm{\theta},\bm{m})$ described in Eq. (8) is likely to include local minimizers. Therefore, the L-BFGS-B algorithm used by ROHSA is also likely to converge toward one of these local minima, making the solution highly dependent on the initialization $\bm{\theta}^{(I)}_{(0)}$. In order to overcome this difficulty, we designed ROHSA on an iterative multi-resolution process from coarse to fine grid (described in the following section) to automatically choose $\bm{\theta}^{(I)}_{(0)}$ and to converge towards a satisfactory local minimum.

2.4.3 ROHSA algorithm

ROHSA is based on an iterative algorithm using a multi-resolution process from coarse to fine grid presented in Algorithm 1. The number of iterations $I$ depends on the size of the fine grid $S$ and is defined by the relation $2^{I}=S$. For example, a grid of size $S^{2}$ = $256^{2}$ requires $I=8$ iterations. Each iteration is made of three steps.

Data are averaged at scale $i$ as

[TABLE]

where $\mathcal{V}_{i}$ defines the neighborhood at scale $i$, as described in Fig. 1 and $K_{i}$ is the number of positions in that neighborhood. For $i=1$, all the spatial information is compressed into a single spectrum: $\langle T_{B}\rangle_{1}=\langle T_{B}(v_{z},\bm{r})\rangle$. 2. 2.

The parameters $\bm{\theta}^{(i)}$, and $\bm{m}^{(i)}$, are estimated on that spatially averaged data version $\langle T_{B}\rangle_{i}$ by minimizing the cost function given in Eq. (8).

The minimization (line 3 of Algorithm 1) is made using L-BFGS, described in the previous section. Note that for scale $i=1$, there is no spatial information $\langle T_{B}\rangle_{1}$ and the result does not depend on the regularization. 3. 3.

Parameters $\bm{\theta}^{(i)}$ are spatially interpolated at nearest neighborhood to serve as initialization for the next scale $i+1$.

The free hyper-parameters $\lambda_{A}$, $\lambda_{\mu}$, $\lambda_{\sigma}$, $\lambda^{\prime}_{\sigma}$ remain constant during the iterations.

3 Evaluation on numerical simulation

To evaluate the performance of ROHSA, we applied it on synthetic 21 cm observations computed from a numerical simulation of thermally bi-stable turbulence flow. That allowed us to compare directly the solution given by ROHSA to the properties of the gas present in the simulation. That direct comparison with numerical reality is an essential test to evaluate the performances of a source separation algorithm like ROHSA.

3.1 Numerical simulation

To test ROHSA we used the hydrodynamical simulation of thermally bi-stable turbulence performed by Saury et al. (2014). We used their 1024N01 simulation ($1024^{3}$ pixels and a physical size of the box of 40 pc) characterized by (1) an initial density $n_{0}=0.1$ cm*-3*, (2) a large scale velocity $v_{S}=12.5$ km s*-1* and (3) a spectral weight $\zeta=0.2$. The initial density corresponds to the typical density of the WNM before condensation, the large scale velocity represents the amplitude given to the field that generates large scale turbulent motions in the box and finally, the spectral weight controls the modes of the turbulent mixing (here a majority of compressible modes). The Mach number of this simulation has been evaluated around $\mathcal{M}=0.85$ for $T>200$ K.

In order to explore the performances of ROHSA we use only a subset of this simulation. We concentrate our analysis on a $256\times 256\times 1024$ pixels region with moderate CNM fraction in order to limit the effect of HI self-absorption (see Sect. 3.2.3).

3.2 21 cm line synthetic observations

The synthetic 21 cm observations were computed using the formalism described by Miville-Deschênes & Martin (2007).

3.2.1 Distribution of velocity fluctuations

In the three dimensional spatial space, the neutral hydrogen can be described by three 3D fields : the temperature $T(\textbf{x})$, the density $\rho(\textbf{x})$, and the z-component of the turbulent velocity field $v_{z}(\textbf{x})$. Here the three dimensional spatial positions are denoted by the vector x while the two dimensional vector expressing the line-of-sight is denoted by $\bm{r}$. The z-axis is taken along the line of sight.

Information about the velocity field is inevitably lost because of the projection along z-axis. This make this description of HI a non-exhaustive one. For each position x, we assume that the velocity dispersion of a given cell is dominated by thermal motions. This is a fair approximation as the turbulent velocity dispersion at the cell size (0.04 pc) is $\sigma_{\rm turb}\sim 0.3$ km s*-1*, which is smaller than the thermal broadening everywhere in the simulation: the smallest thermal broadening for the coldest gas found in the simulation ($T=20$ K) is $\sigma_{\rm therm}=0.4$ km -1. The distribution function of the z-component of the velocity $v_{z}(\textbf{x})$ of a given cell is then given by $\phi_{v_{z}}(\textbf{x})$, a Maxwellian shifted by $v_{z}(\textbf{x})$

[TABLE]

where $\Delta(\textbf{x})=\sqrt{k_{B}T(\textbf{x})/m_{H}}$ is the thermal broadening of the 21 cm line, $m_{H}$ the hydrogen atom mass and $k_{B}$ the Boltzmann constant.

3.2.2 Brightness temperature - general case

The general case for the computation of the 21 cm brightness temperature $T_{b}(v_{z},\bm{r})$ is based on the following radiative transfer equation:

[TABLE]

where $\tau(v_{z},\bm{r},z)$ is the optical depth of the 21 cm line defined as

[TABLE]

and C=$1.82243\times 10^{18}$ cm*-2* (K km s*-1*)-1. In this representation, a gas cell at position $z^{\prime}$ absorbs emission from cell located behind along the line of sight, i.e., at $z>z^{\prime}$.

3.2.3 Optically thin limit

In the optically thin limit, in cases where the self-absorption is negligible (i.e., $\tau(v_{z},r,z)<<1$ everywhere), the 21 cm brightness temperature is proportional to the density $\rho$ :

[TABLE]

where $H$ is the depth of the cloud and $\otimes$ the spatial convolution. Here we consider the case that includes spatial smoothing by a telescope beam $B(\bm{r})$.

The integrated column density $N_{\rm HI}^{\rm thin}(\bm{r})$ and the centroid velocity $\langle v_{z}(\bm{r})\rangle$ of the 21 cm line can be obtain directly integrating $T_{B}^{\rm thin}(v_{z},\bm{r})$ along the velocity axis :

[TABLE]

and

[TABLE]

3.2.4 Synthetic observation

We computed the synthetic Position-Position-Velocity (PPV) data cube in the general case using Eq. (14). Each spectrum has an effective velocity resolution of 0.8 km s*-1* and covers $-40<v_{z}<40$ km s*-1*. We considered the beam $B(\bm{r})$ of the instrument by convolving the synthetic PPV cube with a Gaussian kernel characterized by standard deviations of two pixels along the spatial axis. Then we added a homogeneous Gaussian noise of $0.05$ K to each spectrum.

In order to mimic observation, integrated column density maps shown in the rest of the paper are computed using the optically thin limit presented Eq. (17). The integrated column density map of the synthetic PPV cube is shown in Fig. 2.

3.3 Results

ROHSA was then applied to decompose the synthetic observation computed in Sect. 3.2. In this section we discuss the choice of the free parameters of ROHSA, the global properties of the Gaussian sample and the properties of individual components identified by the decomposition. Then, the mapping of a three-phase coherent model is presented with direct comparisons to the phases extracted directly from the simulation.

3.3.1 Choosing the free parameters of ROHSA

ROHSA has 6 free parameters : the number of Gaussian components $N$, 4 hyper-parameters $\lambda_{i}$ and, the maximum number of iteration of the LBFGS algorithm.

The most important parameter is the number of Gaussian $N$. It has to be sufficiently high to ensure a complete encoding of the signal, i.e., that the residual is dominated by noise. As we will discuss in Sect. 3.3.3, a given number of Gaussian does not imply that all components are used to describe the signal along every lines of sight. Components are allowed to have an amplitude of zero at any position. This is especially relevant for components encoding cold features. Since the CNM clouds occupy a small fraction of the total volume (see also Sect. 3.3.3), we expect the associated amplitude fields to have a large fraction close to zero. This is ensured by the energy term $\|\bm{\sigma}_{n}-m_{n}\|_{2}^{2}$ which minimizes the variance of the dispersion velocity of each component. Amplitudes are brought to zero, if there is no need for a Gaussian to describe the signal at some location, instead of encoding another phase like the WNM for example. Avoiding phase mixing allows to prevent overfitting.

As with the number of Gaussian components $N$, different values of the hyper-parameters $\lambda_{\bm{a}}$, $\lambda_{\bm{\mu}}$, $\lambda_{\bm{\sigma}}$, $\lambda^{\prime}_{\bm{\sigma}}$ can be tested to obtain a satisfactory solution. If the hyper-parameters are null (i.e., no regularization) the signal could be fully encoded, but no spatial coherence will appear in the solution $\bm{\theta}^{(I)}$. On the other hand, if the hyper-parameters are too high, the solution will tend towards a solution that could be too spatially coherent, or even flat, wiping out small scale fluctuations and providing a bad fit to the data. A spatially coherent solution that describes the data well with the smallest value of $N$ is a good criterion to select the values of the hyper-parameters.

The last criterion set by the user is the maximum number of iterations of the LBFGS algorithm computing Eq. (11) (see section 2.4.2). That parameter must be large enough to ensure the convergence of the solution. The convergence of ROHSA from a numerical perspective is presented in Fig. 3 that shows the evolution of the cost function $J(\bm{\theta}(\bm{r}),\bm{m})$ for 730 iterations.

The decomposition of the synthetic observations presented in this section converged to a satisfactory solution with $N=8$, $\lambda_{\bm{a}}$ = 10000, $\lambda_{\bm{\mu}}$ = 10000, $\lambda_{\bm{\sigma}}$ = 10000, and $\lambda^{\prime}_{\bm{\sigma}}$ = 1000. As previously recommended, these values are empirically found to converge towards a noise-dominated residual and a signal that is encoded with a minimum number of Gaussian. To illustrate this, Fig. 4 shows the normalized probability distribution function of the relative difference (N${HI}$ - $\tilde{N}{HI}$) / N${HI}$ between the solution $\tilde{N}{HI}$ inferred with ROHSA and the data NHI for different number of Gaussian components $N$.

As a posterior assessment, the skewness $\mu_{3}$ of the residual is shown to quantify the quality of the encoding. As the emission is only positive and noise has a symmetric distribution centered on zero, the skewness of the residual can be used as a way to evaluate if the emission is well estimated or if the algorithm has over-fitted the data and included some noise in the solution. A positive skewness is usually an indication that some emission is left in the residual. The model has not enough freedom to encode the emission fully; more components are needed or the regularization terms should be lowered. On the other hand a negative skewness indicates that the decomposition is over-fitting the data; positive noise fluctuations are included in the model leaving more negative noise fluctuations than positive ones in the residual. This is usually an indication that the regularization coefficients ($\lambda_{i}$) should be larger to increase the spatial smoothness of the solution. A skewness of the residual close to zero is an indication of a valid solution, one that is not distorted by the regularization and that does not overfit the signal. For the case of the synthetic observations presented here, $N=8$ fully encodes the signal with a relatively low skewness —$\mu_{3}$— = 0.04.

The computation time of ROHSA depends on the maximum number of iterations, the dimension of the PPV cube and the number of Gaussians. For each step of the multi-resolution process from coarse to fine grid, the computation time used by ROHSA for this particular case (purple line) is presented in Appendix B with lower number of Gaussian for the sake of comparison. The evaluation performed here requires nearly two hours of computation time on a single CPU, which makes it difficult to explore a large range of hyper-parameters and number of Gaussian. To overcome this difficulty, A GPU implementation of ROHSA is under development.

3.3.2 Global properties of the Gaussian sample

ROHSA recovered the total emission of the synthetic observation with a relative variation of 0.3 %. An example of the Gaussian decomposition for a representative 4x4 mosaic of the simulation is shown Fig. 5. The spatial coherence of the solution can be seen over the mosaic with a smooth variation of the amplitude, the central velocity and the velocity dispersion of each Gaussian. It is already possible to distinguish in those spectra the convergence of the decomposition toward different velocity dispersions, i.e different temperatures/phases of the gas due to the energy term $\lambda^{\prime}_{\bm{\sigma}}\|\bm{\sigma}_{n}-m_{n}\|_{2}^{2}$. To have a clear view of these different components, let’s take a look at the probability distribution function $\bm{\sigma}$ weighted by the fraction of total emission of each Gaussian $\sqrt{2\pi}\bm{a}_{n}\bm{\sigma}_{n}/\sum_{\bm{r}}N_{HI}(\bm{r})$ presented in Fig. 6. This diagram shows the amount of gas in a given range of velocity dispersion (i.e indirectly a certain range of temperature). It is clear that ROHSA, in this case, converges toward a three-phase model with typical velocity dispersion close to the expected values in the CNM ($\sigma<2$ km s*-1*), LNM ($\sigma\sim 6$ km s*-1*) and WNM ($\sigma\sim 8$ km s*-1*). We will see that a similar behavior is also present in the application to a high Galactic latitude observation presented in Sect. 4. Note also that since eight Gaussian have been used by ROHSA, some phases are encoded by several components. Association of the different components to characterize a three-phase model will be presented in Sect. 3.3.4.

3.3.3 Properties of individual components

Integrated column density maps of each components obtained with ROHSA for $N=8$ Gaussian are presented in Fig. 7. For each component $G_{n}$, the mean velocity $\langle\bm{\mu}_{n}\rangle$ and mean velocity dispersion $\langle\bm{\sigma}_{n}\rangle$ averaged over the field are presented in Table 1. The surface filling factor appears to vary considerably between the eight components. The components with low values of $\langle\bm{\sigma}_{n}\rangle$ are sparsely present, while the component with the largest velocity dispersion is present everywhere on the field.

We recall that the numerical simulations used here were designed to reproduce the WNM-CNM condensation process of the HI through the thermal instability. Because of the factor 100 difference in density between the two phases, and the fact that the mass fraction in each one is about 50%, it implies that the cold phase fills only a few percents of the volume (Saury et al. 2014). This translates directly in the column density maps recovered by ROHSA; the narrow components, corresponding to colder structures, fill only a fraction of the projected field of view while the larger component is present everywhere.

The eight velocity fields and velocity dispersion fields are presented in Figures 8 and 9 respectively. At some location of the fields, when there is no need for a Gaussian to describe the signal over several pixels, the amplitude goes to zero. The corresponding velocity and velocity dispersion fields have then no reason to fluctuate. It turns out that where $a_{n}$ goes to zero, $\bm{\mu}_{n}$ and $\bm{\sigma}_{n}$ are flat. It explains the apparent variation of spatial resolution as seen for example in Fig. 8 (bottom right). We notice that the first component $G_{1}$ (top left) and the second component $G_{2}$ (top right), encoding respectively the WNM and the LNM, are defined everywhere, so $\bm{\mu}_{1}$ and $\bm{\sigma}_{1}$ also have fluctuations everywhere. Implications will be discuss in Sect. 5.

3.3.4 Mapping the three-phase neutral ISM

In order to compare the result of the decomposition with the reality given by the numerical simulations, we grouped the eight components into three fields corresponding to the WNM, the LNM and the CNM. The comparison with the numerical simulations requires that we identify ranges in temperature that demarcate these three phases. As with most numerical simulations that include the classical heating and cooling processes of the ISM (Wolfire et al. 1995), the simulation of Saury et al. (2014) we use here shows a continuum of temperature, with no clear separation and with a significant fraction of the gas present at temperatures corresponding to the thermally unstable regime (see their Figures 14 and 15). To facilitate the comparison with previous studies, we decided to use the canonical values $T_{k\,{\rm lim,CNM/LNM}}=500$ K and $T_{k\,{\rm lim,LNM/WNM}}=5000$ K (Heiles & Troland 2003b) to separate the simulation in three components. The integrated column density maps associated to each phase are computed following the methodology described in Sect. 3.2.

The comparison between the integrated column density maps recovered with ROHSA and those inferred directly from the simulation is presented Fig. 10. The intensity and the morphology of each phase is well recovered. It is nevertheless possible to see some leakages between the phases, in particular between WNM and LNM. This is due partly to the not well defined temperature thresholds used to separate the phases. It is also due to small confusions during the Gaussian decomposition where the intensity, the dispersion velocity and velocity centroid of each component in the PPV space are close to each other. In other words, for similar central velocities, the scales of fluctuations on the velocity axis characterizing each component are too close (see Figs. 5 and 6).

One way to evaluate the quality of the reconstruction is to compare the statistical properties of the cloud and inter-cloud components. Three different fields are used : the integrated column density field of the cloud medium (CNM), the integrated column density field of the inter-cloud medium (LNM+WNM), and because it is fully sampled in the plan of sky, the centroid velocity field of the inter-cloud medium. Figure 11 presents the integrated column density field and the centroid velocity field, computed using Eq. (18), of the inter-cloud medium obtained combining the LNM and the WNM inferred with ROHSA and those obtained directly from the simulation.

In order to compare the estimates from ROHSA to the ones obtained from the simulation we compute the spatial power spectrum (SPS) of each image. The SPS of the integrated column density of the cloud and inter-cloud medium are presented Fig. 12 and the SPS of the centroid velocity field of the inter-cloud medium is presented Fig.  13. In each case, the statistics recovered by ROHSA is consistent with the numerical simulation over all scales. The shape of these power spectra is interesting in itself; the inter-cloud medium is featureless with an almost constant power law, as the cloud phase is more structured, with a break at about 20 pixels, showing a typical scale linked to the condensation process. Interestingly, ROHSA is able to capture all these features very well.

4 Application on high-latitude HI gas

After validating the identification of the HI phases on numerical simulations, in this section we present the application of ROHSA on a 21 cm observation of a high Galactic latitude region.

4.1 North ecliptic pole (NEP)

To avoid the complication of low-latitude observations, where the 21 cm emission is significantly affected by velocity crowding and self-absorption, we chose to apply RHOSA on one of the high Galactic latitude fields of the GHIGLS111http://www.cita.utoronto.ca/GHIGLS/ survey (Martin et al. 2015). We chose the North Ecliptic Pole field (NEP), a 12°$\times$ 12°region, centered on $l=96^{\circ}.40$, $b=30^{\circ}.03$, observed with the Green Bank Telescope, providing a 9’.55 spatial resolution. The HI spectra have an effective velocity resolution of 0.807 km s*-1* and cover -200 ¡ v [km s*-1*] ¡ 50. The integrated column density map computed using Eq. (17) is shown Fig. 14 and a mosaic of representative emission spectra is shown in Fig. 15.

As Fig. 15 shows, high-latitude spectra of HI are more complex that the synthetic observations computed from the numerical simulations of Saury et al. (2014). This is caused by a combination of the longer line of sight in the observation (about 200 pc at $b=30^{\circ}$ compared to the 40 pc box of the simulation) and to the presence of non-local velocity components. The 21 cm emission in NEP indeed exhibits significant emission in the intermediate velocity clouds (IVC) and high velocity cloud (HVC) ranges. In this paper we do not consider the velocity channels with HVC emission. We focus on the phase separation of the local velocity cloud (LVC) and IVC components.

The LVC range, between $-20$ km s*-1* and $+20$ km s*-1* shows rather smooth emission profiles, with a narrow peak around $-3$ km s*-1* on top of a wider feature (see Fig. 15). The latter is rather smooth but when inspected in details it shows faint spectral structures on all scales along the velocity axis. The sensitivity of the GHIGLS data is such that these fluctuations of the emission profiles are not due to noise. In fact, they can be followed from one spectrum to the next quite easily. These fluctuations of the emission spectra, at scales of a few km s*-1* reveal the presence of CNM and LNM features on a range of velocities.

This field was selected because of its representative 21 cm emission for Galactic latitudes of $b\sim 30^{\circ}$. The emission features are not too complex neither too simple. In addition, a first Gaussian decomposition of NEP 21 cm data was performed by Martin et al. (2015) which provides an interesting point of comparison. Unlike ROHSA, Martin et al. (2015) used a method similar to the one described by Haud (2000), that considers only the term $\|L(v_{z},\bm{r},\bm{\theta}(\bm{r}))\|_{2}^{2}$ in the parameter optimization. We highlight some qualitative comparisons in the following sections.

4.2 Results

In order to decompose the 21 cm emission of the NEP field, we used ROHSA with $N=12$ Gaussian components and each hyper-parameter has been set to 1000. As for the previous case, these values have been chosen empirically following the same methodology as described in Sect. 3.3.1, allowing us to converge toward a noise-dominated residual with a minimum number of Gaussian. Note that the hyper-parameters values are not the same as for the first application presented in Sect. 2. The complexity of the underlying signal structure and its signal-to-noise ratio are the main causes of these differences. However, a detailed understanding of the behaviour of these hyper-parameters would require testing different values over a large number of observations. Such exploration is currently complicated by computations limitations. A GPU version of the code is under development to overcome this main limitation.

The combination of twelve Gaussian produces a solution that recovers 99 % of the total emission with spatially coherent components. The total integrated column density encoded by ROHSA and the residual between our model and the data are shown in Fig. 14 (middle and right panels).

4.2.1 Global properties of the Gaussian sample

Like for the application on the synthetic observations, the Gaussian parameters recovered for the NEP field have a strong spatial coherence; ROHSA converges towards a solution with smooth variations of the Gaussian parameters across the field. It turns out that ROHSA converges toward a multiphase model with Gaussian components of various widths, very similar to the application to numerical simulations presented in Section 3.3 but more complex due to the presence of an IVC component in the data.

To have a global view of the thermal state of the gas as a function of velocity, it is useful to look at the two-dimensional dispersion-velocity diagram $\bm{\sigma}-v$ weighted by the fraction of total emission of each Gaussian $\sqrt{2\pi}\bm{a}_{n}\bm{\sigma}_{n}/\sum_{\bm{r}}N_{HI}(\bm{r})$ shown in Fig. 16. This diagram shows isolated complexes of Gaussian components in the $\bm{\sigma}-v$ space. This is a direct result of the way the parameter optimization is done in ROHSA, with a regularization term that favours the minimum variance of $\bm{\sigma}$. Figure 16 highlights the fact that the 21 cm emission in NEP is mainly composed of negative velocity components. A clear trend is visible in the velocity range $-60<v<0$ km s*-1* with $\bm{\sigma}$ decreasing from 10 km s*-1* to 1 km s*-1* going from negative to positive velocities. It likely reflects the radiative condensation of warm intermediate velocity clouds into the local velocity component of the neutral ISM.

At this point it is interesting to compare the results of ROHSA with the Gaussian decomposition of the same data performed by Martin et al. (2015). The $\bm{\sigma}-v$ diagram of Martin et al. (2015) (see their Fig. 7) shows a continuous distribution with arches that bridge together LVC and IVC gas, an effect that is not observed in our results. Similarly to what we have done here, Martin et al. (2015) used numerical simulations to evaluate the performances of their Gaussian decomposition algorithm. Their tests revealed that such arches in the $\bm{\sigma}-v$ diagram are unphysical. They are the result of LVC and IVC gas components that overlap in velocity.

It is important to point out that both solutions are providing as good a representation of the same dataset. The significant difference between the two solutions highlights the challenge of extracting a physically meaningful representation of the data. We recall that Martin et al. (2015) used an algorithm similar to the ones used by Haud (2000); Miville-Deschênes et al. (2017a); Kalberla & Haud (2018) where the spatial coherence of the solution is not enforced through regularization terms in the cost function. In these previous studies spatial coherence is attempted by providing spatially coherent initial guesses. Then each spectrum is fitted independently and no spatial coherence in the solution is enforced. In practice, this method is rather effective for relatively sparse emission data like CO (Miville-Deschênes et al. 2017a), but in the case of the more confused 21 cm data it produces parameter maps that are more affected by small scale noise due to the degeneracy of the solution.

A Gaussian decomposition algorithm that fits each spectrum individually is easily fooled by components that overlap in velocity. In this specific case, such an algorithm would find a solution with a smaller number of components but with larger values of $\bm{\sigma}$. The main innovation in ROHSA is that it is able to cluster different phases even if they are close in velocity. The four energy terms added to the cost function $J(\bm{\theta}(\bm{r}))$ allow ROHSA to find a spatially coherent solution while avoiding the mix of components due to the high confusion present in the emission.

4.2.2 A cloud/inter-cloud medium vision of the North ecliptic pole

Integrated column density maps, centroid velocity fields and dispersion velocity fields obtained with ROHSA are respectively presented in Figs. 17, 18 and 19. Mean velocities $\langle\bm{\mu}_{n}\rangle$ and mean velocity dispersions $\langle\bm{\sigma}_{n}\rangle$ of the twelve Gaussian components Gn are presented in Table 2. In this section we focus on building a coherent cloud/inter-cloud medium vision considering the local component of the emission identified previously. Two of the four components of the local gas, $G_{7}$ and $G_{9}$, are associated with the CNM with $\langle\bm{\sigma}_{7}\rangle$ = 1.6 km s*-1* and $\langle\bm{\sigma}_{9}\rangle$ = 1.9 km s*-1*. The other ones are used to build the inter-cloud medium. Integrated column density fields and centroid velocity fields of the cloud medium and inter-cloud medium are presented Fig. 20.

As noted by Martin et al. (2015) in their two-phase decomposition of the local component, filamentary structures are observed in the narrow component. The associated velocity dispersion fields (see Fig. 19, component $G_{7}$ and $G_{9}$) shows coherent fluctuating values over a large part of the field. The core of these filamentary structures appears narrower than the envelop with velocity dispersion reaching about 0.87 km s*-1* (the spectral resolution) in their centers.

The wider component has an integrated column density field with no particular structure like filaments (see Fig. 20, top-right). Like for the numerical simulation, the sum of the wide components is likely to represent a phase that fills a large fraction of the volume, as would an inter-cloud medium. One interesting aspect of the ROHSA decomposition is that it then allows to extract the velocity field of this volume filling component (Fig. 20, bottom-right), enabling the characterization of the turbulent cascade in a mixture of lukewarm phase and warm phase.

5 Discussion

Historically a large number of studies used a Gaussian basis to model 21 cm data. Different algorithms have been developed, all of them are fitting each spectrum individually, with or without information from the neighbouring solutions to initialize the fit. To constrain further the degeneracy of the fit, solutions with the smallest number of Gaussian components have often been favoured (e.g., Lindner et al. 2015). Because of velocity blending, the solution with the smallest number of components is not necessarily the best one. In some cases, narrow features overlap in velocity, making it impossible to separate them if the environment is not considered. Usually, this confusion breaks apart a few beams away and more components can be recovered. The fundamental idea behind ROHSA is that we are trying to extract diffuse components that have column density, centroid velocity and velocity dispersion with smooth spatial variations. The optimization scheme has been designed with that concept at its core. In order to achieve this, ROHSA fits the whole data cube at once.

The application of ROHSA on both synthetic observations from numerical simulations and observational data converges naturally toward a multiphase model of the neutral ISM. The ability of ROHSA to extract the multiphase nature of the neutral ISM opens a totally new perspective on the study of the nature of the condensation process acting in the ISM. It is clear from the numerical point of view that the formation of cold clouds is the result of the condensation of the warm and diffuse gas through the thermal instability coupled with turbulence. From the observational point of view, spatial correlations between the different phases can now be made in order to quantify more precisely how the CNM emerges from this condensation process. This separation also opens the possibility to describe the properties of the very specific multiphase turbulence of the HI.

ROHSA appears to be efficient at clustering different structures in PPV space, even when there is a high level of confusion. In particular, the separation of the LVC and IVC is known to be particularly challenging as the CNM and WNM of both components overlap significantly in velocity (Martin et al. 2015). As shown in Fig. 16, ROHSA limits significantly the ”Arch effect” typical of this confusion (see Sect. 4). Globally the performance in ROHSA on high Galactic latitude regions opens a large range of possibilities regarding the study of infalling neutral clouds from the galactic halo.

We would also like to point out that no a-priori informations about the number of phases present in the neutral ISM is provided to ROHSA. The algorithm rests only on the hypothesis of the existence of components with similar line-width through the energy term $\lambda^{\prime}_{\bm{\sigma}}\|\bm{\sigma}_{n}-m_{n}\|_{2}^{2}$. In that respect, ROHSA is perfectly adapted to decompose hyper-spectral observations of any type, not only 21 cm emission.

At this time the main limitations of ROHSA are computational. Firstly, as the whole PPV cube is fitted at once, the use of ROHSA is limited to cubes that can fit in memory. Secondly, the current computation time of ROHSA is not negligible (e.g., about two hours on a single CPU for a $256\times 256\times 100$ PPV cube with 8 Gaussians). This limits the possibility to make a deep exploration of the hyper-parameters $\lambda_{i}$. In particular it would be interesting to explore various balance of the hyper-parameters for the deduced quantities ($a$, $\mu$ and $\sigma$). It is expected that the amplitude of the spatial variations of these quantities are not the same. For instance, in a multi-phase medium like the HI the density field (represented by $a$) might vary more strongly at small scale than the velocity field. This might require different values of $\lambda_{a}$ compared to the others. A GPU verion of the code is under development that would allow such an exploration.

6 Summary

In this study we have presented a new Gaussian decomposition algorithm named ROHSA. Energy terms have been added to the classical cost function to take into account simultaneously the spatial coherence of the emission and the multiphase nature of the gas. In order to identify a solution with spatially smooth parameters, the fit is performed on the whole hyper-spectral cube at once.

The performance of ROHSA has been evaluated using a synthetic 21 cm observation computed from a numerical simulation of thermally bi-stable turbulence. It was then tested on a 21 cm observation of a high Galactic latitude field observed with the GBT. The main conclusions are the following :

ROHSA is able to highlight naturally the physics of any multiphase medium without a-priori regarding the number of phase. 2. 2.

Evaluation on numerical simulation of thermally bi-stable turbulence shows that the sum of Gaussian components is a good approximation to model the multiphase nature of the neutral ISM. 3. 3.

The multiphase model inferred with ROHSA provides a spatially coherent vision of the integrated column density map, the centroid velocity field and the velocity dispersion field of each component. 4. 4.

The power spectra of the integrated column density and centroid velocity fields are well recovered with ROHSA. Statistical properties of turbulence in the multiphase neutral ISM now become accessible. 5. 5.

The decomposition of a high-latitude HI gas observation shows the wide range of applications enabled with ROHSA, for instance to study the radiative condensation of the WNM and the nature of the ISM at the disk-halo interface.

Appendix A Optimization algorithm

Terms used to compute the gradient $\nabla J(\bm{\theta},\bm{m})$ of the cost function $J(\bm{\theta},\bm{m})$ are detailed in this appendix. It involves the Jacobian of the residual $\nabla L(\bm{\theta})$ and the gradients of the regularization term $\nabla R^{t}(\bm{\theta},\bm{m})=[\nabla_{\bm{\theta}}R^{t}(\bm{\theta},\bm{m}),\nabla_{\bm{m}}R^{t}(\bm{\theta},\bm{m})]$

{ceqn}

[TABLE]

{ceqn}

[TABLE]

{ceqn}

[TABLE]

Appendix B Computation time

The computation time needed to perform the Gaussian decomposition of the synthetic PPV cube presented in Sect. 3 is described in Fig. 21. The computation time depends on the number of Gaussian component $N$, on the size of the cube (number of spectra and number of velocity channels) and on the maximum number of iteration used in the optimization.

For a given number of component $N$, the computation time scales linearly with the number of spectra and number of velocity channels. Therefore, as for each step of the multi-resolution process, from coarse to fine grid, the size of the grid is multiplied by a factor of four, the computation time also increases by a factor four at each step. Note that each step has the same maximum number of iteration (here 800) that impact also linearly the computation time.

Finally, we observed that the computation time depends non-linearly on the number of Gaussian $N$, as seen in Fig. 21.

Acknowledgements.

Part of this work was supported by Hyperstars, a project funded by the MASTODONS initiative of the CNRS mission for inter-disciplinarity and by the Programme National “Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES. 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 Lucie Riu for enlightening conversations. We thank the anonymous referee whose comments and suggestions have improved this manuscript.