Probing the structure of the gas in the Milky Way through X-ray high resolution spectroscopy

We have developed a new X-ray absorption model, called {\tt IONeq}, which computes the optical depth $\tau(E)$ simultaneously for ions of all abundant elements, assuming ionization equilibrium and taking into account turbulent broadening. We use this model to analyze the interstellar medium (ISM) absorption features in the Milky Way for a sample of 18 galactic (LMXBs) and 42 extragalactic sources (mainly Blazars). The absorbing ISM was modeled as a combination of three components/phases - neutral ($T\lesssim1\times10^{4}$ K), warm ($T\sim5\times 10^{4}$ K) and hot ($T\sim2\times10^{6}$ K). We found that the spatial distribution of both, neutral and warm components, are difficult to describe using smooth profiles due to nonuniform distribution of the column densities over the sky. For the hot phase we used a combination of a flattened disk and a halo, finding comparable column densities for both spatial components, in the order of $\sim 6-7\times10^{18}\;{\rm cm^{-2}}$, although this conclusion depends on the adopted parametrization. If the halo component has sub-solar abundance $Z$, then the column density has be scaled up by a factor $\frac{Z_\odot}{Z}$. The vertically integrated column densities of the disk components suggests the following mass fractions for these three ISM phases in the Galactic disk: neutral $\sim~89\%$, warm $\sim 8\%$ and hot $\sim 3\%$ components, respectively. The constraints on the radial distribution of the halo component of the hot component are weak.


