Issue 
A&A
Volume 626, June 2019



Article Number  A101  
Number of page(s)  19  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/201935335  
Published online  21 June 2019 
ROHSA: Regularized Optimization for HyperSpectral Analysis^{★}
Application to phase separation of 21 cm data
^{1}
AIM, CEA, CNRS, Université ParisSaclay, Université Paris Diderot,
Sorbonne Paris Cité,
91191
GifsurYvette,
France
^{2}
Institut d’Astrophysique Spatiale, CNRS UMR 8617, Université ParisSud 11,
Batiment 121,
91405
Orsay, France
email: antoine.marchal@ias.upsud.fr
^{3}
Laboratoire des Signaux et Systèmes (CNRS, CentraleSupélec, University of ParisSud), Université ParisSaclay,
91192
GifsurYvette,
France
^{4}
LIP6, Université Pierre et Marie CurieParis 6, UMR7606,
4 Place Jussieu Paris cedex 05,
75252
Paris,
France
^{5}
Instituto de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de México,
58089
Morelia,
Mexico
Received:
21
February
2019
Accepted:
29
April
2019
Context. Extracting the multiphase structure of the neutral interstellar medium is key to understanding 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 leading the atomictomolecular (HItoH_{2}) transition and the formation of molecular clouds. Up to now, the mapping of these phases out of 21 cm emission hyperspectral cubes has remained elusive mostly due to the velocity blending of individual cold structures present on a given line of sight. As a result, most of the current knowledge about the HI phases rests on a small number of absorption measurements on lines of sight crossing radio sources.
Aims. The goal of this work is to develop a new algorithm to perform separation of diffuse sources in hyperspectral 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 means to map each phase separately, but the algorithm developed here should be generic enough to extract diffuse structures in any hyperspectral cube.
Methods. We developed a new Gaussian decomposition algorithm named ROHSA based on a multiresolution process from coarse to fine grid. ROHSA uses a regularized nonlinear leastsquare criterion to take into account the spatial coherence of the emission and the multiphase nature of the gas simultaneously. 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 bistable turbulence. We apply ROHSA to a 21 cm observation of a region of high Galactic latitude from the GHIGLS survey and present our findings.
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, this test reveals that a Gaussian decomposition of HI emission is able to recover physically meaningful information about the underlying threedimensional fields (density, velocity, and temperature). The application on a real 21 cm observation of a field of high Galactic latitude produces a picture of the multiphase HI, with isolated, filamentary, and narrow (σ ~ 1−2 km s^{−1}) structures, and broader (σ ~ 4−10 km s^{−1}), diffuse, and spacefilling components. The testcase field used here contains significant intermediatevelocity 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 introduce ROHSA, a new algorithm that performs a separation of diffuse sources in hyperspectral data on the basis of a Gaussian decomposition. The algorithm makes no assumption about the nature of the sources, except that each one has a similar line width. The tests we made shows that ROHSA is well suited to decomposing complex 21 cm line emission of regions of high Galactic latitude, but its design is general enough that it could be applied to any hyperspectral data type for which a Gaussian model is relevant.
Key words: ISM: clouds / ISM: kinematics and dynamics / ISM: structure / methods: data analysis / methods: numerical / methods: observational
ROHSA is available in free access via the following web page: https://github.com/antoinemarchal/ROHSA
© A. Marchal et al. 2019
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Star formation in galaxies is strongly linked to the physical processes that govern the evolution of the interstellar medium (ISM). Stars form by gravitational collapse of dense (n > 10^{4} cm^{−3}) and cold (T ~ 10 K) structures in molecular clouds but the process that leads to the formation of these overdensities 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 radiosources, 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 cloudintercloud medium in pressure equilibrium.Very rapidly Field (1965) and 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 states: a dense cold neutral medium (CNM – T ~ 50 K, n ~ 50 cm^{−3}) immersed in a diffuse warm neutral medium (WNM – T ~ 8000 K, n ~ 0.3 cm^{−3}). This vision was later complemented by Wolfire et al. (1995, 2003) considering updated heating (dominated by the photoelectric effect on small dust grains) and cooling (dominated by CII: 158 μm, OI: 63 μm, Lα and electron recombinations onto positive charged grains) processes of the ISM.
This description of the diffuse neutral gas complemented a parallel hypothesis that emerged in the 1950s (e.g., von Weizsäcker 1951; von Hoerner 1951; Chandrasekhar & Münch 1952) and considered the ISM as a multiscale turbulent medium. In this case, the density and velocity structures are the result of a highly dynamical and outofequilibrium medium. In order to reconcile the “static/twophase” and the “turbulent” hypotheses, 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 transsonic turbulent flow, while the CNM shows a much more contrasted density structure, in accordance with the cloudintercloud 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 midway 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, these latter authors 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 an 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, 2018a). Based on a coherent modeling of emission and absorption spectra, Heiles & Troland (2003b), Murray et al. (2015, 2018a) 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 ~1% to more than 50% (Murray et al. 2018a).
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 simulations that are currently underconstrained 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; McClureGriffiths 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 multiscale nature of the HI contained in these large hyperspectral data cubes has remained elusive due to the difficulty in separating 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.
The paper is organized as follows. In Sect. 2, we describe the methodology used to develop our Gaussian decomposition algorithm. An evaluation of numerical simulation is presented in Sect. 3. In Sect. 4 we present an application to observations. 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). Very few spectra at high Galactic latitudes do not comply withthat 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 seven Gaussian components or less. 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 (σ ≤ 3 km s^{−1}; 3 < σ < 7 km s^{−1}; σ ≥ 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. These latter authors concluded that the narrow component, also seen in 21 cm absorption spectra, is 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 (σ) of 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 of these phases on the sky. Being able to separate the different phases on each line of sight would allow us to study the structure and kinematics 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, nonGaussian profiles, noise peaks, plurality of solutions. 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 velocity blending affects both the emission and absorption spectra; it is present in the data, and 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 loss ofinterest, 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, limitingthe 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 has been 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, 2018b). Indeed, key information about the nature of the HI came 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. 2018b).
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 transsonic 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 (MivilleDeschênes et al. 2003). If the HI at high Galactic latitude is indeed represented by a twophase medium with small, cold, and transsonic structures immersed in a relatively lowMachnumber 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 MivilleDeschênes et al. (2017b) have implemented Gaussian decomposition methods that use some information about their neighbors in order to favor spatially coherent solutions. Formally nevertheless, 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 surrounding area. 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 nonnegativity constraints on the amplitude. This algorithm, called ROHSA, is described below.
2.4 ROHSA
ROHSA performs a regression analysis using a regularized nonlinear leastsquare 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 the spatial coherence of the emission and the multiphase nature of the gas simultaneously. The quasiNewton algorithm, LBFGSB, used to perform the optimization is then briefly described. Finally, we formulate the algorithm performed by ROHSA based on a multiresolution process from coarse to fine grid.
2.4.1 Model
The data arethe measured brightness temperature T_{B}(v_{z}, r) at a given projected velocity v_{z} across sky coordinates r. The proposed model is a sum of N Gaussian (1)
is parametrized by with a_{n} ≥0 being the amplitude, μ_{n} the position, and σ_{n} the standard deviation 2D maps of the nth Gaussian profile across the plan of sky. The residual is (3)
The estimated parameters are defined as the minimizer of a cost function that includes the sum of the squares of the residual, (4)
where Σ is the standard deviation 2D map of the noise assumed without spatial correlation. In practice this term is estimated using a sequence of empty velocity channels of T_{B}(v_{z}, 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 smallscale spatial fluctuations of each parameter, measured by the energy at high spatial frequencies. The considered highpass filter is the secondorder differences, that is the Laplacian filtering, defined by the 2D convolution kernel, (5)
The following regularization term, itself containing energy terms, is added to the cost function given in Eq. (4), (6)
where D is a matrix performing the 2D convolution using the kernel d and λ_{a}, λ_{μ}, and λ_{σ} are hyperparameters 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 σ_{n} contains information about the gas thermodynamics, we design an additional term in the cost function to group Gaussians with similar σ_{n}. This is implemented in order to favor any solution that would produce components ascribable to each of the phases (WNM, LNM or CNM). In order to do that we add another term, , which constrains σ_{n} to be close to an unknown scalar value m_{n}. The full regularization term is then (7)
with m = (m_{1}, …, m_{N}) and a_{n} ≥ 0, ∀ n ∈ [1, …, N]. The last two terms in Eq. (7), representative of a joined constraint imposed on σ_{n}, allow 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 (8)
2.4.2 Optimization algorithm
Unlike Q(θ), each energy term proposed in Eq. (6) involves linear dependences on the parameter θ. The cost function in Eq. (8) is therefore a regularized nonlinear leastsquare criterion. The minimizer, (9)
has no closed form expression and is not directly tractable because of the complexity of the model and the size of the unknown and data. The proposed solution relies instead on an iterative optimization algorithm that uses the gradient (10)
which is tractable since it involves the residual, the Jacobian of the residual ∇L(θ), and 2D convolutions with the kernel d for D and D^{t}. The gradient ∇R^{t}(θ, m) = [∇_{θ}R^{t}(θ, m), ∇_{m}R^{t}(θ, m)] and ∇J(θ) are detailed in Appendix A.
For the optimization, RHOSA relies on LBFGSB (for Limitedmemory BroydenFletcherGoldfarbShanno with Bounds), a quasiNewton iterative algorithm described by Zhu et al. (1997) which allows for the positivity constraints of the amplitudes to be taken into account. In this algorithm, after an initialization θ_{(0)}, the solution is approached iteratively by (11)
where is approximated with the LBFGS formula. The iterations are repeated until one of the two following criteria is met: (1) the total number of evaluations of J(θ, m) and ∇J(θ, m) exceeds a maximum number of iterations defined by the user; (2) the projected gradient is sufficiently small (i.e., proj ∇J(θ, m)/(1+J(θ, m)) < 10^{−10}).
Due to its nonlinearity, the leastsquare criterion J(θ, m) described in Eq. (8) is likely to include local minimizers. Therefore, the LBFGSB algorithm used by ROHSA is also likely to converge toward one of these local minima, making the solution highly dependent on the initialization . In order to overcome this difficulty, we based the design of ROHSA on an iterative multiresolution process from coarse to fine grid (described in the following section) to automatically choose and to converge towards a satisfactory local minimum.
2.4.3 ROHSA algorithm
ROHSA is based on an iterative algorithm using a multiresolution 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.
 1.
Data are averaged at scale i as (12)
where 𝒱_{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: .
 2.
The parameters θ^{(i)} and m^{(i)} are estimated on that spatially averaged data version by minimizing the cost function given in Eq. (8).
The minimization (line 2 of Algorithm 1) is made using LBFGS, described in the previous section. We note that for scale i = 1, there is no spatial information and the result does not depend on the regularization.
 3.
Parameters θ^{(i)} are spatially interpolated at nearest neighborhood to serve as initialization for the next scale i + 1.
The free hyperparameters λ_{A}, λ_{μ}, λ_{σ}, remain constant during the iterations.
Fig. 1 Graphic visualization of neighborhoods 𝒱_{2} and 𝒱_{3} used to obtain the spatially averaged data versions and . 
3 Evaluation on numerical simulation
To evaluate the performance of ROHSA, we applied it to synthetic 21 cm observations computed from a numerical simulation of thermally bistable turbulence flow. This allowed us to directly compare 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 sourceseparation algorithm like ROHSA.
3.1 Numerical simulation
To test ROHSA we used the hydrodynamical simulation of thermally bistable 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 largescale velocity v_{S} = 12.5 km s^{−1} and (3) a spectral weight ζ = 0.2. The initial density corresponds to the typical density of the WNM before condensation, the largescale velocity represents the amplitude given to the field that generates largescale 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 to be around ℳ = 0.85 for T > 200 K.
In orderto explore the performances of ROHSA we use only a subset of this simulation. We concentrate our analysis on a region of 256 × 256 × 1024 pixels with a moderate CNM fraction in order to limit the effect of HI selfabsorption (see Sect. 3.2.3).
3.2 21 cm line synthetic observations
The synthetic 21 cm observations were computed using the formalism described by MivilleDeschênes & Martin (2007).
3.2.1 Distribution of velocity fluctuations
In the 3D spatial space, the neutral hydrogen can be described by three 3D fields: the temperature T(x), the density ρ(x), and the zcomponent of the turbulent velocity field v_{z}(x). Here the 3D spatial positions are denoted by the vector x while the 2D vector expressing the lineofsight is denoted by r. The zaxis is taken along the line of sight.
Information about the velocity field is inevitably lost because of the projection along zaxis. This makes this description of HI a nonexhaustive 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 σ_{turb} ~ 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 σ_{therm} = 0.4 km^{−1}. The distribution function of the zcomponent of the velocity v_{z} (x) of a given cell is then given by , a Maxwellian shifted by v_{z}(x), (13)
where , which is the thermal broadening of the 21 cm line, m_{H} is the hydrogen atom mass, and k_{B} the Boltzmann constant.
3.2.2 Brightness temperature: general case
The generalcase for the computation of the 21 cm brightness temperature T_{B} (v_{z}, r) is based on the following radiative transfer equation. (14)
where τ(v_{z}, r, z) is the optical depth of the 21 cm line defined as (15)
and C = 1.82243 × 10^{18} cm^{−2} (K km s^{−1})^{−1}. In this representation, a gas cell at position z′ absorbs emission from the cell located behind it along the line of sight, i.e., at z > z′.
3.2.3 Optically thin limit
In the optically thin limit, in cases where the selfabsorption is negligible (i.e., τ(v_{z}, r, z) ≪ 1 everywhere), the 21 cm brightness temperature is proportional to the density ρ: (16)
where H is the depth of the cloud and ⊗ the spatial convolution. Here we consider the case that includes spatial smoothing by a telescope beam B(r).
The integrated column density and the centroid velocity ⟨v_{z}(r)⟩ of the 21 cm line can be obtained directly by integrating along the velocity axis: (17)
Fig. 2 Integrated column density N_{HI} (optically thin approximation) of the 21 cm synthetic observation computed from the thermally bistable numerical simulation of Saury et al. (2014). 
3.2.4 Synthetic observation
We computed the synthetic positionpositionvelocity (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(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. We then added a homogeneous Gaussian noise of 0.05 K to each spectrum.
In orderto mimic observation, integrated column density maps shown in the rest of the paper are computed using the optically thin limit presented in 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. Subsequently, the mapping of a threephase 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 six free parameters : the number of Gaussian components N, four hyperparameters λ_{i} and, the maximum number of iterations of the LBFGS algorithm.
The most important parameter is N. This has to be sufficiently high to ensure a complete encoding of the signal, that is, to ensure that the residual is dominated by noise. As we discuss in Sect. 3.3.3, a given number of Gaussian components does not imply that all components are used to describe the signal along every line 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 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 for overfitting to be prevented.
As with N, different values of the hyperparameters λ_{a}, λ_{μ}, λ_{σ}, can be tested to obtain a satisfactory solution. If the hyperparameters are null (i.e., no regularization) the signal could be fully encoded, but no spatial coherence will appear in the solution θ^{(I)}. On the other hand, if the hyperparameters are too high, the solution will tend towards a solution that could be too spatially coherent, or even flat, wiping out smallscale 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 hyperparameters.
The lastcriterion set by the user is the maximum number of iterations of the LBFGS algorithm computing Eq. (11); see Sect. 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, which shows the evolution of the cost function J(θ (r), m) for 730 iterations.
The decomposition of the synthetic observations presented in this section converges to a satisfactory solution with N = 8, λ_{a} = 10 000, λ_{μ} = 10 000, λ_{σ} = 10 000, and = 1000. As previously recommended, these values are empirically found to converge towards a noisedominated residual anda signal that is encoded with a minimum number of Gaussian components. To illustrate this, Fig. 4 shows the normalized probability distribution function of the relative difference (N_{HI} – Ñ_{HI})/N_{HI} between the solution Ñ_{HI} inferred with ROHSA and the data N_{HI} for different values of N.
As an a posterior assessment, the skewness μ_{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 overfitted the data and included some noise in the solution. A positive skewness is usually an indication that some emission is left in the residual. If 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 overfitting 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 (λ_{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 μ_{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 Gaussian components. For each step of the multiresolution process from coarse to fine grid, the computation time used by ROHSA for this particular case (purple line) is presented in Appendix B with a lower number of Gaussian components 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 hyperparameters and number of Gaussian components. To overcome this difficulty, a GPU implementation of ROHSA is under development.
Fig. 3 Evolution of the cost function J(θ(r), m) as a function of the number of iterations performed by ROHSA on the synthetic observation computed in Sect. 3.2. 
Fig. 4 Normalized probability distribution function of the relative difference (N_{HI} – Ñ_{HI})/N_{HI} between the solution Ñ_{HI} inferred withttROHSA and the data N_{HI} for different numbers of Gaussian components N = [1, 3, 5, 7, 8]. The norm of the skewness μ_{3} is shown in the legend to quantify the quality of the encoding. 
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 4 × 4 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, that is, different temperatures/phases of the gas due to the energy term . To have a clear view of these different components, let us take a look at the probability distribution function σ weighted by the fraction of total emission of each Gaussian 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 threephase model with typical velocity dispersion close to the expected values in the CNM (σ < 2 km s^{−1}), LNM (σ ~ 6 km s^{−1}), and WNM (σ ~ 8 km s^{−1}). We see that a similar behavior is also present in the application to an observation of high Galactic latitude presented in Sect. 4. We note also that since eight Gaussian components have been used by ROHSA, some phases are encoded by several components. Association of the different components to characterize a threephase model is presented in Sect. 3.3.4.
3.3.3 Properties of individual components
Integrated column density maps of each component obtained with ROHSA for N = 8 Gaussian are presented in Fig. 7. For each component G_{n}, the mean velocity ⟨μ_{n} ⟩ and mean velocity dispersion ⟨σ_{n}⟩ 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 ⟨σ_{n} ⟩ 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. The factor 100 difference in density between the two phases, and the fact that the mass fraction in each one is about 50%, implies that the cold phase fills only a few percent 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 Figs. 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 then have no reason to fluctuate. It turns out that where a_{n} goes to zero, μ_{n} and σ_{n} are flat. This 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 the WNM and the LNM, respectively, are defined everywhere, meaning that μ_{1} and σ_{1} also have fluctuations everywhere. The implications of this are discussed in Sect. 5.
3.3.4 Mapping the threephase 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 Figs. 14 and 15). To facilitate the comparison with previous studies, we decided to use the canonical values T_{k lim,CNM/LNM} = 500 K and T_{k 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 poorly 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 to one another (see Figs. 5 and 6).
One way to evaluate the quality of the reconstruction is to compare the statistical properties of the cloud and intercloud components. Three different fields are used: the integrated column density field of the cloud medium (CNM), the integrated column density field of the intercloud medium (LNM + WNM), and, because it is fully sampled in the plan of sky, the centroid velocity field of the intercloud medium. Figure 11 presents the integrated column density field and the centroid velocity field, computed using Eq. (18), of the intercloud medium obtained combining the LNM and the WNM inferred with ROHSA and those obtained directly from the simulation.
In orderto compare the estimates from ROHSA to the ones obtained from the simulation, we compute the spatial power spectrum (SPS) of each image. The SPSs of the integrated column density of the cloud and intercloud medium are presented Fig. 12 and the SPS of the centroid velocity field of the intercloud medium is presented in 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 intercloud 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 scalelinked to the condensation process. Interestingly, ROHSA is able to capture all these features very well.
Fig. 5 Example of the Gaussian decomposition obtained by ROHSA for a random 4 × 4 mosaic of the synthetic observation. The original signal is shown in blue and the total brightness temperature encoded by ROHSA is shown in black. The other lines show the individual Gaussian components. The spatial coherence of the solution canbe seen over the mosaic with a smooth variation of the amplitude, the central velocity and the dispersion velocity of each component. 
Fig. 6 Probability distribution function σ weighted by the fraction of total emission of each Gaussian of the simulated field. ROHSA converges toward three distinguishable phases associated to the WNM, LNM, and CNM. 
4 Application on highlatitude 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 region with high Galactic latitude.
4.1 North ecliptic pole
To avoid the complication of lowlatitude observations, where the 21 cm emission is significantly affected by velocity crowding and selfabsorption, we chose to apply RHOSA to one of the high Galactic latitude fields of the GHIGLS^{1} survey (Martin et al. 2015). We chose the North ecliptic pole (NEP) field, a 12° × 12° region centered on l = 96°.40, b = 30°.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 in Fig. 14 and a mosaic of representative emission spectra is shown in Fig. 15.
As Fig. 15 shows, highlatitude 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° compared tothe 40 pc box of the simulation) and to the presence of nonlocal velocity components. The 21 cm emission in NEP indeed exhibits significant emission in the intermediate velocity cloud (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 and + 20 km s^{−1} shows relatively smooth emission profiles, with a narrow peak around − 3 km s^{−1} on top of a broader feature (see Fig. 15). The latter is rather smooth but when inspected in detail 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 kilometres per second 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 ~ 30°. The emission features are not particularly complex, nor are they especially 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 in the parameter optimization. We highlight some qualitative comparisons in the following sections.
Fig. 7 Integrated column density maps (left: G_{1}, G_{3}, G_{5 }, G_{7}); (right: G_{2}, G_{4}, G_{6}, G_{8}) obtained by ROHSA on the synthetic observation computed in Sect. 3.2.. Mean velocity ⟨μ_{n} ⟩ and mean velocity dispersion ⟨σ_{n}⟩ are presented in Table 1. The surface filling factor varies considerably between components, depending on their ⟨σ_{n}⟩ value. 
Fig. 8 Centroid velocity fields μ (left: μ_{1}, μ_{3}, μ_{5}, μ_{7}); (right: μ_{2}, μ_{4}, μ_{6}, μ_{8}) obtained by ROHSA using the synthetic observation computed in Sect. 3.2. 
Fig. 9 Dispersion velocity fields σ (left: σ_{1}, σ_{3}, σ_{5}, σ_{7}); (right: σ_{2}, σ_{4}, σ_{6}, σ_{8}) obtained by ROHSA using the synthetic observation computed in Sect. 3.2. 
Fig. 10 Left: integrated column density maps of the threephase model extracted by ROHSA. Right: integrated column density maps of the threephase model inferred directly from the simulation using the canonical values T_{k lim,CNM/LNM} = 500 and T_{k lim,LNM/WNM} = 5000 K. The phases WNM, LNM, and CNM are presented from top to bottom. 
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 hyperparameter 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 noisedominated residual with a minimum number of Gaussian components. We note that the hyperparameter values are not the same as for the first application presented in Sect. 2. The complexity of the underlying signal structure and its signaltonoise ratio (S/N) are the main causes of these differences. However, a detailed understanding of the behavior of these hyperparameters 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 12 Gaussian components 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).
Fig. 11 Synthetic observation of the integrated column density field (top) and the centroid velocity field (bottom) associated to the intercloud medium (WNM + LNM). Left: inferred with ROHSA. Right: inferred directly from the simulation using all cells with T_{k lim,LNM/WNM} > 500 K. For statistical comparison, the spatial power spectra of each one are shown in Figs. 12 and 13. 
Fig. 12 Spatial power spectrum of the column density. The intercloud medium (WNM + LNM) is represented by the orange dotted line (simulation) and the red line (ROHSA). The CNM is shown as a cyan dotted line (simulation) and blue line (ROHSA). 
Fig. 13 Spatial power spectrum of the centroid velocity field for the intercloud medium (WNM + LNM). The orange dotted line indicates the intercloud medium inferred directly from the simulation using all cells with T_{k lim,LNM/WNM} > 500 K (bottomright panel of Fig. 11). The red line shows the intercloud medium inferred with ROHSA (bottomleftpanel of Fig. 11). 
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 Sect. 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 twodimensional dispersionvelocity diagram σ − v weighted by the fraction of total emission of each Gaussian shown in Fig. 16. This diagram shows isolated complexes of Gaussian components in the σ − v space. This is a direct result of the way the parameter optimization is done in ROHSA, with a regularization term that favors the minimum variance of σ. Figure 16 highlights the fact that the 21 cm emission in NEP is mainly composed of negativevelocity components. A clear trend is visible in the velocity range − 60 < v < 0 km s^{−1} with σ decreasing from 10 to 1 km s^{−1} going from negative to positive velocities. This 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 σ − v diagram of Martin et al. (2015, see their Fig. 7) shows a continuous distribution with arches that bridge together the 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 σ − 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 provide 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); MivilleDeschê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. Each spectrum is then 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 (MivilleDeschênes et al. 2017a), but in the case of the more confused 21 cm data it produces parameter maps that are more affected by smallscale 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 σ. 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(θ(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/intercloud medium vision of the North ecliptic pole
Integrated column density maps, centroid velocity fields, and dispersion velocity fields obtained with ROHSA are presented in Figs. 17–19, respectively. Mean velocities ⟨μ_{n} ⟩ and mean velocity dispersions ⟨σ_{n}⟩ of the 12 Gaussian components G_{n} are presented in Table 2. In this section we focus on building a coherent cloud/intercloud 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 ⟨σ_{7}⟩ = 1.6 km s^{−1} and ⟨σ_{9}⟩ = 1.9 km s^{−1}. The other components are used to build the intercloud medium. Integrated column density fields and centroid velocity fields of the cloud medium and intercloud medium are presented in Fig. 20.
As noted by Martin et al. (2015) in their twophase 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}) 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 broader component has an integrated column density field with no particular structure like filaments (see Fig. 20, topright). Like for the numerical simulation, the sum of the broad components is likely to represent a phase that fills a large fraction of the volume, as would an intercloud medium. One interesting aspect of the ROHSA decomposition is that it then allows to extract the velocity field of this volumefilling component (Fig. 20, bottomright), enabling the characterization of the turbulent cascade in a mixture of lukewarm phase and warm phase.
Fig. 14 Left: integrated column density N_{HI} of the NEP field which is part of the GHIGLS survey. N_{HI} was computed in the optically thin approximation (see Eq. (17)). Middle: integrated column density Ñ_{HI} of NEP inferred with ROHSA. Right: residual Ñ_{HI}–N_{HI} between the integrated column density field inferred with ROHSA and the original integrated column density field computed with the data. 
Fig. 15 Example of the Gaussian decomposition obtained by ROHSA (colored line) for a random 4 × 4 mosaic of NEP. The original signal is show by the blue histogram and the total brightness temperature encoded by ROHSA isshown in black. The other lines detail the components of the Gaussian model. The spatial coherence of the solution canbe seen over the mosaic with a smooth variation of the amplitude, central velocity, and velocity dispersion of each Gaussian component. 
Fig. 16 Twodimensional probability distribution function σv weighted by the fraction of total emission of each Gaussian of NEP. The NEP is mainly composed of negative intermediate velocity components. 
Mean velocity ⟨μ_{n}⟩ and mean velocity dispersion ⟨σ_{n}⟩ of the 12 Gaussian components G_{n} inferred by ROHSA in NEP.
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 neighboring solutions to initialize the fit. To further constrain the degeneracy of the fit, solutions with the smallest number of Gaussian components have often been favored (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 a 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 an observational point of view, spatial correlations between the different phases can now be made in order to more precisely quantify 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 significantly limits the “Arch effect” typical of this confusion (see Sect. 4). Globally, the performance in ROHSA on regions of high Galactic latitude 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 information 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 . In that respect, ROHSA is perfectly adapted to decomposing hyperspectral observations of any type, not only 21 cm emission.
At this time the main limitations of ROHSA are computational. First, as the whole PPV cube is fitted at once, the use of ROHSA is limited to cubes that can fit in memory. Second, the current computation time of ROHSA is not negligible (e.g., about two hours on a single CPU for a 256 × 256 × 100 PPV cube with 8 Gaussians). This limits the possibility to make a deep exploration of the hyperparameters λ_{i}. In particular it would be interesting to explore various weights of the hyperparameters for the deduced quantities (a, μ and σ). It is expected that the amplitude of the spatial variations of these quantities are not the same. For instance, in a multiphase medium like the HI, the density field (represented by a) might vary more strongly on smaller scales than the velocity field. This might require different values of λ_{a} compared to λ_{μ}, λ_{σ} and . A GPU version of the code is under development that would allow such an exploration.
Fig. 17 Integrated column density maps (left: G_{1}, G_{4}, G_{7}, G_{10} ; middle: G_{2}, G_{5}, G_{8}, G_{11} ; right: G_{3}, G_{6}, G_{9}, G_{12}) obtained by ROHSA applied on NEP. are presented in Table 2. 
Fig. 18 Centroid velocity fields μ (left: μ_{1}, μ_{4}, μ_{7}, μ_{10} ; middle: μ_{2}, μ_{5}, μ_{8}, μ_{11} ; right: μ_{3}, μ_{6}, μ_{9}, μ_{12}) obtained by ROHSA applied on NEP. 
Fig. 19 Velocity dispersion maps σ (left: σ_{1}, σ_{4}, σ_{7}, σ_{10} ; middle: σ_{2}, σ_{5}, σ_{8}, σ_{11} ; right: σ_{3}, σ_{6}, σ_{9}, σ_{12}) obtained by ROHSA applied on NEP. 
Fig. 20 Left: CNM. Right: intercloud medium (WNM+LNM) in NEP inferred with ROHSA. Top and bottom: column density and centroid velocity fields, respectively. 
6 Summary
Here we present a new Gaussian decomposition algorithm named ROHSA. Energy terms have been added to the classical cost function to take into account the spatial coherence of the emission and the multiphase nature of the gas simultaneously. In order to identify a solution with spatially smooth parameters, the fit is performed on the whole hyperspectral cube at once.
The performance of ROHSA has been evaluated using a synthetic 21 cm observation computed from a numerical simulation of thermally bistable turbulence. It was then tested on a 21 cm observation of a field of high Galactic latitude observed with the GBT. The main conclusions are as follows.
 1.
ROHSA is able to naturally highlight the physics of any multiphase medium without a priori information regarding the number of phases.
 2.
Evaluation on numerical simulation of thermally bistable turbulence shows that the sum of Gaussian components is a good approximation to model the multiphase nature of the neutral ISM.
 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.
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.
The decomposition of a highlatitude 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 diskhalo interface.
Acknowledgements
Part of this work was supported by Hyperstars, a project funded by the MASTODONS initiative of the CNRS mission for interdisciplinarity and by the Programme National “Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU with INC/INP cofunded by CEA and CNES. This work took part under the program MilkyWayGaia of the PSI2 project funded by the IDEX ParisSaclay, ANR11IDEX000302. We gratefully acknowledge Lucie Riu for enlightening conversations. We thank the anonymous referee whose comments and suggestions have improved this manuscript.
Appendix A Optimization algorithm
Terms used to compute the gradient ∇J(θ, m) of the cost function J(θ, m) are detailedin this appendix. Incorporated are the Jacobian of the residual ∇L(θ) and the gradients of the regularization term ∇R^{t}(θ, m) = [∇_{θ}R^{t}(θ, m), ∇_{m}R^{t}(θ, m)] (A.1) (A.2) (A.3)
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. B.1. The computation time depends on the number of Gaussian components N, the size of the cube (number of spectra and number of velocity channels), and the maximum number of iterations used in the optimization.
For a given N, the computation time scales linearly with the number of spectra and the number of velocity channels. Therefore, as for each step of the multiresolution 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. We note that each step has the same maximum number of iterations (here 800) that also linearly impact the computation time.
Finally, we observed that the computation time depends nonlinearly on N, as seen in Fig. B.1.
Fig. B.1 Computation time used by one CPU to perform the Gaussian decomposition of the simulated PPV cube used in Sect. 3, for N = 1, 2, 4, 6, and 8, as function of the size grid. The maximum number of iterations in each case has been set to 800. 
References
 Audit, E., & Hennebelle, P. 2005, A&A, 433, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Begum, A., Stanimirovic, S., Goss, W. M., et al. 2010, ApJ, 725, 1779 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S., & Münch, G. 1952, ApJ, 115, 103 [NASA ADS] [CrossRef] [Google Scholar]
 Clark, B. G. 1965, ApJ, 142, 1398 [NASA ADS] [CrossRef] [Google Scholar]
 Davis, R. J. 1957, ApJ, 125, 391 [NASA ADS] [CrossRef] [Google Scholar]
 Dickey, J. M., & Lockman, F. J. 1990, ARA&A, 28, 215 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Dickey, J. M., McClureGriffiths, N. M., Gaensler, B. M., & Green, A. J. 2003, ApJ, 585, 801 [NASA ADS] [CrossRef] [Google Scholar]
 Dieter, N. H. 1964, AJ, 69, 288 [NASA ADS] [CrossRef] [Google Scholar]
 Dieter, N. H. 1965, AJ, 70, 552 [NASA ADS] [CrossRef] [Google Scholar]
 Field, G. B. 1965, ApJ, 142, 531 [NASA ADS] [CrossRef] [Google Scholar]
 Field, G. B., Goldsmith, D. W., & Habing, H. J. 1969, ApJ, 155, L149 [NASA ADS] [CrossRef] [Google Scholar]
 Haud, U. 2000, A&A, 364, 83 [NASA ADS] [Google Scholar]
 Haud, U., & Kalberla, P. M. W. 2007, A&A, 466, 555 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heeschen, D. S. 1955, ApJ, 121, 569 [NASA ADS] [CrossRef] [Google Scholar]
 Heiles, C., & Troland, T. H. 2003a, ApJS, 145, 329 [NASA ADS] [CrossRef] [Google Scholar]
 Heiles, C., & Troland, T. H. 2003b, ApJ, 586, 1067 [NASA ADS] [CrossRef] [Google Scholar]
 Hennebelle, P., & Pérault, M. 1999, A&A, 351, 309 [NASA ADS] [Google Scholar]
 Hennebelle, P., Banerjee, R., VázquezSemadeni, E., Klessen, R. S., & Audit, E. 2008, A&A, 486, L43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kalberla, P. M. W., & Haud, U. 2018, A&A, 619, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kalberla, P. M. W., Burton, W. B., Hartmann, D., et al. 2005, A&A, 440, 775 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kanekar, N., Subrahmanyan, R., Chengalur, J. N., & Safouris, V. 2003, MNRAS, 346, L57 [NASA ADS] [CrossRef] [Google Scholar]
 Koyama, H., & Inutsuka, S.i. 2002, ApJ, 564, L97 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, M.Y., Stanimirovic, S., Murray, C. E., Heiles, C., & Miller, J. 2015, ApJ, 809, 56 [NASA ADS] [CrossRef] [Google Scholar]
 Lindblad, P. O. 1966, BAN Suppl., 1, 177 [Google Scholar]
 Lindner, R. R., VeraCiro, C., Murray, C. E., et al. 2015, AJ, 149, 138 [NASA ADS] [CrossRef] [Google Scholar]
 Martin, P. G., Blagrave, K. P. M., Lockman, F. J., et al. 2015, ApJ, 809, 153 [NASA ADS] [CrossRef] [Google Scholar]
 Matthews, T. A. 1957, AJ, 62, 25 [NASA ADS] [CrossRef] [Google Scholar]
 McClureGriffiths, N. M., Pisano, D. J., Calabretta, M. R., et al. 2009, ApJS, 181, 398 [NASA ADS] [CrossRef] [Google Scholar]
 Mebold, U. 1972, A&A, 19, 13 [NASA ADS] [Google Scholar]
 MivilleDeschênes, M.A., & Martin, P. G. 2007, A&A, 469, 189 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 MivilleDeschênes, M.A., Levrier, F., & Falgarone, E. 2003, ApJ, 593, 831 [NASA ADS] [CrossRef] [Google Scholar]
 MivilleDeschênes, M.A., Murray, N., & Lee, E. J. 2017a, ApJ, 834, 57 [NASA ADS] [CrossRef] [Google Scholar]
 MivilleDeschênes, M.A., Salomé, Q., Martin, P. G., et al. 2017b, A&A, 599, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Muller, C. A. 1957, ApJ, 125, 830 [NASA ADS] [CrossRef] [Google Scholar]
 Muller, C. A. 1959, IAU Symp., 9, 360 [NASA ADS] [Google Scholar]
 Murray, C. E., Lindner, R. R., Stanimirovic, S., et al. 2014, ApJ, 781, L41 [NASA ADS] [CrossRef] [Google Scholar]
 Murray, C. E., Stanimirović, S., Goss, W. M., et al. 2015, ApJ, 804, 89 [NASA ADS] [CrossRef] [Google Scholar]
 Murray, C. E., Stanimirović, S., Goss, W. M., et al. 2018a, ApJS, 238, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Murray, C. E., Peek, J. E. G., Lee, M.Y., & Stanimirović, S. 2018b, ApJ, 862, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Ostriker, E. C., McKee, C. F., & Leroy, A. K. 2010, ApJ, 721, 975 [NASA ADS] [CrossRef] [Google Scholar]
 Peek, J. E. G., Babler, B. L., Zheng, Y., et al. 2018, ApJS, 234, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Poppel, W. G. L., Marronetti, P., & Benaglia, P. 1994, A&A, 287, 601 [NASA ADS] [Google Scholar]
 Roy, N., Kanekar, N., Braun, R., & Chengalur, J. N. 2013a, MNRAS, 436, 2352 [NASA ADS] [CrossRef] [Google Scholar]
 Roy, N., Kanekar, N., & Chengalur, J. N. 2013b, MNRAS, 436, 2366 [NASA ADS] [CrossRef] [Google Scholar]
 Saury, E., MivilleDeschênes, M.A., Hennebelle, P., Audit, E., & Schmidt, W. 2014, A&A, 567, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stanimirović,S., & Heiles, C. 2005, ApJ, 631, 371 [NASA ADS] [CrossRef] [Google Scholar]
 Stanimirović,S., Murray, C. E., Lee, M.Y., Heiles, C., & Miller, J. 2014, ApJ, 793, 132 [NASA ADS] [CrossRef] [Google Scholar]
 Stil, J. M., Taylor, A. R., Dickey, J. M., et al. 2006, AJ, 132, 1158 [NASA ADS] [CrossRef] [Google Scholar]
 Takakubo, K. 1967, BAN, 19, 125 [NASA ADS] [Google Scholar]
 Takakubo, K., & van Woerden, H. 1966, BAN, 18, 488 [NASA ADS] [Google Scholar]
 Taylor, A. R., Gibson, S. J., Peracaula, M., et al. 2003, AJ, 125, 3145 [NASA ADS] [CrossRef] [Google Scholar]
 Verschuur, G. L. 2004, AJ, 127, 394 [NASA ADS] [CrossRef] [Google Scholar]
 Verschuur, G. L., & Schmelz, J. T. 1989, AJ, 98, 267 [NASA ADS] [CrossRef] [Google Scholar]
 Verschuur, G. L., & Magnani, L. 1994, AJ, 107, 287 [NASA ADS] [CrossRef] [Google Scholar]
 von Hoerner, S. 1951, Zeitschrift fur Astrophysik, 30, 17 [NASA ADS] [Google Scholar]
 von Weizsäcker, C. F. 1951, ApJ, 114, 165 [NASA ADS] [CrossRef] [Google Scholar]
 Winkel, B., Kerp, J., Flöer, L., et al. 2016, A&A, 585, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wolfire, M. G., Hollenbach, D., McKee, C. F., Tielens, A. G. G. M., & Bakes, E. L. O. 1995, ApJ, 443, 152 [Google Scholar]
 Wolfire, M. G., McKee, C. F., Hollenbach, D., & Tielens, A. G. G. M. 2003, ApJ, 587, 278 [NASA ADS] [CrossRef] [Google Scholar]
 Zhu, C., Byrd, R. H., Lu, P., & Nocedal, J. 1997, ACM Trans. Math. Softw., 23, 550 [CrossRef] [MathSciNet] [Google Scholar]
All Tables
Mean velocity ⟨μ_{n}⟩ and mean velocity dispersion ⟨σ_{n}⟩ of the 12 Gaussian components G_{n} inferred by ROHSA in NEP.
All Figures
Fig. 1 Graphic visualization of neighborhoods 𝒱_{2} and 𝒱_{3} used to obtain the spatially averaged data versions and . 

In the text 
Fig. 2 Integrated column density N_{HI} (optically thin approximation) of the 21 cm synthetic observation computed from the thermally bistable numerical simulation of Saury et al. (2014). 

In the text 
Fig. 3 Evolution of the cost function J(θ(r), m) as a function of the number of iterations performed by ROHSA on the synthetic observation computed in Sect. 3.2. 

In the text 
Fig. 4 Normalized probability distribution function of the relative difference (N_{HI} – Ñ_{HI})/N_{HI} between the solution Ñ_{HI} inferred withttROHSA and the data N_{HI} for different numbers of Gaussian components N = [1, 3, 5, 7, 8]. The norm of the skewness μ_{3} is shown in the legend to quantify the quality of the encoding. 

In the text 
Fig. 5 Example of the Gaussian decomposition obtained by ROHSA for a random 4 × 4 mosaic of the synthetic observation. The original signal is shown in blue and the total brightness temperature encoded by ROHSA is shown in black. The other lines show the individual Gaussian components. The spatial coherence of the solution canbe seen over the mosaic with a smooth variation of the amplitude, the central velocity and the dispersion velocity of each component. 

In the text 
Fig. 6 Probability distribution function σ weighted by the fraction of total emission of each Gaussian of the simulated field. ROHSA converges toward three distinguishable phases associated to the WNM, LNM, and CNM. 

In the text 
Fig. 7 Integrated column density maps (left: G_{1}, G_{3}, G_{5 }, G_{7}); (right: G_{2}, G_{4}, G_{6}, G_{8}) obtained by ROHSA on the synthetic observation computed in Sect. 3.2.. Mean velocity ⟨μ_{n} ⟩ and mean velocity dispersion ⟨σ_{n}⟩ are presented in Table 1. The surface filling factor varies considerably between components, depending on their ⟨σ_{n}⟩ value. 

In the text 
Fig. 8 Centroid velocity fields μ (left: μ_{1}, μ_{3}, μ_{5}, μ_{7}); (right: μ_{2}, μ_{4}, μ_{6}, μ_{8}) obtained by ROHSA using the synthetic observation computed in Sect. 3.2. 

In the text 
Fig. 9 Dispersion velocity fields σ (left: σ_{1}, σ_{3}, σ_{5}, σ_{7}); (right: σ_{2}, σ_{4}, σ_{6}, σ_{8}) obtained by ROHSA using the synthetic observation computed in Sect. 3.2. 

In the text 
Fig. 10 Left: integrated column density maps of the threephase model extracted by ROHSA. Right: integrated column density maps of the threephase model inferred directly from the simulation using the canonical values T_{k lim,CNM/LNM} = 500 and T_{k lim,LNM/WNM} = 5000 K. The phases WNM, LNM, and CNM are presented from top to bottom. 

In the text 
Fig. 11 Synthetic observation of the integrated column density field (top) and the centroid velocity field (bottom) associated to the intercloud medium (WNM + LNM). Left: inferred with ROHSA. Right: inferred directly from the simulation using all cells with T_{k lim,LNM/WNM} > 500 K. For statistical comparison, the spatial power spectra of each one are shown in Figs. 12 and 13. 

In the text 
Fig. 12 Spatial power spectrum of the column density. The intercloud medium (WNM + LNM) is represented by the orange dotted line (simulation) and the red line (ROHSA). The CNM is shown as a cyan dotted line (simulation) and blue line (ROHSA). 

In the text 
Fig. 13 Spatial power spectrum of the centroid velocity field for the intercloud medium (WNM + LNM). The orange dotted line indicates the intercloud medium inferred directly from the simulation using all cells with T_{k lim,LNM/WNM} > 500 K (bottomright panel of Fig. 11). The red line shows the intercloud medium inferred with ROHSA (bottomleftpanel of Fig. 11). 

In the text 
Fig. 14 Left: integrated column density N_{HI} of the NEP field which is part of the GHIGLS survey. N_{HI} was computed in the optically thin approximation (see Eq. (17)). Middle: integrated column density Ñ_{HI} of NEP inferred with ROHSA. Right: residual Ñ_{HI}–N_{HI} between the integrated column density field inferred with ROHSA and the original integrated column density field computed with the data. 

In the text 
Fig. 15 Example of the Gaussian decomposition obtained by ROHSA (colored line) for a random 4 × 4 mosaic of NEP. The original signal is show by the blue histogram and the total brightness temperature encoded by ROHSA isshown in black. The other lines detail the components of the Gaussian model. The spatial coherence of the solution canbe seen over the mosaic with a smooth variation of the amplitude, central velocity, and velocity dispersion of each Gaussian component. 

In the text 
Fig. 16 Twodimensional probability distribution function σv weighted by the fraction of total emission of each Gaussian of NEP. The NEP is mainly composed of negative intermediate velocity components. 

In the text 
Fig. 17 Integrated column density maps (left: G_{1}, G_{4}, G_{7}, G_{10} ; middle: G_{2}, G_{5}, G_{8}, G_{11} ; right: G_{3}, G_{6}, G_{9}, G_{12}) obtained by ROHSA applied on NEP. are presented in Table 2. 

In the text 
Fig. 18 Centroid velocity fields μ (left: μ_{1}, μ_{4}, μ_{7}, μ_{10} ; middle: μ_{2}, μ_{5}, μ_{8}, μ_{11} ; right: μ_{3}, μ_{6}, μ_{9}, μ_{12}) obtained by ROHSA applied on NEP. 

In the text 
Fig. 19 Velocity dispersion maps σ (left: σ_{1}, σ_{4}, σ_{7}, σ_{10} ; middle: σ_{2}, σ_{5}, σ_{8}, σ_{11} ; right: σ_{3}, σ_{6}, σ_{9}, σ_{12}) obtained by ROHSA applied on NEP. 

In the text 
Fig. 20 Left: CNM. Right: intercloud medium (WNM+LNM) in NEP inferred with ROHSA. Top and bottom: column density and centroid velocity fields, respectively. 

In the text 
Fig. B.1 Computation time used by one CPU to perform the Gaussian decomposition of the simulated PPV cube used in Sect. 3, for N = 1, 2, 4, 6, and 8, as function of the size grid. The maximum number of iterations in each case has been set to 800. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.