---
abstract: 'We present a method to include the effects of gravity in the plasma physics code [Cloudy]{}. More precisely, a term is added to the desired gas pressure in order to enforce hydrostatic equilibrium, accounting for both the self-gravity of the gas and the presence of an optional external potential. As a test case, a plane-parallel model of the vertical structure of the Milky Way disk near the solar neighbourhood is considered. It is shown that the gravitational force determines the scale height of the disk, and it plays a critical role in setting its overall chemical composition. However, other variables, such as the [[shape of incident continuum and the intensity of the]{}]{} Galactic magnetic field, strongly affect the [[predicted]{}]{} structure.'
author:
- |
    Y. Ascasibar[^1] and A. I. Díaz\
    Departamento de Física Teórica, Universidad Autónoma de Madrid, 28049 Madrid, Spain
bibliography:
- 'references.bib'
date: 'Draft version 2.0 ()'
title: 'Photoionized gas in hydrostatic equilibrium: the role of gravity'
---

ISM:structure – methods:numerical – [[gravitation – radiative transfer]{}]{}

Introduction {#secIntro}
============

Many astrophysical objects, from the scale of individual stars to the hot plasma of the intracluster medium, can be approximately described as being in hydrostatic equilibrium with the underlying gravitational potential. The density and temperature of the gas are related to the enclosed mass by the condition that the net acceleration vanishes. Thus, the hydrostatic equilibrium equation reads $$\frac{1}{\rho(x)}{\frac{\ddP(x)}{\ddx}} = g(x)
\label{eqHydro_usual}$$ where $P$ denotes the total (gas, magnetic, turbulent, cosmic-ray, and radiation) pressure, and $x$ is the relevant spatial coordinate. In spherical coordinates, $x \equiv r$, and the gravitational acceleration $$g(r) = - \frac{GM(r)}{r^2} = - \frac{4\pi G}{r^2} \int_0^r x^2 \rho(x)\ {{\rm d}}x
\label{eqGrav_sphere}$$ is determined by the total mass $M(r)$ contained within radius $r$. For a symmetric plane-parallel atmosphere, the relevant coordinate is the distance to the mid plane, $z$, and the gravitational acceleration is $$g(z) = - 2\pi G\,\Sigma(z) = - 4\pi G \int_0^z \rho(x)\ {{\rm d}}x
\label{eqGrav_plane}$$ In either case, the solution to equation (\[eqHydro\_usual\]) is $$P(x) = P_0 + \Delta P_{\rm grav}(x) \equiv P_0 + \int_0^x \rho(u)\ g(u)\ {{\rm d}}u
\label{eqHydro}$$ where $P_0$ is the central pressure. There are many situations of astrophysical interest where $\Delta P_{\rm grav}(x) \ll P_0$, and the hydrostatic equilibrium condition results in an almost isobaric atmosphere. However, depending on the total mass and the characteristic temperature of the gas, the gravitational field may be strong enough to have a significant effect on the structure of the system.

One example of a photoionized gas where the gravitational acceleration is not negligible is the Galactic interstellar medium (ISM). Although interstellar matter only accounts for a relatively small fraction of the total disk mass, it plays a crucial role on the formation and evolution of the Galaxy. New stars are formed from the cold molecular phase of the ISM, and the energy input from the stellar population [[influences]{}]{} the structure and composition of the interstellar gas. The complex interplay between gas, stars, cosmic rays, and magnetic fields is responsible for the rich structure, dynamics, and observational phenomenology of our Galaxy [see e.g. @Ferriere01].

The present work describes the implementation of gravitational hydrostatic equilibrium in the plasma physics code [Cloudy]{}, last described by @Ferland+98. Numerical details are given in Section \[secCloudy\] and Appendix \[secUserGuide\]. Section \[secResults\] discusses the results obtained for the vertical structure of the Milky Way disk in the solar neighbourhood, and a brief summary is given in Section \[secConclusions\].

Numerical implementation {#secCloudy}
========================

The [Cloudy]{} photoionization code computes the ionization structure, level populations, electron temperature, and observable spectrum of a gas, given its chemical composition and the spectrum of an ionizing source. The program allows the user to set the gas density at the illuminated face and to choose from several prescriptions to specify how it varies with distance.

One of these prescriptions is [[constant pressure]{}]{}, which may refer to the gas or the total pressure $$P(x) = P_{\rm gas} + P_{\rm ram} + P_{\rm turb} + P_{\rm mag} + P_{\rm lines} + \Delta P_{\rm rad}$$ where $P_{\rm gas}=nkT$ is the thermal pressure, $P_{\rm ram}=\rho v_{\rm wind}^2$ and $P_{\rm turb}=\rho v_{\rm turb}^2/2$ arise from uniform and turbulent motion, respectively, $P_{\rm mag}=\frac{B^2}{8\pi}$ is the magnetic pressure, $P_{\rm lines}$ is the radiation pressure of trapped emission lines [@FerlandElitzur84; @ElitzurFerland86] and $\Delta P_{\rm rad} \equiv \int \rho\,a_{\rm rad}\,{{\rm d}}x$ compensates for the acceleration $a_{\rm rad}$ due to the absorption of the incident radiation. This term is analogous to $\Delta P_{\rm grav}$ and, to some extent, one could argue that the command should already be named “hydrostatic” rather than “constant pressure” even if it did not account for the effect of gravity. Since $a_{\rm rad}$ is positive (the momentum of the absorbed photons points outwards), $\Delta P_{\rm rad}$ induces a positive pressure gradient (see Figure \[figConstantP\]), whereas gravity must be balanced by a negative gradient.

We have modified the [[constant pressure]{}]{} option in version 08.00 of [Cloudy]{} to account for the gravitational term $\Delta P_{\rm grav}$. Source code is available upon request, and a brief user guide is provided in Appendix \[secUserGuide\]. The required modifications are relatively minor, and they are not expected to incur in any side effect on the functionality of the algorithm. [[They are currently being implemented in the development version of the code and will be available in the next official release[^2].]{}]{}

First, new variables are created to store the values of $\Delta P_{\rm grav}$ (initialized to 0 at the illuminated face), the symmetry of the problem (spherical or plane-parallel) and the mass (0 by default) of an additional component located at the centre.

As [Cloudy]{} moves away from the illuminated face, $\Delta P_{\rm grav}$ is updated according to $$\Delta P_{i+1} = \Delta P_i + \rho_i\, g_i\, h_i$$ where $\Delta P_i$ and $\Delta P_{i+1}$ denote the old and new values, while $\rho_i$, $g_i$, and $h_i$ are the density, gravitational acceleration, and thickness of the last zone, respectively. The gravitational acceleration is computed according to expression (\[eqGrav\_sphere\]) or (\[eqGrav\_plane\]) depending on the geometry. [Cloudy]{} keeps track of the cumulative mass $M(r)$ and the gas column density $\Sigma(z)$ in each case. Additional mass components or any external (gravitational or not) potential can be trivially accounted for, simply adding the appropriate contribution to the local acceleration $g_i$.

Finally, one has to modify the pressure loop so that the the term $\Delta P_{\rm grav}$ is added right after $\Delta P_{\rm rad}$ to balance the gravitational force. Then, the code proceeds exactly as in the original constant pressure case, trying different values of the gas density until the deviation from the desired pressure is smaller than the requested tolerance.

The Milky Way disk {#secResults}
==================

In order to test the proposed implementation of hydrostatic equilibrium, we consider a simple model of the ISM in the solar neighbourhood, assuming a plane-parallel geometry that is symmetric with respect to the Galactic mid plane, $z=0$. The gas chemical composition [@CowieSongaila86; @SavageSembach96; @Meyer+98], including dust grains [@Mathis+77], is set by the [[abundances ISM]{}]{} command, while cosmic ray heating and ionization are treated following @FerlandMushotzky84 with the instruction [[cosmic ray background]{}]{} (see Hazy[^3], the [Cloudy]{} documentation, for details).

![ Molecular, neutral, and ionized Hydrogen density as a function of height $z$ with respect to the disk mid plane (dotted, solid, and dashed lines, respectively) for different values of the mid-plane density $n_0$. []{data-label="figN"}](figs/Fig1.eps){width="50.00000%"}

![ Gas temperature (top) and pressure (bottom) along the $z$-axis for values of the mid-plane density $n_0=1.0$, $0.5$, $0.2$, and $0.1$ cm$^{-3}$. []{data-label="figP"}](figs/Fig2.eps){width="50.00000%"}

Given a mid-plane density $n_0$, the vertical structure of the disk is self-consistently computed by the program under the assumptions of ionization, thermal, and hydrostatic equilibrium. An important limitation is that the gas is described by a single density and temperature that may vary with distance, but the coexistence of different phases at the same point is not implemented [see e.g. @BoularesCox90; @SilichTenorio-Tagle98; @Ferriere01; @Cox05 for examples of multiphase models of galaxy disks in hydrostatic equilibrium]. Another simplification is that the radiation field is assumed to strike the system from one side (in our case, the mid plane), whereas in reality the stars are distributed across a scale height comparable to that of the gaseous component [e.g. @MillerCox93; @Robin+03]. With these caveats in mind, it is not surprising that our simple model does not provide a perfect match to the available observational data. It should be able, though, to give a rough estimate of the *volume-averaged* densities and scale heights of the different ISM phases. More importantly, this kind of models offers some quantitative insight on the role played by the different physical mechanisms that may be at stake. Although the present study is mostly concerned with gravity, we will also explore the influence of the spectral energy distribution of the incident radiation, as well as the presence of an external potential and a magnetic field.

The role of gravity
-------------------

To illustrate the effect of gravitational acceleration on the equilibrium structure of a generic photoionized region, we ran a first model where the incident continuum consists of the cosmic microwave background and the interstellar radiation field, interpolated from the measurements of @Black87 by means of the [[table ISM]{}]{} command.

![ Density of the different Hydrogen phases (top), electron temperature (middle) and total pressure (bottom) in the standard [[constant pressure]{}]{} case. []{data-label="figConstantP"}](figs/Fig3.eps){width="50.00000%"}

Results for different values of the central density $n_0$ are plotted in Figure \[figN\]. In all cases, there is an inner region where the gas is almost fully ionized, a thin layer of molecular and atomic Hydrogen at high densities, and an outer region of low-density gas where Hydrogen is predominantly neutral, becoming more and more ionized at large distances as density decreases. Gas temperature and pressure as a function of $z$ are shown in Figure \[figP\]. The temperature of the cold dense phase is of the order of $10-100$ K, while the diffuse gas ranges from $5000$ to $10000$ K. Gas pressure is roughly constant up to the ionization front, but it decreases rapidly through the molecular layer due to the high gas densities (above $100$ cm$^{-3}$). The pressure stabilizes again in the outer region, although the gravitational acceleration imposes a negative pressure (and density) gradient. The effect increases with distance; at some point, the cumulative gas mass and therefore the gravitational force are so high that both density and pressure decay extremely fast. Gas temperature, on the other hand, increases steadily towards the Galactic halo.

This structure does not provide a realistic description of the local ISM. On the one hand, the approximation that all the radiation comes from an infinitely thin disk at the mid plane creates a photoionized region at low $z$. Although there is indeed a thin layer of ionized gas near the Galactic plane, composed of both discrete [Hii]{} regions and diffuse gas [see e.g. @Paladini+05], its contribution near the solar neighbourhood is negligible [e.g. @TaylorCordes93; @CordesLazio_02]. On the other hand, neither the position nor the thickness of the neutral phases agree with observations of the Milky Way ISM [@BoularesCox90; @Ferriere01; @NakanishiSofue03; @NakanishiSofue06]. Nevertheless, a central ionized region may be present in several astrophysical scenarios[^4], and therefore it is interesting to test the influence of gravity in this particular example.

The role of gravity can be assessed by comparison with the standard [[constant pressure]{}]{} solution, depicted in Figure \[figConstantP\] for a mid-plane density $n_0=0.5$ cm$^{-3}$. The vertical structure of the disk is basically the same (i.e. gravity does not play a significant role) up to the ionization front. There, the temperature decreases sharply, and the density has to rise accordingly to compensate. In the standard implementation, pressure increases in order to balance the radiative acceleration ($\Delta P_{\rm rad}>0$). The effect is very significant near the ionization front, but it rapidly decreases as photons become absorbed. Although the pressure keeps increasing monotonically with $z$, the slope is quite shallow, and the physical state of the gas in the outer region (density, chemical composition, and temperature) remains approximately constant. Thus, the main effect of gravity on a plane-parallel atmosphere is to set the extent of the molecular layer and the total scale height of the system. Without gravity, the neutral (mostly molecular) layer would simply extend to infinity, as can be readily seen in Figure \[figConstantP\].

The incident radiation field
----------------------------

![ Incident continuum given by the [[table ISM]{}]{} command (dashed line), absorption of ionizing photons by the [[extinguish]{}]{} command with a Hydrogen column density of $10^{21}$ cm$^{-2}$ (solid line), and propagation through a uniform gas slab with constant density of $1$ cm$^{-3}$ and the same Hydrogen column density (dotted line). []{data-label="figContinua"}](figs/Fig4.eps){width="50.00000%"}

![ Effects of the incident continuum on the vertical structure of a model with $n_0=1$ cm$^{-3}$. The [[table ISM]{}]{} case is depicted in Figures \[figN\] and \[figP\]. The two top panels show the molecular, neutral and ionized Hydrogen densities for the incident radiation fields obtained by applying the [[extinguish]{}]{} command and the transmission through a uniform slab. The third and fourth panels show the gas temperature and the total pressure in each case. []{data-label="figExtTrans"}](figs/Fig5.eps){width="50.00000%"}

As discussed in the documentation of the [[table ISM]{}]{} command, ionizing radiation between 1 and 4 Rydberg is heavily absorbed by the neutral Hydrogen in the ISM so that few such photons exist, at least in the Galactic plane. In order to describe the actual radiation field incident on a typical region in the Galactic plane, the use of the [[extinguish]{}]{} command is recommended. However, it is also stated that this command should not be used except as a quick test, and a more physical extinction of the continuum may be accomplished by transmitting it through a model of the absorbing slab.

For the sake of simplicity, we assume a uniform density of $1$ cm$^{-3}$ and a total column density of $10^{21}$ cm$^{-2}$, appropriate for the Milky Way disk [@DickeyLockman90]. The main difference between the resulting continua, shown in Figure \[figContinua\], is that the [[extinguish]{}]{} command does not affect the energy band between 0.1 and 1 Rydberg, responsible for molecular Hydrogen photodissociation. When the incident radiation field is propagated through the ISM, many of these ultraviolet photons become absorbed by the dust component and re-emitted in the far infrared.

The details of the interstellar radiation field have a dramatic impact on the vertical structure computed by the code. Figure \[figExtTrans\] compares the results obtained for $n_0=1$ cm$^{-3}$ using the [[extinguish]{}]{} command and those obtained by propagating the [[table ISM]{}]{} field through the uniform slab before reaching the Galactic mid plane. The main difference between either method and the unextinguished case represented on the top panel of Figure \[figN\] is the absence of the inner layer of ionized gas. Central temperatures and pressures are almost two orders of magnitude lower, and, since there is no phase transition associated to any ionization front, the density of the cold neutral material is equal to $n_0$, also orders of magnitude lower than in the model without attenuation.

The absorption of the ultraviolet photodissociating radiation turns out to be very important as well. The location of the transition from the cold to the warm neutral phase may change by a large factor, but the most relevant issue is the predicted fraction of Hydrogen in molecular form. In particular, the [[extinguish]{}]{} command overestimates the photodissociation rate so much that there is virtually no molecular Hydrogen at all. Our results can only strengthen the advice given in Hazy of modelling the absorbing material rather than using the [[extinguish]{}]{} command.

Stellar mass and Galactic magnetic field
----------------------------------------

![ Total stellar mass density (solid line) and contributions of individual populations with different ages (dotted lines) according to @Robin+03. The approximation given in equation (\[eqRhoStar\]) is shown as a dashed line. []{data-label="figRhoStar"}](figs/Fig6.eps){width="50.00000%"}

![ Densities of the different Hydrogen phases (top panels), temperatures (third panel) and pressures (bottom panel) for two models including a stellar mass component of the form (\[eqRhoStar\]). In addition, one of them also includes a tangled magnetic field with an intensity $B_0=5~\mu$G at the Galactic plane. []{data-label="figRhoStar_B"}](figs/Fig7.eps){width="50.00000%"}

We have shown that photoionized regions must have a finite extent because of self-gravity. The exact value of the characteristic scale length depends, however, on many other factors. For instance, the gas may represent only a small fraction of the total mass of the system, and there can be a significant amount of non-thermal pressure support from turbulent motions, cosmic rays, and magnetic fields [@AbelFerland06; @Pellegrini+07; @Pellegrini+09].

In the case of the Milky Way, the distribution of stellar mass in the Galactic disk can be inferred by combining number counts of individual stars and dynamical arguments. We approximate the results obtained by @Robin+03 with the expression $$\rho_*(z) = \rho_0\ \left[ 1 - 0.9 \left( \frac{2z}{h_*} - \frac{z^2}{h_*^2} \right) \right]
\label{eqRhoStar}$$ where $\rho_0 = 0.042$ M$_\odot$ pc$^{-3}$ and $h_*=520$ pc. This function, plotted as a dashed line in Figure \[figRhoStar\], provides a reasonable approximation for $z<h_*$, well beyond the half-width half-maximum of the distribution[^5]. At larger distances, we take $\rho_* \approx 0$. The adopted value of the central stellar density $\rho_0$ implies the same mass as $1.7$ cm$^{-3}$ of pure Hydrogen, or $\sim1.2$ cm$^{-3}$ in the case of the ISM, where Hydrogen accounts for $\sim72$ per cent of the gas mass. The contribution of stars to the local density near the Sun is thus comparable to that of the gaseous component. The cumulative surface density $$\Sigma_*(z) = 2\rho_0\ \left[ z - 0.9 \left( \frac{z^2}{h_*} - \frac{z^3}{3h_*^2} \right) \right]
\label{eqSigmaStar}$$ increases up to $\sim 17.5$ M$_\odot$ pc$^{-2}$ at $z=h_*$, about 3 times the observed surface density of Hydrogen (in all forms) at the solar radius, and 2 times that of the ISM gas, adding up the mass in Helium and heavier elements.

When the contribution of the stellar component to the gravitational acceleration is taken into account, the pressure gradient required to maintain hydrostatic equilibrium increases, and the transition between the cold and the hot phase moves inwards by about a factor of 2 (top panel in Figure \[figRhoStar\_B\]). In general terms, the presence of an external gravitational potential will always push the gas distribution inwards, leading to a reduction of the relevant scale lengths and the total gas mass obtained for the system.

On the other hand, magnetic fields tend to counteract the action of gravity. Due to the effect of line dragging, the intensity of a tangled field increases with gas density as $B \propto n^{2/3}$ [observationally, $B \propto n^\kappa$, with $0.5<\kappa<1$, see e.g. @Crutcher99; @Henney+05], and the associated magnetic pressure varies as $P_{\rm mag} \propto B^2 \propto n^{4/3}$. This equation of state is stiffer than that of an ideal gas, and therefore a much milder change in gas density is necessary in order to achieve hydrostatic equilibrium.

The second panel of Figure \[figRhoStar\_B\] shows how the inclusion of a tangled magnetic field with mid-plane intensity $B_0=5~\mu$G [@Beck01; @HeilesCrutcher05] changes the structure of our ISM model. As can be inferred from the bottom panels, the magnetic field provides most of the pressure support against gravity. The extent of the cold phase increases by roughly an order of magnitude, reaching now hundreds of pc. In the outer, hot layer, pressure is dominated by the thermal component, and the magnetic field (for our adopted value of $B_0$) does not play a significant role.

Comparison with observations
----------------------------

![ Comparison between the model predictions and the observed densities of molecular (top), atomic (middle) and ionized (bottom) Hydrogen. The observational data are represented as dashed lines, with the grey areas delimited by dotted lines corresponding to 20 per cent error bands. Black solid lines display the results of the model with $n_0=1$ cm$^{-3}$, $B_0=5~\mu$G, and stellar density (\[eqRhoStar\]). A model with $n_0=0.75$ cm$^{-3}$, $B_0=4~\mu$G, $\kappa=0.5$, and $N({\rm H})=10^{22}$ cm$^{-2}$ is plotted in solid blue (see text for details). []{data-label="figObs"}](figs/Fig8.eps){width="50.00000%"}

Our single-phase model can only describe the average properties of the ISM. More precisely, it should be able to predict the density of molecular, atomic, and ionized Hydrogen as a function of the distance $z$ from the mid plane of the Galaxy, averaged over some scale larger than the typical size of the inhomogeneities.

Much observational work has been devoted to the determination of the average vertical structure of the Milky Way ISM near the solar neighbourhood. For the ionized component, we will adopt the distribution derived by @Reynolds91 $$\frac{ n_{\rm HII}(z) }{\rm cm^{-3}}
= 0.015 \,\exp\! \left(\! -\frac{z}{\rm 70~pc}  \right)
+ 0.0025\,\exp\! \left(\! -\frac{z}{\rm 900~pc} \right)$$ from pulsar dispersion measures, where the first term represents the contribution from discrete [Hii]{} regions near the Galactic plane, and the second term corresponds to the diffuse ionized gas. For the neutral Hydrogen, we use the formula $$n_{\rm HI, H_2}(z) = n_0 \,{\rm sech}^2 \left[ \frac{z}{z_{1/2}}\ln(1+\sqrt{2}) \right]$$ with $n_0=0.6$ cm$^{-3}$ and $z_{1/2}=125$ pc for the atomic component [@NakanishiSofue03], whereas $n_0=0.09$ cm$^{-3}$ and $z_{1/2}=93$ pc describe the molecular gas [@NakanishiSofue06]. The corresponding density distributions are shown in Figure \[figObs\] by dashed lines, and the effect of 20 per cent errors in each normalization and scale length is represented by shaded areas enclosed by dotted lines. This estimate roughly corresponds to the actual errors quoted in @Reynolds91, and it is substantially smaller than those in @NakanishiSofue06. No information about the observational errors is provided in @NakanishiSofue03, but, given the large dispersion in the Hydrogen column densities along different lines of sight [@DickeyLockman90], as well as the amplitude of the fluctuations in the azimuthally-averaged surface density profiles, 20 per cent is arguably a realistic estimate.

The results of a model with central density $n_0=1$ cm$^{-3}$, central magnetic field $B_0=5~\mu$G, and a stellar distribution given by equation (\[eqRhoStar\]) are also shown in Figure \[figObs\] as black solid lines. Although the model gives a reasonable order-of-magnitude estimate (in agreement with the observed values within a factor of two) of the scale height of the Galactic disk, as well as the densities and total masses of both atomic and molecular Hydrogen, there are some important discrepancies: first, the molecular fraction near the mid plane is severely underestimated; second, the central density is almost a factor of 2 above the observed one; and third, the amount of diffuse ionized gas is extremely low over the whole radial range, especially in the inner regions.

The first two problems are closely related. One could easily specify a lower mid plane density, but then the amount of molecular Hydrogen would decrease to unacceptably low levels. It seems that, even after transmission through a column density of $10^{21}$ cm$^{-2}$, the intensity of the photodissociating continuum is still too high for molecules to form at the observed rate. The third problem is also related to the simple prescription adopted for the incident continuum, since almost all the ionizing radiation becomes absorbed before reaching the illuminated face. In reality, young stars are spread around the Galactic plane. Adding any thin source of ionizing photons at the mid plane results in an [Hii]{} region qualitatively similar to the one observed in Figure \[figN\]; in order to reproduce the observed ionization fraction, one would need to specify a distributed radiation source, but such feature is unfortunately not currently implemented in [Cloudy]{}.

The values of the model parameters can be tuned to better match the observed total density profiles and the molecular Hydrogen fraction. Blue solid lines in Figure \[figObs\] shows the results of a model where the mid-plane density is set to $n_0=0.75$ cm$^{-3}$, the intensity of the magnetic field is $B_0=4~\mu$G at the centre and varies with gas density as $(B/B_0) = (n/n_0)^{0.5}$, and the incident continuum has been propagated through a column density of $N({\rm H})=10^{22}$ cm$^{-2}$. Nevertheless, we strongly advise caution when interpreting these results. Although the fair agreement obtained for the atomic and molecular data provides encouraging support for the validity of our implementation of gravity, a more realistic model of the ISM is clearly required in order to derive robust quantitative constraints on its physical properties. Based on our crude model, we would simply argue that reasonable values of the input parameters (in particular, central gas density and magnetic field) yield reasonable scale heights for the disk. Atomic and molecular Hydrogen densities typically agree with observations within a factor of two, but a prescription for distributed sources is necessary in order to model the ionized component.

Conclusions {#secConclusions}
===========

The present [[work]{}]{} describes an implementation of gravity in the well-tested photoionization code [Cloudy]{}. For an atmosphere in hydrostatic equilibrium, the gravitational force has to be balanced by a pressure gradient according to equation (\[eqHydro\]).

By way of example, we consider a plane-parallel model of the Galactic ISM in the solar neighbourhood to illustrate the influence of the gravitational acceleration on the physical properties of the gas. The main result, valid in the general case, is that gravity sets the scale height of the system. If self-gravity is not accounted for, there is nothing that prevents the gas (mostly molecular at large radii) from extending to infinity. When gravity is present, the pressure has to decrease in order to maintain hydrostatic equilibrium, and the gas becomes less dense. The precise structure depends on many other variables; we have investigated the effect of [[the incident radiation field, the presence of an external mass, and]{}]{} a tangled magnetic field, but other ingredients, such as turbulent motions, cosmic rays, dust grains, or variations in the atomic rate coefficients [see e.g. @Abel+08], as well as the external pressure of the intergalactic medium [@SilichTenorio-Tagle01] may also play an important role.

To summarize, the present results show that gravity (in particular, the self-gravity of the gas) should not be neglected in studies of the structure of the atomic and molecular layers of a photoionized region, especially when one is interested in their total extent.

Concerning the specific problem of the vertical structure of the Milky Way ISM, the models discussed here are able to reproduce the observed distribution of atomic and molecular Hydrogen with an accuracy of about a factor of two, but the amount of diffuse ionized gas is severely underestimated at all heights. Much more realistic results could be obtained by providing support for distributed radiation sources, as well as a treatment for the coexistence of different gas phases at the same point.

Acknowledgments {#acknowledgments .unnumbered}
===============

We would like to thank G. Ferland and the [Cloudy]{} team for developing, maintaining, and making the software publicly available. It is also a pleasure to acknowledge useful comments from G. Tenorio-Tagle and E. Pérez on a preliminary version of this manuscript. Implementation in the development version of [Cloudy]{} would not be possible without the help from G. Ferland, R. Williams, R. Porter, and P. van Hoof. [[We also thank the referee, W. Henney, for his insightful comments and suggestions, leading to major improvements in our Milky Way model.]{}]{} Financial support for this work has been provided by the Spanish *Ministerio de Educación y Ciencia* through project AYA2007-67965-C03-03.

User Guide {#secUserGuide}
==========

In order to specify a hydrostatic equilibrium configuration, the [[constant pressure]{}]{} command must be included in the [Cloudy]{} input script. In addition, a [[gravity]{}]{} command with syntax

     gravity <symmetry> [<factor> [_log]]

has been implemented. The compulsory argument [[symmetry]{}]{} must be either *spherical* or *plane-parallel* in order to tell the code whether to use equation (\[eqGrav\_sphere\]) or (\[eqGrav\_plane\]) to compute the gravitational acceleration[^6]. An optional factor can be introduced that multiplies the gas mass, mimicking the presence of an external component with identical distribution. A more flexible way to specify an external potential is provided by the command

     gravity external <mass> [<extent> <index>] [_log]

that models a distribution of mass around the centre of symmetry of the system. In the spherical case, the first number represents the mass (in solar masses) of a pointlike object located at $r=0$; in the plane-parallel case, it corresponds to the total surface density (in M$_\odot$ pc$^{-2}$) of a uniform thin sheet at the mid plane $z=0$. In both cases, the total [[mass]{}]{} $m_0$ is assumed to be distributed over an [[extent]{}]{} $x_0$ as a power law with [[index]{}]{} $\alpha$, i.e. $m(x) = m_0 (x/x_0)^\alpha$. All numbers are interpreted as linear, unless the keyword [[\_log]{}]{} is specified.

The following is an example of a typical input script used in the preparation of this work:

    # -----------------------------------------------
    #  Cloudy input file for hydrostatic equilibrium
    #                Yago Ascasibar (UAM, Fall 2009)
    # -----------------------------------------------
    #
    c -----------------------------------------------
    title - Hydrostatic equilibrium -
    c -----------------------------------------------
    c Incident continuum
    c -----------------------------------------------
    # CMB
    # table ism
    # extinguish column = 21 leak = 0
    table read "attenuated_ism.dat"
    nuF(nu) = -2.0453234 at 7.541e-05 Ryd
    c -----------------------------------------------
    c Density and chemical composition
    c -----------------------------------------------
    hden 1 linear
    abundances ism
    c -----------------------------------------------
    c Gravity
    c -----------------------------------------------
    constant pressure
    gravity plane-parallel
    gravity external Sigma = 43.680, z0 = 520, a = 1
    gravity external Sigma =-39.312, z0 = 520, a = 2
    gravity external Sigma = 13.104, z0 = 520, a = 3
    c -----------------------------------------------
    c Non-thermal components
    c -----------------------------------------------
    cosmic ray background
    magnetic field -5.3
    c -----------------------------------------------
    c Stopping conditions
    c -----------------------------------------------
    stop temperature off
    stop thickness 1e4 linear parsec
    stop eden -4.5
    iterate to convergence
    c -----------------------------------------------
    c Output files
    c -----------------------------------------------
    punch overview "overview.txt" last
    punch pressure "pressure.txt" last
    # punch transmitted continuum "trans.txt" last
    # punch hydrogen conditions "Hydrogen.txt" last
    c -----------------------------------------------
    c                         ... Paranoy@ Rulz! ;^D
    c -----------------------------------------------

The first blocks set a label to identify the run, the incident radiation field, the central Hydrogen density $n_0$ and the gas chemical composition. Commands related to hydrostatic equilibrium come next. Note, in particular, how several [[gravity external]{}]{} commands have been combined in order to specify the stellar mass distribution given by expression (\[eqSigmaStar\]). The integration continues until one of the stopping conditions is met (the disk thickness reaching 10 kpc or the electron density falling below $\sim3\times10^{-5}$ cm$^{-3}$), and it is iterated until the optical depths converge. Finally, the main physical properties of the gas (e.g. densities, temperatures) as a function of distance, as well as the different contributions to the total pressure, are printed to output files on the disk.

[^1]: E-mail: yago.ascasibar@uam.es

[^2]: The current stable version of [Cloudy]{} can be downloaded from http://wiki.nublado.org/wiki/Download

[^3]: Available at http://www.nublado.org

[^4]: Note that such an equilibrium configuration is, in any case, Rayleigh-Taylor unstable. Cold gas blobs and filamentary structures would slowly detach from the dense, neutral material and fall towards the low-density, ionized region, whereas hot gas bubbles would tend to rise buoyantly through the neutral phase.

[^5]: In fact, $\rho_*(h_*)=0.1\rho_0$

[^6]: Actually, the command parser looks for the character sequences “sphe” and “plan”, so expressions like “sphere” or “planar” are also valid.