INTRODUCTION
The interstellar medium (ISM) is an important constituent of galaxies, affecting star formation and their overall evolution. This complex environment shows multiple phases characterized by different gas temperatures (e.g. McKee & Ostriker 1977;Draine 2011). Using the high-resolution X-ray spectroscopy technique the physical conditions that characterize each phase can be analyzed. In order to carry out such an analysis bright X-ray sources with featureless spectra are required, allowing the identification of the absorption features that are imprinted in the spectra by the ISM located on the line of sight between the observer and the source. In this way, X-ray spectra provide information about basic properties of the ISM such as ionization degree, elemental abundances, column densities and temperatures.
A method commonly used to study the absorption lines due to the different ionic species observed in X-ray spectra is through fitting a Gaussian profile to individual lines. In this way the measurement of the equivalent widths (EWs) allows On the other hand, there have been attempts to model the ISM with physically motivated models, which include ionization equilibrium conditions, to study the physical conditions of the gas directly (Pinto et al. 2010(Pinto et al. , 2013Gatuzz et al. 2016;Nicastro et al. 2016a). As has been shown by Gatuzz et al. (2016), the use of physical models helps to prevent misidentification of absorption features compared with traditional phenomenological fitting methods. For in- stance, using such a physical model Nicastro et al. (2016b) concluded that the absorption line located at ∼ 22.275Å, which was previously attributed to the warm-hot intergalactic medium (WHIM), is produced by the O ii Kβ transition associated with the local ISM. In this sense, a detailed modeling of the neutral, warm and hot absorption features using the most up to date atomic data available is crucial to understand the Milky Way contribution to the atomic and ionic column densities.
In the recent years Gatuzz et al. (2013aGatuzz et al. ( ,b, 2014Gatuzz et al. ( , 2015Gatuzz et al. ( , 2016 have performed a comprehensive analysis of the galactic ISM using Chandra and XMM-Newton X-ray highresolution spectra, finding an accurate modeling of the neutral, singly and double ionized absorbers. In this paper, we present an extension of such analysis for which we have developed a new X-ray absorption model, called IONeq, which includes ionization equilibrium, velocity broadening and predicts the absorption features due to all astrophysically abundant elements. The outline of this paper is as follows. In Section 2 we describe the IONeq X-ray absorption model. In Section 3 the data selection and the reduction are summarized. Results obtained from the spectral fits are reviewed in Section 4, followed by a detailed discussion of the density profiles modeling in Section 5. In Section 6 we briefly compare our results with previous studies while possible caveats in our analysis are listed in Section 7. Finally, we summarize the results of our analysis in Section 8.

IONIZATION EQUILIBRIUM MODEL: IONEQ
We have developed a new X-ray absorption model called IONeq, which combines the information on the absorption features for astrophysically abundant elements and calculates the ionization balance for them. The latter part is based on the model developed in Churazov et al. (2001). The model assumes ionization equilibrium: where ni,j is the number density of ions j of an element i and the two terms on the r.h.s. represent creation and destruction rates of this ion. R i,j →j is the ionization/recombination rate (per ion j), leading to transition from ion j to j . The physical processes included in the creation and destruction rates are photoionization, Auger ionization, direct collisional ionization, radiative recombination and dielectronic recombination. We include the photoionization cross sections from , which are convolved with the radiation field to obtain the photoionization rates R ph i,j according to the equation where E0 is the ionization potential, J(E) is the ionizing spectrum (in units of erg s −1 cm −2 sr −1 keV −1 ) and σi,j,s(E) is the photoionization cross-section for the shell s . Ionization of electrons from inner shells can lead to ejection of one or more electrons from higher shells, which are less bounded. This process, known as Auger effect, was included by using Auger probabilities of the various possible final ionization states from Kaastra & Mewe (1993).
In the current IONeq version the radiation field J(E) is assumed to have a power law with the photon-index Γ = 2, although future versions of the model will allow the use of different ionizing radiation fields. The model includes the ionization parameter ξ that regulates the photoionization rates and is defined as n is the hydrogen density (see Tarter et al. 1969, for definitions of the ionization parameter). The range of energies used in Equation 3 is from 1 Rydberg to 1000 Rydberg, as defined in xstar (Kallman & Bautista 2001). For the collisional ionization processes, we include the ionization rates from Voronov (1997). Radiative recombination rates were taken from . In the case of dielectronic recombination we use the atomic data from Arnaud & Rothenflug (1985). Elements up to zinc are incorporated, with their respective ions. It is important to emphasize that IONeq does not solve the thermal equilibrium equations, that is the balance heating rate=cooling rate, from which one obtain the gas temperature. Instead, we treat both, the gas temperature T and the ionization parameter ξ, as independent (free) parameters of the model.
Once the ion fractions are known they can be used to compute the total optical depth τ defined as where Ai is the abundance of the i-element with respect to hydrogen, ξi,j = ni,j/ j ni,j is the ion fraction, σi,j(E) is the photoabsorption cross-section of the ion (including photoionization) and N (H) is the hydrogen column density. We use the standard solar abundances from Grevesse & Sauval (1998) and the photoabsorption cross-section used by xstar.
We included broadening effects in the expression for τ (E), assuming that the broadening consists of two components: v th and v turb . The former is due to the thermal velocities of ions where k is the Boltzmann constant, Ai is the atomic mass number and mp is the proton mass, while the latter component v turb , accounts for non-thermal effects such as turbulent motions in the gas.
The broadening associated with v turb does not depend on the mass of the ion and the convolution of the optical depth τ (E) with the corresponding Gaussian is done on-fly when fitting the data. v turb is a free parameter of the model.
The value of v th , on the other hand, depends on the gas temperature and the ion mass. We provide two versions of the IONeq model. For one version the value of v th was pre-computed for each ion by setting T in Equation (5) to the temperature at which the ion fraction ξi,j is maximal in CIE. The cross sections in Equation (4) for individual elements/ions were convolved with the corresponding Gaussians and stored in the model. This was done in order to avoid the need to convolve the optical depths for individual ions when fitting the data with the IONeq model. For the collisional ionization equilibrium this approximation is reasonable for most of ions, although for He-like and Ne-like ions, which are present over very broad range of temperatures, a more accurate approach would be to convolve the optical depth for each ion using the value of T provided as a parameter of the model. However, for the purpose of this study the main goal is not to correctly split the line broadening into thermal and turbulent components, but rather have reasonable overall broadening v 2 = v 2 th + v 2 turb . As we discuss in the next section we adjust v turb to get the right line ratios for the brightest objects in our sample and then use these values for fainter objects in our sample, assuming that the line broadening is the same for all objects. The second version of the IONeq model performs the proper optical depth convolution for each ion using the T parameter.
Finally, the total photoabsorption model takes the form where I obs (E) is the observed spectrum, Isource(E) the spectrum emitted by the source, and exp(−τ (E)) is the absorption coefficient, which is calculated by the IONeq model. In summary, the parameters of the model include the temperature Te, ionization parameter log(ξ), hydrogen column density N (H), turbulence velocity v turb and redshift z. IONeq can be downloaded 1 for use in the X-ray data analysis package xspec (Arnaud 1996).

DATA SAMPLE AND FITTING PROCEDURE
We have compiled a sample of galactic (LMXBs) and extragalactic X-ray sources from both the XMM-Newton Science Archive (XSA 2 ) and the Chandra Source Catalog (CSC 3 ) in order to study the X-ray absorption contribution from various components of the galactic ISM. In order to build our sample we have selected those sources for which we have at least 1000 counts in the oxygen edge absorption region (21-24Å). In order to get an unbiased sample we did not impose additional constraints such as, for example, significant detection of a particular line (e.g. O vii Kα). A total of 18 galactic and 42 extragalactic sources were selected. Tables 1 and 2 show the specifications of the sources, including the galactic coordinates, distances (in the  Willingale et al. (2013). N (H)-Neutral, N (H)-Warm and N (H)-Hot values are obtained from the best-fit described in Section 3.
case of LMXBs) and the 21 cm hydrogen column density measurements obtained by Willingale et al. (2013). Figure 1 shows an Aitoff projection of the location of the sources in galactic coordinates. It is clear that our final sample allows the analysis of X-ray absorption along multiple lines of sight, including regions in the galactic disk and also in the galactic halo.
A total of 257 XMM-Newton observations were reduced following the standard Scientific Analysis System (SAS, version 16.0.0) threads 4 , while 165 Chandra observations were reduced following the standard Chandra Interactive Analy-4 http://xmm.esac.esa.int/sas/current/documentation/ threads/ sis of Observations (CIAO, version 4.9) threads 5 . In the case of Chandra observations we used the findzo algorithm 6 to estimate the zero order position for the spectra. The spectral fitting was carried out with the xspec software (version 12.9.0s 7 ). Finally, we use the χ 2 statistics in combination with the Churazov et al. (1996) weighting method.
For the purpose of this study we choose to model the absorbing ISM as a combination of three components/phases -neutral (Te(neutral) 1 × 10 4 ); this component includes all neutral and molecular gas, warm (Te(warm) ∼ few × 10 4 K) and hot Te(hot) ∼ 10 6 K. We have carried out a fit in the 11-24Å wavelength region for each source listed in Tables 1 and 2 using a simple IONeq(neutral)*IONeq(warm)*IONeq(hot)*Bknpower model. The broken power-law component provides a better fit, in general, than a simple powerlaw. For every source all individual observations were fitted simultaneously and the normalization was allowed to vary between them to account for possible variations of the source flux. For each IONeq component we assume collisional ionization equilibrium (i.e. ξ = 0). The value of Te is the parameter that determines the ionic fractions, while N (H) and v turb affects the line equivalent widths.
Many objects in our sample are too faint to allow accurate determination of both Te and v turb . To mitigate this problem we decided to fix these two parameters at the bestfitting values for the brightest sources in the sample. To this end, we use Cygnus X-2 for the galactic sample and PKS 2155-304 for the extragalactic sample. Figure 2 shows the best-fitting IONeq model. For each source all observations were combined for illustrative purposes. In this plot, black data points correspond to the observations while the solid red lines represent the best-fit model for each case. The bottom panels show the fit residuals in units of σ. For the neutral component we set log(Te) = 4.0. Our best fit parameters for other two components are log(Te) = 4.7 (for the warm component) and log(Te) = 6.3 (for the hot component). Table 4 lists the ionic fractions for O and Ne corresponding to these components. In the case of unresolved and saturated lines, there is a degeneracy between the velocity broadening (i.e. Doppler parameter b = √ 2v turb ) and the column density, which can leads to under(over)estimation of the column densities (Draine 2011;Nicastro et al. 2016a).
However, if the spectrum shows absorption lines of the same ion, this degeneracy can be broken. The same approach can be extended to lines of different ions or even different elements, although in this case the results become sensitive to the assumed ionization state of the gas and the abundance ratios. This approach has been used in IONeq which fits simultaneously all lines. Later in this section we compare the broadening predicted by IONeq for the hot component with the curve of growth analysis. Figure 2 shows the spectra of the brightest sources in our sample. They contain several lines characteristic for each ISM component, allowing one to estimate the broadening. For example, for the hot component the absorption lines from both O vii Kα and O vii Kβ are dominating, providing almost model independent constraints on the broadening. The IONeq best fit parameters for the turbulence broadening are v turb = 75 +15 −55 km/s (for neutral and warm components) and v turb = 60 +28 −9 km/s (for the hot component) for Cygnus X-2. In the case of PKS 2155-304 we found v turb = 50 +12 −37 km/s (for neutral and warm components) and v turb = 110 +33 −15 km/s (for the hot component). Note that the errors quoted above reflect only the photon counting noise and do not include possible systematic errors.
The degeneracy between v turb and N (H) for saturated lines leads to large uncertainties in the v turb values. Figure 3 shows curves of growth for the hot component absorption lines (solid lines for O vii Kα, dashed lines for O vii Kβ and doted lines for Ne ix Kα). The ionic column densities have been rescaled to equivalent hydrogen column densities assuming solar abundances and CIE for Te = 2 × 10 6 K. The vertical shaded regions indicate the N (H) values obtained from the IONeq best fit for Cygnus X-2 (gray) and PKS 2155-304 (red). As the curve of growths show, if only one line is analyzed (e.g. O vii Kα) the v turb and column densities will be essentially unconstrained, however, a multiple lines analysis leads to a better determination of both parameters. We also have measured EWs of individual lines in the spectra of both sources by fitting them with Gaussians and the obtained values, including their uncertainties, are indicated as horizontal regions (solid lines for O vii Kα, dashed lines for O vii Kβ and doted lines for Ne ix Kα). We found that the column densities derived from the curve of growth method are in good agreement with the values obtained from the IONeq fit, which uses v turb and N (H) as model parameters. Table 3 shows a comparison between different v turb val- ues assumed in previous works for the neutral, warm and hot components. Due to the degeneracy between v turb and column densities in the saturated line regime it is difficult to estimate uncertainties for the v turb , specially for sources for which multiple transitions for the same ion can not be identified. In this sense, the results listed in Table 3 agree within a factor of ∼ 2, which might lead to substantial changes in the derived column densities. We decide to use our best-fit parameters for Cygnus X-2 as reference for the LMXBs and PKS 2155-304 for the extragalactic sources by fixing Te and v turb to the values obtained from the IONeq modeling. This essentially means that only the column densities of each ISM component and the broken power law continuum are the remaining free parameters of the model. It is important to note that this is the first analysis when all three ISM components are being fitted simultaneously using both Chandra and XMM-Newton high-resolution X-ray spectra.    Figure 4 shows column densities of the neutral-warmhot N (H) components from the IONeq model as function of latitude and longitude for all analyzed sources. Galactic and extragalactic sources are shown as blue filled and red empty squares, respectively. The plots show very strong correlation of the column densities with the Galactic latitude and, less prominently, with the longitude. Figure 5 shows a comparison between the N (H) values for the different ISM components. It is clear that i) the neutral component has the largest column density and ii) there is a positive correlation between column densities of different components. Such good correlation suggest that the contribution from intrinsic (to the sources) absorption is subdominant. Figure 7 shows a comparison between 21 cm survey measurements from Willingale et al. (2013), which include molecular gas, and the column densities derived from the IONeq(neutral) model. Galactic and extragalactic sources are indicated with blue filled and red empty squares, respectively. The column densities obtained for galactic sources tend to be underestimated by our model. This is expected given that 21 cm surveys provide N (H) measurements over the entire line of sight while some of the galactic sources are only few kpc away from us. For example, blue stars on . Hydrogen column densities obtained from the IONeq fits for the neutral (left panels), warm (middle panels) and hot (right panels) components as function of the galactic latitude (upper panels) and galactic longitude (bottom panels, parametrized as 360 − l for l > 180). Blue filled squares correspond to galactic sources while red empty squares correspond to extragalactic sources. Figure 7 show the N (H) value for galactic sources predicted from the density model for the neutral component described by Equation 15 and assuming a large distance (See Section 5 below). Also, it was pointed out by Gatuzz et al. (2016) that discrepancies between 21 cm survey measurements and values obtained from fitting procedure may arise from intrinsic absorption (e.g. strong stellar winds) or due to the continuum modeling. For the sake of consistency and simplicity we have used a Bknpower component to model the continuum for all sources.
[H] Figure 5. Comparison between hydrogen column densities obtained for the neutral-warm-hot components. Blue filled squares correspond to galactic sources while red empty squares correspond to extragalactic sources.

SPATIAL DISTRIBUTION OF THE GAS
For a given ionic specie i, the column density Ni depends on its density distribution ni(x) as where x is the distance along the line of sight. This relation can be used to constrain the 3D distribution of the gas. As a staring point we consider implications of the most simple possible model, consisting of a very thin disk. In this case, the expected column density for each source is expected to follow a simple law N (H) = Nz/| sin b|, where b is the galactic latitude and Nz is the normalization (vertically integrated disk column density from the disk mid-plane to infinity). To this end, we plot in the Figure 8 the values of Nz = N (H)| sin b| as a function of |b|. In the infinitely thin disk model the Nz should be the same for all sources. Figure 8 shows that for the Neutral and Warm components, the galactic sources have on average larger values of Nz. For the Hot component the trend is in opposite direction. We use this plot as a guide line to build slightly more sophisticated models below.

Neutral component
Given strong dependence of Neutral component column density on b ( Figure 4) and hints on the difference in Nz for galactic and extragalactic sources (Figure 8), a disk model with radially declining density seems appropriate. To this end we used the exponential profile model: where n0 is the normalization, Rc is the core radius, zc is the core height and (R, z) are the Galactocentric cylindrical coordinates, which are related to l, b and x as where Rsun = 8.5 kpc is the adopted distance from the galactic center to the Sun. Given that our sample lacks the objects in close vicinity of the Galactic Center, as well as objects closer to the Earth than ∼2 kpc, one can expect relatively weak constraints on the thickness of the disk zc (as long as it is small) and on the central density. There are two other issues associated with the analysis of cold component. Firstly, due to the large range of column densities of the neutral component, varying from ∼ 10 20 to ∼ 6 × 10 21 cm −2 , the signal-to-noise ratio also varies strongly, reaching ∼ 100 for galactic sources at low b. Secondly, the above model (see Equation 8) predicts smooth distribution of the column densities over the sky, while the observed 21-cm maps do show substantial variations of the intensity on different angular scales. The combination of these two factors makes the results of a straightforward fit of the model to the data dubious, since the sources with the highest column densities (or with longest exposure time) completely dominate the fit, while the entire population of sources with low column densities (extragalactic objects) or low exposure time is essentially ignored. As a result, even very sophisticated models have unacceptably large χ 2 (more than 200 per d.o.f.). To mitigate both problems, we have added quadratically to all measurements errors i) a systematic error of 30% and ii) a constant error of 6 × 10 19 cm −2 , which is the median value of the error in our sample. With these modifications the impact of galactic and extragalactic sources (as well as sources with different exposure times) is better balanced. Of course, with this approach the best-fitting values and corresponding uncertainties do not have clear statistical meaning.  Fitting all free parameters yields zc ∼ 40 pc and Rc ∼ 2.5 kpc. Fixing zc to 140 pc (see, e.g. Robin et al. 2003;Kalberla & Kerp 2009) gives comparable quality fit with Rc ∼ 6 kpc. The parameters of the latter model are given in Table 5. A comparison of the predicted and measured column densities is shown in Figure 9. According to this model the local (to the Earth) vertically integrated column density of the neutral gas is N z,N eutral ∼ 2 × 10 20 cm −2 . This number is expected to be more robust than individual parameters of the model.

Warm component
The column density of the warm component is linked to the O ii, O iii, Ne ii and Ne iii lines. The first two lines fall close to the much more prominent OI absorption features making the column density measurements less robust. In addition, for several objects the O ii Kα line overlaps with the RGS instrumental features (see, e.g. Nicastro et al. 2016a). Two such objects (Swift J1753.5-0127 and HETEJ1900.1-2455) have been excluded from the analysis of the Warm component. For other objects the data from Chandra observatory have been used. As for the Neutral component the scatter in the data points is larger than the statistical uncertainties, with a typical value of the χ 2 about 5-6 (per d.o.f.). As in Section 5.1 we used the same recipe to modify the errors: a systematic error of 30% and a medial error value (∼ 3.6 × 10 19 cm −2 ) have been quadratically added to the measurement errors. Motivated by Figures 4 and 8 we decided to use the same model as for the Neutral component, letting the normalization be the only free parameter (see Nicastro et al. 2016a, for more sophisticated models). The vertically-integrated column density of the warm component is Nz,W arm ∼ 1.8 × 10 19 cm −2 .

Hot component
For the hot component, we did not assign any additional systematic uncertainties on top of the measurements errors. The behavior of the Nz with b (see Figure 8) suggests that a single disk model cannot fit both galactic and extragalactic objects. Accordingly, the 3D distribution of the hot gas was modeled as a combination of two functions: a β-flattened profile, to account for the contribution from the galactic disk and a spherically symmetric profile to account for the possible galactic halo contribution. The total density profile for the hot component is described as where ns(x) = ns,0 1 + ( R 2 + z 2 /Rs) 2 −3βs/2 For our sample we found that the constraints on R f are weak, once R f 5 − 10 kpc. In the analysis below we fix R f to a large value (150 kpc), emulating an extended flattened component (essentially a thick disk). If we use only the flattened component, then the fit yields the disk with the vertical scale height z f ∼ 120 pc and β f ∼ 0.46 with χ 2 = 100.1 for 57 degrees of freedom (d.o.f.). Adding a spherical component with characteristic radius Rs = 13 kpc and βs = 1.1 and small βs ≈ 0. Such version of the spherical component is essentially equivalent to the assumption that the gas is distributed over a large region, so that its contribution to the column densities of galactic sources is small, while for extragalactic objects an additional (constant) value of N (H) is added. In this simplified model (with large R f and Rs), the density in the flattened component is simply which is integrated along the line of sight to get column density, while the contribution of the spherical component to the column density is where d is the distance to the source and 200 kpc is an arbitrary chosen (large) size of the spherical component. As one can see from the above expression only the total column density N (H)s,0 is an important parameter. The best-fitting parameters of this simplified model are given in Table 5. The best-fitting scale height of the flattened component is measured to be z f ≈ 350 pc, albeit with large uncertaintiesthe χ 2 increases by less than 1 with respect to the minimum when z f is in the range between 0.14 and 1.1 kpc. The vertically integrated column densities of two components of this model are comparable: ∼ 6.8 × 10 18 cm −2 and ∼ 6.3 × 10 18 cm −2 for the flattened and spherical components, respectively. Note that these are the effective column densities obtained assuming solar abundance of heavy elements. Should the abundance be different, which is plausibly the case for the spherical component, the column density of the corresponding component should be rescaled as A comparison of the predicted and measured column densities is summarized in Figure 9. While global trends are captured by the models, the data certainly show more complicated behavior than the models described in this section. We reiterate here that above models do not necessarily provide a comprehensive description of the ISM distribution in the entire Galaxy, given that our sample probes only part of the Galaxy volume. For instance, none of the sources has Figure 9. Comparison of the model-predicted and measured column densities for the Neutral, Warm and Hot components. Blue filled squares correspond to galactic sources while red empty squares correspond to extragalactic sources. The errorbars shown correspond to the modified errors for the Neutral and Warm components, and to pure statistical errors for the Hot component. The models broadly follow the trends in the data, but they clearly cannot reproduce full complexity of the measurements. line of sights passing at a 3D distance of less than ∼ 1 kpc from the Galactic center, as illustrated in Figure 6. From this point of view, our modeling of the density distribution is most sensitive to the local environment near the Sun, while the behavior of the gas within few kpc from the Galactic Center is not constrained. Often, LMXB spectra near the galactic center are so heavily affected by the neutral gas absorption that it is difficult to obtain reliable column densities for the hot component (e.g. they have few counts in the oxygen edge absorption region). Different ions, although less abundant, such as Si may be used to trace this component but such an analysis is beyond the scope of this paper.

Local ISM fractions
As is mentioned before, the vertically integrated column densities Nz at the location of the Sun are expected to be more robust than the individual parameters of the full 3D density distribution models. The values of Nz obtained for each ISM phase are given in the Table 6. The comparison of these values suggests the following mass fractions for these three ISM phases in the Galactic disk: neutral ∼ 89%, warm ∼ 8% and hot ∼ 3% components. Of course, these values reflect local (at the location of the Sun) properties. These results broadly agree with previous findings. Piontek & Ostriker (2007) performed numerical simulations including the effects of atomic heating-cooling, magnetic fields, turbulence and vertical gravity, finding that up to ∼ 91% of the ISM mass correspond to the neutral component. Yao & Wang (2006) analyzed the multiphase ISM components using Xray high-resolution spectra of the LMXB 4U 1820-303 concluding that that the neutral component accounts for ∼ 86% of the total ISM mass along this line of sight. Pinto et al. (2013), in their analysis of the ISM using 9 LMXB highresolution X-ray spectra, found that ∼ 90%, ∼ 8% and ∼ 2% of the total mass correspond to neutral, warm and hot components, respectively. Finally, Berkhuijsen & Fletcher (2008) estimated a degree of ionization of 14% by deriving gas densities from dispersion measures and H i column densities towards 393 stars and 191 pulsars. It is important to note that their sample include high-latitude sources, which can lead to a higher fraction of the ionized gas.

Neutral and warm components
For the neutral and warm components, our column densities are within 40% of those obtained by Pinto et al. (2013) for their LMXBs included in our sample. Pinto et al. (2013) used a collisionally-ionized equilibrium model included in the spex analysis package, although it is important to note that Pinto et al. (2013) assumes that the neutral component is at rest. Also, the abundances in Pinto et al. (2013) vary for each source, while in our case we are assuming solar abundances. Finally, the atomic data (i.e. photoabsorption cross-sections) used by Pinto et al. (2013) do not include Auger-damping, a crucial effect required to model accurately atomic absorption edges (García et al. 2005;Gatuzz et al. 2015).
Our results for the neutral component are within 35% of the column densities obtained by Gatuzz et al. (2016) for those LMXBs that are included in our sample. However, for the warm component we found an average of ≈ 2.5 times their column densities. In this sense, it is worth mentioning that Gatuzz et al. (2016) used the ISMabs X-ray absorption model which does not include broadening effects and the column densities are fitted independently for each ion, without considering ionization equilibrium conditions.
We obtained column densities within 25% of those obtained by Nicastro et al. (2016a), who analyzed O i Kα and O ii Kα, Kβ absorption lines and derived column densities using the curve of growth method. We have computed curve of growths for O i and O ii ions using Nicastro et al. (2016a) equivalent widths and we found that, in the case of saturated lines, the column densities can be ∼ 15% higher than those quoted by Nicastro et al. (2016a) when including Auger radiative widths.

Hot component
A similar analysis of the hot gas component associated with the Milky Way has been performed by Miller & Bregman (2013; Nicastro et al. (2016b). Main differences between our procedure and their results include (1) Miller & Bregman (2013) assumed a turbulent velocity of b = 150 km/s, consistent with the sound speed of hydrogen while our Doppler parameter was obtained directly from the spectra; (2) our sample includes both Chandra and XMM-Newton spectra leading to a better constraints on the column densities; (3) our fitting procedure consider a broad wavelength region (11-24Å), which include absorption lines from multiple ions. For example, we are modeling at same time O vii and Ne ix Kα transitions for the LMXBs spectra and if we exclude the Ne absorption edge region the column density N (H)warm can differ by a factor of ∼ 3 in some cases; (4) our density distribution model does not include offsets of the distribution from the Galaxy's center or plane, as implemented in Nicastro et al. (2016a) and (5) we are assuming solar abundance ratios and constant metallicities along all line of sights. In this sense, although our column densities can be rescaled to different solar standard abundances, we did not included abundance gradients in this analysis.

CAVEATS AND LIMITATIONS
The main caveats in our analysis are as following.
(i) Some of the extragalactic sources included in our sample are know to have intrinsic absorption. For instance, due to the presence of warm absorbers (e.g. NGC 5548, Mrk 279, Mrk 290) or ultra-fast outflow signatures (e.g. IRAS13349+2438, Parker et al. 2017). We did tests on selected sources by modeling their intrinsic absorber together with our baseline ISM model. We found that there is no strong impact on the ISM column densities obtained with the IONeq model. We intend to carry out in the near future a detailed analysis of intrinsic X-ray absorption features in extragalactic sources.
(ii) We assume fixed (solar) abundance ratios when deriving v turb . Different abundance ratios will affect the broadening velocity obtained from the fits. Also we assume that the broadening effect is independent on the line-of-sight direction.
(iii) The lack of dense coverage over sky/distances implies that our ability to place comprehensive constraints on the density distribution is limited (e.g. the absence of sources near the Galactic center).
(iv) The systematic errors are important, in particular those caused by homogeneity in the true column density distribution which is not captured by our models. This makes the estimates of the fit uncertainties for the cold and warm components useless. A large sky coverage is required to address this issue.
(v) We assume CIE, deferring the study of photoionization effects for subsequent studies.

CONCLUSIONS
High-resolution X-ray spectroscopy is a powerful technique to study the physical properties and spatial distribution of the absorbing medium in the foreground to bright X-ray sources, such as, AGNs or galactic binaries. The examples of absorbing medium range from the gas local to the source e.g., outflows from an AGN, to the ISM in the Milky Way. The absorbing properties of the absorbing medium depend on variety of parameters, including abundance of heavy elements, ionization state of the medium, turbulent velocities. To this end, we have developed a model that can account for these parameters, having in mind most common astrophysical applications. In this paper we do some basic tests of the model in application to the Milky Way ISM. In this section we briefly summarize our findings.
(i) We have developed a new X-ray absorption model, called IONeq, which computes the optical depth τ (E) taking into account ions of all astrophysically-abundant elements, assuming (collisional and photo) ionization equilibrium. The parameters of the model are the gas temperature e, ionization parameter log(ξ), hydrogen column density N (H) of the absorber, redshift (z) and turbulent broadening velocity (v turb ). The model assumes standard solar abundances from Grevesse & Sauval (1998). The model is publicly available at https://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/ models/ioneq.html.
(ii) While the model includes photoionization equilibrium, here we assume pure collisional ionization equilibrium and measure the absorption column densities in the Milky Way. The model was applied to 18 galactic (LMXBs) and 42 extragalactic sources (mainly Blazars). The velocity broadening v turb was estimated from fitting the spectra of the two brightest sources in our sample (Cygnus X-2 and PKS 2155-304) and then used for all other objects in the sample. For each source we fit a broken-power-law continuum model modified by the ISM absorption, provided by three distinct components: neutral, warm and hot. The "neutral" component, with a characteristic temperature below log(Te) 4, is primarily traced by O i Kα, the Ne i absorption edge and metallic iron absorption features. The "warm" component, with a characteristic temperature log(Te) = 4.7, is traced by O ii, O iii, Ne ii and Ne iii Kα absorption lines. A hot component, with temperature log(Te) = 6.30, can be identified by O vii, O viii, Ne ix Kα and O vii Kβ absorption lines in the spectra. For each of these components we have determined an equivalent hydrogen column density needed to explain the observed absorption features, assuming solar abundances Z . For different abundance Z, the column densities obtained from IONeq can be easily rescaled as N (H) = N (H)IONeq Z Z .
(iii) For the objects with the largest column densities, the correlation between column densities for all three ISM phases is evident in the data, as well as the correlation with the galactic latitude. These results show that (at least for these objects) the measured absorption features are due to the ISM, rather than being intrinsic to the studied objects.
(iv) The derived column densities were used to constrain the ISM density distribution models in the Galaxy. For the neutral and warm components any smooth model gives unacceptably large χ 2 , at least partly (and, possibly, predominantly) due to inhomogeneities of the true column density distribution over the sky. For that reason we decided to use a a minimalistic density distribution model n = n0 exp(−R/Rc) exp(−|z|/zc), where Rc ∼ 5.7 kpc and zc ∼ 14 pc for both components. With these parameters, the disk midplane densities at the Galactic Center are 2.1 and 0.18 cm −3 for the neutral and warm components, respectively. The local (to the disk at the Sun galactocentric radius R ∼ 8 kpc) vertically integrated column densities are N z,Neutral = 2 × 10 20 and Nz,Warm = 1.8 × 10 19 cm −2 . These quantities are expected to be more robust than other parameters of the model separately.
(v) For the hot component, the disk model is clearly insufficient and we adopted a model that includes both the contributions of a disk and a halo. The halo component is characterized by its total column density N halo = 6.3 ×10 18 cm −2 , which is comparable to the vertically integrated disk column density N z,disk = 6.8 ×10 18 cm −2 . The constraints on the radial distribution of the halo component are weak.
(vi) The comparison of the vertically integrated column densities suggests the following mass fractions for these three ISM phases in the Galactic disk: neutral ∼ 89%, warm ∼ 8% and hot ∼ 3% components, respectively. If the comparison is done using the 21 cm measurements for the neutral component, the mass fractions are similar to these values.
(vii) The absence of sources for which we can measure O and Ne lines near the galactic center requires that some assumptions about the size and thickness of the galactic disk be made. Observations focused on other ions, such as Si, may provide insights on the inner part of the galactic density profile. Figures A1 and A2 show the best-fitting models obtained for the Chandra and XMM-Newton LMXBs samples, respectively. For each source the spectra measured in individual observations were fitted simultaneously. In the plots these spectra are combined for illustrative purposes. The bestfitting models are shown with the solid red line. In the case of Chandra all spectra were taken with the Medium Energy Grating (MEG) on board of the High-Energy Transmission Grating (HETG) instrument. Figure A3 shows the bestfitting models obtained for Chandra extragalactic sources. Black points correspond to observations obtained with the Low-Energy Transmission Grating (LETG) in combination with the High-Resolution Camera (HRC). Red points correspond to observations obtained with the LETG in combination with the Advanced CCD Imaging Spectrometer (ACIS). Finally, blue points correspond to MEG observations. Figure A4 shows the best-fit obtained for XMM-Newton extragalactic sources. For all sources we obtained a good fit with χ 2 /d.o.f < 1.10 Figure A3. IONeq best fits for extragalactic Chandra LETG-HRC spectra (black points), LETG-ACIS spectra (red points) and MEG spectra (blue points). Observations are combined for illustrative purposes. Figure A4. IONeq best fits for extragalactic XMM-Newton spectra. Observations are combined for illustrative purposes.