% revised version sept 1995

%\documentstyle[12pt,aasms]{article}
\documentstyle[12pt,aaspp]{article}
% remove \tightenlines before submitting
\tightenlines

\def\nh{$n_{\rm H}$}
\def\nhh{$n_{{\rm H}_2}$}
\def\h2{H$_2$}
\def\hh2{H~+~H$_2$}
\def\molh2{H$_2$~+~H$_2$}
\def\vj{$(v,j)$}
\def\etal{et al.}
\def\delj{$\Delta j$}
\def\delv{$\Delta v$}
\def\jmax{$j_{max}$}
\def\cms{cm$^3$~s$^{-1}$}
\def\cm3{cm$^{-3}$}
\def\ecps{erg~cm$^3$~s$^{-1}$}
\def\ecms{erg~cm$^{-3}$~s$^{-1}$}

% macros specifically for the table
\def\phss{\phantom{1234$-$}}
\def\phaa{\phantom{123$-$}}
\def\phbb{\phantom{123}$-$}
\def\phcc{\phantom{1}}
\def\arm{rms \tablenotemark {c}}
\def\afx{max \tablenotemark {c}}
\def\afn{min \tablenotemark {c}}
\def\pha{\phantom{$-$}}
\def\phb{$-$}
\def\phc{$+$}
\def\phd{$-$}
\def\emp{\phantom{12345}$\cdots$}
\def\sla{\\}

\begin{document}

\title{MASTER EQUATION STUDIES OF COLLISIONAL EXCITATION AND
DISSOCIATION OF H$_2$ MOLECULES BY H ATOMS} 

\author{P. G. Martin}

\affil{Canadian Institute for Theoretical Astrophysics, University of
Toronto, Toronto, ON, Canada M5S~1A7.  pgmartin@cita.utoronto.ca}

\author{D. H. Schwarz}

\affil{Institute of Astronomy, The Observatories, Madingley Road,
Cambridge, England CB3~OHA.  schwarz@ast.cam.ac.uk}

\author{and}

\author{M.\ E.\ Mandy}

\affil{Programme in Chemistry, University of Northern British
Columbia, Prince George, BC, Canada V2N~4Z9. mandy@unbc.edu}

\begin{abstract}

The master equation for collisional excitation of \h2\ by H atoms is
solved over a temperature range of $T$~=~450 to 45,000~K and a density
range of \nh~=~10$^{-3}$ to 10$^9$~\cm3.  Radiative transitions due to
quadrupole emission and dissociative tunneling are included
explicitly. The resulting non-equilibrium steady-state populations of
the \vj\ states of \h2\ are determined and then used to calculate the
rate coefficients for ortho-para interconversion, collisional
dissociation, dissociative tunneling, and cooling. The density and
temperature dependences of these rates are discussed and parameterized.
Comparisons with previous estimates are made.

\bigskip\noindent {\it Subject headings:}

ISM: molecules -- molecular processes -- molecular data -- line:
formation

\end{abstract}

\section{INTRODUCTION}
\label{sec-int}

Evaluation of several fundamental processes involving \h2\ requires a
detailed knowledge of the populations of the vibration-rotation [\vj]
states of the molecule; prominent among these are prediction of the
quadrupole emission spectrum, radiative and dissociative cooling in the
relaxation of an interstellar shock (e.g., Hollenbach \& McKee 1989),
and heating by ``quenching'' of \h2\ that is fluorescently excited in
photodissociation regions (e.g., Sternberg \& Dalgarno 1989).  Since the
interstellar medium is not generally in local thermodynamic equilibrium,
the relative populations of the \vj\ states of molecular \h2\ cannot be
described adequately with the Boltzmann distribution.  Instead, the
competing processes which can affect the population such as radiative
emission and absorption, and collisional excitation, deexcitation, and
dissociation must be treated explicitly.  The equation which describes
the evolution of state populations is often called the master equation
(\S~2).

In this paper we use a new set of rate coefficients for collisions of
\h2\ with H, describing both state-to-state transitions and collision
induced dissociation (CID) (Mandy \& Martin 1993, 1995).  The
significant differences in these rate coefficients with respect to
earlier estimates based on various more approximate calculations have
been discussed by Mandy \& Martin (1993), pointing to a need for revised
cooling and CID rates. A novel feature of the new rate coefficients is
that the set is complete; examination of the complete results has called
into question the assumption that only a few types of collisional
transitions are important (Mandy \& Martin 1993).  Therefore in our
master equation calculations no restriction has been placed on the
number of \vj\ states or on the types of collisional transitions
connecting them.

The calculations reported here, and their earlier counterparts, are
based on finding steady-state populations at a particular temperature
and density. Solution of the master equation is also an integral part of
simulations of more complex situations where the populations and
physical conditions are changing in time and/or space (e.g.,
interstellar shocks and photodissociation regions), in which case the
steady-state populations for the local temperature and density sometimes
might not be achieved (Lepp \& Shull 1983, Chang \& Martin 1991).  In
any case, incorporation of the new rate coefficients directly into such
calculations should bring quantitative improvement and the steady-state
solutions at particular temperatures and densities will provide useful
insight into more complex situations.

In \S~2 we discuss the details of our master equation calculation, the
data used in its implementation, and the numerical method used to solve
this stiff set of equations.  From fractional populations on a grid of
varying density and temperature, rate coefficients for ortho-para
interconversion, dissociation, and cooling are calculated, with a focus
on how they depend on density and temperature (\S~3).  Analytic
expressions encapsulating these dependences are summarized in the
Appendix. In \S~4 we conclude with a discussion which includes
comparisons to earlier rate coefficient estimates and, where possible,
to laboratory experiments.

\section{IMPLEMENTING THE MASTER EQUATION}
\label{sec-me}

The master equation provides the means to determine the populations of
all the \vj\ states of \h2\ by considering explicitly all processes of
interest that can populate and depopulate these states.  While, in
theory, it is possible to include all processes that can affect the
population of the states, in practice only selected ones are included.
In this study we consider radiative de-excitation by quadrupole
emission, dissociative tunneling, and the collisional processes of
excitation, de-excitation, and dissociation.  This allows a
straightforward comparison with the results of earlier studies (Lepp \&
Shull 1983, Roberge \& Dalgarno 1982, Dalgarno \& Roberge 1979,
Hollenbach \& McKee 1979) to assess the effects of the new rate
coefficients and inclusion of all transitions among the \vj\ states.

For the collisional processes, we consider only H atoms as collision
partners, corresponding to dilute \h2\ in a bath of H atoms.
Interstellar formation of \h2\ occurs on grain surfaces, with rate
coefficient $\sim 10^{-17} $ \cms\ over a range of gas temperatures
(Hollenbach \& McKee 1979; Sternberg \& Dalgarno 1989, 1995); simple two
and three body recombination is much slower.  For the interstellar warm
gas conditions in which we are interested \h2\ would not be {\it
infinitely} dilute, and we find that the fractional change in the amount
of \h2\ caused by recombination would be small on the timescale in which
steady-state fractional populations of the \vj\ states are established.
We have ignored recombination.  Nevertheless, as discussed in \S~3, it
is possible that inclusion of recombination might alter both the
relative populations of highly excited \vj\ states and any derived bulk
quantities that are sensitive to these populations.

Since interstellar \h2\ would normally be optically thin in quadrupole
lines, we have excluded reabsorption of the quadrupole emission.  We
have also excluded ultraviolet fluorescent pumping, a fundamental
process in photodissociation regions (Tielens \& Hollenbach 1985,
Sternberg \& Dalgarno 1989).  The influence of recombination and
fluorescent pumping on the rotational excitation seen via ultraviolet
absorption measurements in {\it cold} interstellar gas is an interesting
problem (e.g., Wagenblast \& Hartquist 1989) beyond the scope of this
paper.

The present calculations have relevance in assessing the conditions in
clouds in which \h2\ is not the predominant form of neutral hydrogen,
including some chemistry in the early universe (e.g., Shapiro 1992).
There are situations in which it would be desirable to consider other
collision partners for \h2, especially another \h2. This would permit an
assessment of the chemistry of molecular clouds in which \h2\ is the
dominant species.  Chemical reaction networks are sensitive to the
fractional abundances of H, \h2, and H$^+$ and hence to the extent of
dissociation of \h2.  In the case of internal energy transfer being
studied here, the behavior with collision partner \h2\ is expected to
differ significantly from that with H.  A collision with atomic H with
its open electronic shell is more likely to promote large changes in
\vj\ in the \h2\ target molecule than is a collision with the closed
shell \h2\ molecule.  Also the barrier to exchange is much lower.  Thus
it is anticipated that even when \h2\ is not infinitely dilute in H
atoms and H$_2$~+~H$_2$ collisions have to be considered, \hh2\
collisions might still dominate important aspects of the overall
kinetics.

%F Figure 1
In our model all 348 \vj\ states of \h2\ in its electronic ground state
are included.  As displayed in Figure~1, these consist of manifolds of
$j$ states for $v = 0$ to 14.  Within each manifold, quantum numbers $j$
can range from 0 to \jmax, with the actual value of \jmax\ depending on
the value of $v$ (\jmax = 38 for $v = 0$, \jmax = 4 for $v = 14$).
These states include 47 that are ``quasibound'' (Blais \& Truhlar 1979).
A quasibound diatomic molecule has an internal energy in excess of what
is required for dissociation; see, e.g., (0,32) in Figure~1.  It is
confined for a finite lifetime in this state by the rotational barrier,
through which tunneling can occur, leading to spontaneous dissociation.
When this lifetime is less than the collision time, this phenomenon is
sometimes referred to as rotational predissociation.  Dissociative
tunneling from the quasibound states can make a significant contribution
to the total dissociation rate, especially in low density regimes where
the collisions times are large.  The tunneling probability (\S~2.2)
varies strongly with the \vj\ state and can be quite large, further
increasing the stiffness of the master equation and posing special
challenges in numerical modeling.

Let $n_k$ be the population of a \vj\ state with energy $E_k$.  The sum
of these populations is \nhh\ and the fractional population $p_k =
n_k$/\nhh.  Let $n_{349}$ represent the population that has become
unbound (because \h2\ is assumed to be infinitely dilute in H,
$2n_{349}$ does not modify the H atom density \nh).  For the bound
states, the master equation for our restricted problem then has the form
%E1
\begin{equation}
\label{eq:master}
{dn_k \over dt} \ = \
\sum_{l:\ E_l > E_k} { A_{kl}n_l} \ - \ \sum_{l:\ E_k > E_l} { A_{lk}n_k} 
\ + \ n_{\rm H}\sum_{l \ne k} ( \gamma_{kl}n_l - \gamma_{lk}n_k )  
\ - \ r_{349,k}n_k \;,
\end{equation}
%
where we adopt the matrix-motivated convention that subscript $kl$
represents a transition {\it to} state $k$ {\it from} state $l$ (i.e.,
$kl \equiv k \leftarrow l$).  $A_{kl}$ is a radiative transition
probability (\S~2.3), $\gamma_{kl}$ is a collisional rate coefficient
(\S~2.1) and $r_{349,k}$ is a dissociative tunneling transition
probability (zero except for the quasibound states; \S~2.2). The
availability of electronic versions of these data is described in
Appendix~C.

Equation~1 supplemented by
%E2
\begin{equation}
\label{eq:diss}
{dn_{349} \over dt} \ = \
\sum_{l=1,348} ( n_{\rm H} \gamma_{349,l} \ + \ r_{349,l} ) \ n_l  \;.
\end{equation}
%
This coupled set of 349 equations can be expressed readily in the matrix
form required by the numerical method for solution.  We wish to
integrate the master equation to a steady state, that is until the rates
of change of the fractional populations of the \vj\ states are (close
to) zero:
%E3
\begin{equation}
\label{eq:steady}
dp_k/dt \ \simeq \ 0  \; .
\end{equation}
%
Note that as temperature and density increase, $d$\nhh$/dt$ (given by
the negative of equation 2) increasingly departs from zero; in this
study this leads to a significant depletion of the molecules at the
highest temperature and density conditions considered (\S~2.4).

\subsection{Collisional Rate Coefficients}

The only collisional processes considered in this study are those of
molecular \h2\ interacting with H atoms.  The interaction potential
energy surface for \hh2\ is well known for the ground electronic state
(Liu 1973, Siegbahn \& Liu 1978, Truhlar \& Horowitz 1978; Varandas
\etal\ 1987; Boothroyd \etal\ 1991, 1995).  The rate coefficients
$\gamma$ used here were determined from cross sections which had been
calculated with the quasiclassical trajectory (QCT) method using the
LSTH potential energy surface (see Mandy \& Martin 1993 and references
therein for details).  Limited comparisons of various revised surfaces
and (low energy) cross sections based on them indicate that the LSTH
calculations should be quite accurate (Mandy \& Martin 1991, 1992; Keogh
\etal\ 1992, Sun \& Dalgarno 1994; Lepp, Buch \& Dalgarno 1995;
Boothroyd \etal\ 1995).

Outcomes of \hh2\ collisions fall into three classes: non-reactive, in
which an inelastic collision leaves the original \h2\ molecule in a
different \vj\ state which must be of the same nuclear spin parity;
exchange, in which the atomic H collision partner combines with one of
the H atoms of the target diatomic to form an new \h2\ molecule in some
\vj\ state (among such collisions are those which change the parity,
which we refer to as ortho-para interconversion); and dissociation, in
which the molecular bond is severed leaving three H atoms.  All of these
processes are included here.  Rate coefficients have been evaluated at
15 temperatures between 450 and 45,000~K (about 180,000 coefficients per
temperature, 120,000 exchange plus 60,000 non-reactive which were
evaluated separately in the QCT method).

Some caution must be exercised in using the rate coefficients at the
high end of this temperature range.  Although most excited electronic
state potential energy surfaces for \hh2\ are well separated from the
ground electronic state potential energy surface considered here, there
is one that has a conical intersection with the ground state surface for
$D_{3h}$ geometries, the lowest of which is at 2.7 eV, about 2/3 of the
way to the dissociation limit (see Fig.~1).  The possibility that
collisions are not electronically adiabatic is expected to affect
collisional cross sections at high energy, and hence rate coefficients
at higher temperatures, above $\sim 10^4$~K.  While this has not been
explored in any detail, available information suggests that the effect
will be less than 50\% for state-to-state transitions (Blais, Truhlar,
\& Mead 1988; D. Truhlar 1991, private communication).  Wu \& Kuppermann
(1993) have shown that quantum interference related to the conical
intersection can have a measurable (and oscillating) effect on the cross
section at total energies as low as 1.8~eV.  Although the net effect on
non-reactive or reactive rate coefficients has not been assessed
quantitatively it is unlikely to be large at the temperatures in which
we are interested (A. Kuppermann 1995, private communication).

There are other unmodeled complications that could arise in a high
temperature gas, including electronic excitation of the H collision
partner, and increasing fractional ionization in the ambient gas which
would introduce electrons and protons as other collision partners.  It
can also be noted that at high temperatures collisional dissociation is
rapid and \h2\ is short-lived. When we carried out calculations at
45,000~K, we found \h2\ largely destroyed by the time steady-state
populations were achieved (\S~2.4).  Furthermore, because of rapid
dissociative cooling, a simulation at a constant high temperature might
not be relevant in many astrophysical applications.

Caution must also be exercised at the low end of the temperature range
because the rate coefficients become increasingly sensitive to the the
energy dependence of the cross section near threshold.  The statistical
noise in the QCT cross sections is largest when the cross sections are
small (as they usually are near the threshold for a transition).  Also
the QCT approach, by its very nature, does not include quantum effects
which are expected to become increasingly important at lower T.  To
gauge the magnitude of this effect, where available the quantum cross
sections for reactive transitions between (0,0) and a few of the lower
energy \vj\ states were substituted in the calculation of the rate
coefficients (Mandy \& Martin 1993).  At 1,000~K the changes in the
values of the rate coefficients caused by the substitution of the
quantum results amounted to 25\% rms and 53\% maximum; at 10,000~K these
decreased to 8\% and 14\% respectively.  As mentioned, master equation
solutions were obtained for $T$ as low as 450~K.

For physically meaningful results from the master equation, the
collisional rate coefficients used must obey detailed balance
(microscopic reversibility; Mandy \& Martin 1993) to numerical accuracy.
Calculating the upward coefficient $\gamma_{lk}$ ($E_l > E_k$) to
machine accuracy from the downward coefficient $\gamma_{kl}$ using
detailed balance ensures that this is the case.

\subsection{Dissociative Tunneling}

The \vj\ states which are quasibound (or rotationally predissociated)
have a finite probability of tunneling through the rotational barrier to
form an unbound dissociated state. The dissociative tunneling transition
probabilities have been found from
%E4
\begin{equation}
\label{eq:tunn}
r_{349,k} = 2 \pi c \Gamma_k \;,
\end{equation}
%
where $\Gamma_k$ is the energy width of state $k$ (measured here in
cm$^{-1}$) as calculated by LeRoy \& Bernstein (1971) and Schwartz \&
LeRoy (1987).

It is possible to classify the quasibound \vj\ states into two broad
groups based on the magnitude of the dissociative tunneling transition
probability relative to that for radiative de-excitation or collision.
For most quasibound states, dissociative tunneling is the most probable
of the processes unless the density is very high (we have included the
collisional [$\gamma_{kl}$] coefficients in the calculation for all
quasibound states, although their role is significant only at higher
densities).  In the interstellar medium and in some laboratory
conditions, a molecule excited to one of these \vj\ states is
effectively dissociated.

The second group is composed of the few states with high $j$ and
internal energies only slightly above the dissociation energy. These
states have much smaller tunneling transition probabilities that are
comparable to or smaller than the Einstein A-values and much smaller
than the probability of collisional energy transfer at the densities of
interest, and so to avoid an artificial buildup of their populations it
is necessary to treat them in the same manner as the classically bound
states by including their radiative transition probabilities (see
\S~2.3). The six \vj\ states so treated were (0,32), (0,33), (1,31),
(2,29), (4,26), and (5,24).

\subsection{Radiative Transitions}

Radiative emission from a \vj\ state of \h2\ obeys the quadrupole
selection rule: $\Delta j = 0, \pm 2$, except in the case of $j$ = 0
where only $\Delta j$ = 2 is allowed.  Thus radiative transitions cannot
lead to ortho-para interconversion.  There is no restriction on the
change in vibrational quantum number save through the requirement that
the energy of the final state be below that of the emitting state.  We
used the $A_{kl}$ of Turner, Kirby-Docken, \& Dalgarno (1977),
supplemented for some higher states by further calculations (A. Dalgarno
1982, private communication).  Values were required for several \vj\
states with internal energy close to the dissociation limit, namely
(0,31), (3,27) and the six bound-like quasibound states mentioned above
(\S~2.2).  Their $A_{kl}$ values were extrapolated from adjacent
smoothly-behaved values.

\subsection{Numerical Solution}

This system of first-order differential equations is stiff in that the
transition probabilities involved span many orders of magnitude as do
the populations of the \vj\ states.  The routine D02EJF in the NAG
Fortran Library (Fox \& Wilkinson 1990) provides an appropriate method
to solve this system of equations.  This routine uses a variable-order,
variable-step method implementing the backward differentiation formulae.
The starting conditions are unimportant to the steady-state result.  In
practice, we initialized the populations to the corresponding Boltzmann
populations if both $n_H > 100$~\cm3\ and $T > 10000$~K; otherwise we
populated only the lowest ortho and para states at 0.75 and 0.25,
respectively.  The integration proceeds until a user-specified function
of the solution is zero. The function we used takes the maximum of the
fractional changes in the non-quasibound populations after each time
step and divides it by the time step (see equation 3).  A steady state
is deemed to have been reached when this maximum is below some specified
value, and the integration is terminated. Sensitivity tests showed that
by using this criterion, and for the integration tolerances chosen, the
steady-state populations of the states could normally be determined to
better than four significant figures.

However, at higher temperatures and densities it was necessary to adopt
another criterion to stop the integration.  Under these physical
conditions the probability of dissociation of \h2\ was sufficiently high
that our solution without recombination could approach $p_{349}=1$
without the above precise termination criterion being met; negative
values for the populations of some of the 348 states could then occur as
numerical artifacts.  Instead, the integration was stopped when 99\% of
the initial \h2\ had dissociated.  For the range (grid) of physical
conditions studied (\S~2.4.1), this affected only the 45,000~K model at
$n_H = 10$~\cm3, those for $T \ge 20,000$~K at $n_H = 10^{3}$~\cm3\ and
those for $T \ge 10,000$~K at $n_H \ge 10^{6}$~\cm3.  In these cases the
solutions were examined to assess how closely they appeared to be
approximating steady-state conditions and it was found that the
fractional populations were converged to better than 0.5\%.

Since the steady-state populations of the excited states for all models
were determined to better than 0.5\%, the cooling, dissociation and
ortho-para interconversion rates which were calculated directly from
these populations have similar (or better) accuracy.  Thus, our solution
of the master equation did not add significantly to the uncertainties
already implicit in the calculation due to the uncertainties in the
collisional rate coefficients, the dissociative tunneling probabilities,
or the radiative emission probabilities.

The large values of some of the dissociative tunneling probabilities
contributed significantly to the stiffness of the system.  To facilitate
the numerical solution it was possible to cap their values at a factor
10$^6$ times larger than the greater of the A-value or the collisional
probability (at that density) and still not affect significantly the
solution.  A consequence is that the quasibound states with large (true)
values of the probability for dissociative tunneling are overpopulated
relative to the true steady-state solutions.  The dissociative tunneling
rate is the product of the smaller adopted coefficients and the
artificially larger populations of the quasibound states. Thus we were
able to reproduce the true rate in our calculations; in each of the test
cases we considered agreement was significantly better than the 0.5\%
from our worst-case convergence.  Radiative cooling or CID from these
overpopulated states was insufficient to affect populations of other
states; since the affected quasibound states are still underpopulated
relative to the other excited states, there is no detrimental effect on
the assessment of cooling or dissociation from the whole system
(contributions to these global quantities are again much less than
0.5\%).

\subsubsection{Range of Physical Conditions}

The range of physical conditions that can be considered is constrained
by the validity of the collisional rate coefficients (450~K -- 45000~K;
see \S~2.1) and, practically, by whether the solution of the master
equation converges to a steady state.  Both of these factors serve to
restrict the temperature range, but are compatible.  We found that below
450~K a steady state could not be attained within a reasonable amount of
computer time. Above 45000~K solutions indicate that the molecular
dissociation is virtually complete before the fractional populations of
the \vj\ states come to within 0.5\% of the steady state.

The range of densities was chosen to be \nh~=~10$^{-3}$ to 10$^9$~\cm3.
We found that this range adequately spanned the regime where we wished
to assess the competing effects of collisional processes and radiative
decay.  Outside this range, steady-state populations did not change
significantly with density (unless the density becomes extremely high,
$\approx 10^{20}$~\cm3\ and greater; there the collisional rates become
comparable to the rate of dissociative tunneling and there is again a
strong dependence of populations on the density).

\section{RESULTS} 

Steady-state populations were found for a grid of physical conditions:
15 temperatures from 450~K to 45,000~K (450, 600, 780, 1000, 1400, 2000,
3000, 4500, etc.), and 13 densities spanning the range \nh~=~10$^{-3}$
to 10$^9$~\cm3 at order of magnitude intervals.  A selection of the
results from the total of 195 simulations will be presented, making use
of graphical displays to support our qualitative interpretation and
discussion.  The fundamental data are the population distributions which
are described in the first subsection.  We are also especially
interested in such bulk quantities as the rates of dissociation,
cooling, and ortho-para interconversion.  These are calculated once the
populations are determined. Observing how the population distribution
changes with density and temperature is essential to understanding the
behavior of the bulk quantities.  To provide quantitative results, these
rates have been parameterized using interpolating functions of
temperature and density (see Appendix).

\subsection{The Populations}

%F Figure 2

Figure~2 shows the energy levels of Figure~1 contoured as functions of
$v$ and $j$. This basic behavior serves as background and orientation
for the displays of \vj\ population distributions which will be
contoured in similar diagrams.  States to the lower left are more
tightly bound and so will have larger populations.  Because of the
roughly three to one ratio of ortho to para hydrogen (\S~3.2), it is
necessary to plot their distributions separately, but other than this
scaling they are very similar in shape.

%F Figure 3a

At high density, where collisional energy transfer is the dominant
process, the populations are closest to the Boltzmann distribution.
This is illustrated in Figure 3a for the steady state at $T=2000$~K and
$n_H=10^9$~\cm3.  The contours of the distribution roughly parallel the
energy contours (Fig.~2), except at low $j$ where the states have low
statistical weight.  However, the distribution is not quite Boltzmann.
As a result of dissociative tunneling, most quasibound states are highly
depleted; note the normal behavior of the ortho bound-like quasibound
states (0,33), (1,31), (2,29) which have low tunneling rates.
Furthermore CID depletes, less dramatically, those bound states near the
dissociation limit, especially at high $v$ and low $j$ where the CID
rate coefficients are largest (see below).  Note that in our simulations
neither mode of depletion is compensated by recombination.  Because
these relative (fractional) populations within \h2\ sum to unity there
is a corresponding slight overpopulation of the lowest energy \vj\
states; this is most pronounced at high temperature.  A consequence of
the underpopulation of the highest states is that the average internal
vibrational-rotational energy of the \h2\ is also less than the thermal
equilibrium value $kT$ for each mode.  The rotational mode can be
substantially excited toward this thermal equilibrium value by 1000~K
but the vibrational mode is never fully excited, even at the higher
values of $T$ (see also Chang \& Martin 1991, Fig.~6).

At low density, where the probability of collision is sufficiently low
that radiative emission is the dominant process, any collisionally
excited \h2\ molecule is likely to lose all internal energy (in excess
of the zero point energy) by a radiative cascade before undergoing
another collision.  Thus the populations are far from Boltzmann, mostly
occupying the two ground \vj\ states (0,0) and (0,1). In the slightly
excited \vj\ states there are {\it trace populations which scale
linearly with density}; the range of these excited states becomes more
extensive at higher temperatures.

%F Figure 3b

At densities intermediate between these extremes, there is a competition
between radiative deexcitation and collisional redistribution.
Figure~3b shows the population distribution for $T=2000$~K and
$n_H=10^3$~\cm3.  It has the same characteristic shape as Figure~3a, but
the excited-state populations are very much depressed (hence the
contours are more closely spaced). For example the contours $10^{-4}$
and $10^{-9}$ of Figure~3a more or less overlie the contours $10^{-6}$
and $10^{-12}$ in Figure~3b.  Comparison of the horizontal spacing of
the contours at fixed $v$ in Figure~3b relative to Figure~3a shows that
the ratios of populations within a given $v$ manifold are not generally
equal to the Boltzmann ratios (Lepp \& Shull 1983); particularly for the
higher $j$ states, which are at increasingly higher energies (Fig.~1),
the populations are relatively depressed.

%F Figure 4

Figure~4 illustrates how there is a (temperature-dependent) ``critical
density'' $n_{{\rm H},k}$ for each particular \vj\ state $k$, at which
the probability of depopulation by radiative emission (Fig.~4a; no $T$
dependence) is equalled by that from collisional energy transfer
(combined rate coefficients contoured in Fig.~4b for 2000~K):
%E5
\begin{equation}
\label{eq:ncrit}
n_{{\rm H},k}\sum_{l=1,349} \gamma_{lk} = 
\sum_{l:\ E_k > E_l} A_{lk} \; ,
\end{equation}
%
(or more compactly $n_{{\rm H},k} \Gamma_k = A_k$).  Thus the ratio of
the values in Figures~4a ($A_k$) and 4b ($\Gamma_k$) gives the critical
density shown in Figure~4c (we ignore dissociative tunneling for the
quasibound states).  In the derived (bulk) quantities to be discussed
below, behavior related to critical density is exhibited as a region of
sharp change in the dependence on density. Since individual states have
different critical densities, and since the population distributions (as
well as the critical densities) depend on temperature, the position and
sharpness of this transition from low to high density asymptotic
behavior exhibits temperature dependence.

The uncertain contribution of recombination to the excitation of
\h2\ is often excluded (e.g., Burton, Hollenbach, \& Tielens 1990),
as we have done here.  The process of \h2\ formation on grain surfaces
is subject to considerable uncertainty (Sternberg \& Dalgarno 1989,
1995), with the rate coefficient estimated to be $R \sim 10^{-17}$ \cms\ for
the conditions of interest here.  While the newly recombined \h2\ will
be internally excited, various conflicting suggestions have been made as
to the details of this excitation distribution.  For example, Sternberg
\& Dalgarno (1989) adopt the simple assumptions made by Black \&
Dalgarno (1976) that there is equipartition with one third of the
binding energy available for \h2\ internal excitation and that this
energy is statistically distributed. With these assumptions, the
branching coefficients $\beta_k$ into state $k$ follow a Boltzmann
distribution at about 17000~K; for (ortho) states within about 1~eV of
the dissociation limit $\beta_k \sim 3 \times 10^{-3}$ (relatively
higher for high $j$ at the same energy) leading to entry rate
coefficients $R \beta_k \sim 3 \times 10^{-20}$~\cms. On the other hand,
Wagenblast \& Hartquist (1989) adopt a model advanced by Duley \&
Williams (1986) in which there is very little rotational excitation ($j$
is 0 or 1) and $v$ is intermediate ($\sim 7$); these states are about
1.5~eV or more below the dissociation limit (Fig.~1) and might have
$\beta_k$ as large as $10^{-1}$.  Clearly there is a lot of parameter
space to explore to see what effects might arise. This is beyond the
scope of this paper, but we can make an estimate to illustrate the
potential importance.

The rate coefficient due to recombination, $R \beta_k$, which affects
$dn_k/dt$ can be compared directly to the exit rate coefficient from
state $k$ due to collisions and radiative emission, $(n_{\rm H_2}/n_{\rm
H})p_k (\Gamma_k + A_k/n_{\rm H})$, which for steady-state values of
$p_k$ balances the entry rate coefficient.  Consider states about 1~eV
below the dissociation limit, for which $A_k \sim 3 \times
10^{-6}$~s$^{-1}$ (Fig.~4a) and $\Gamma_k \sim 10^{-9}$~\cms\ at 2000~K
(Fig.~4b).  Above the critical density ($\sim 3 \times 10^3$~\cm3;
Fig.~4c), $p_k \sim 10^{-9}$ (Fig.~3a). In the case used by Sternberg \&
Dalgarno, we then compare $R \beta_k \sim 3 \times 10^{-20}$~\cms\ to
$10^{-18} (n_{\rm H_2}/n_{\rm H})$~\cms, indicating that recombination
contributions would be competitive if $(n_{\rm H_2}/n_{\rm H}) < 3
\times 10^{-2}$.  At lower densities like $10^3$~\cm3, $A_k$ comes into
play, $p_k$ becomes smaller ($\sim 10^{-12}$ in Fig.~3b) and the limit
on $(n_{\rm H_2}/n_{\rm H})$ rises somewhat (note that $A_k p_k/n_{\rm
H_2}$ stabilizes at low densities).

Qualititatively, the states likely to be most affected are those with
the lowest populations (high energy states, low temperature and/or low
density conditions); their total population would be negligible compared
to the population of lower energy states.  Thus we expect that only
those quantities that depend sensitively on the populations of these
highly excited states might be affected by recombination. These would
include line emission from these states (not addressed here) and
dissociation (\S~3.3 and 4.2), but not ortho-para interconversion
(\S~3.2) or radiative cooling (\S~3.4).

There are of course other complications that could be considered too.
When recombination is in a steady state with ultraviolet
photodestruction of \h2, the accompanying fluorescent excitation of
molecules escaping photodestruction can also affect the populations of
the high-lying energy states. In a plasma, production of excited \h2\ by
associative detachment of H$ + $H$^{-}$ might become interesting too
(Black, Porter, \& Dalgarno 1981). Thus model assumptions should be
examined critically before the results are widely applied: {\it caveat
emptor}.

\subsection{Ortho-para Interconversion}

Nuclear spin considerations determine whether a molecule of \h2\ is in
an ortho or para state.  Para \h2\ has antiparallel nuclear spins which
restricts the rotational quantum number $j$ to even values.  Parallel
nuclear spins are characteristic of ortho \h2\ with odd values of $j$.
These states have statistical weight $g = (2s+1)(2j+1)$, where the
nuclear spin factor $s$ is 0 for para and 1 ortho.  In our modeling of
\h2\ infinitely dilute in H, ortho-para interconversion occurs only as
the result of the exchange reaction.  Quadrupole emission and
non-reactive energy transfer do not change the nuclear spin and thus can
change $j$ by only even values.

The rate coefficient for exchange from ortho to para hydrogen,
$\gamma_{po}$, is simply obtained:
%E6
\begin{equation}
\label{eq:gpo}
\gamma_{po} = \sum_{k:\ ortho} \ \sum_{l:\ para}  \gamma_{lk}p_k \, . 
\end{equation}
%
A similar equation holds for para to ortho conversion and for
steady-state populations $\gamma_{po} = \gamma_{op}$. Independent
evaluation of the two rate coefficients (including for Boltzmann
populations) provides useful numerical checks.  In Figure~5 we show loci
of $\gamma_{po}$ as a function of density $n_{\rm H}$ for fixed values
of temperature $T$.  Changes in the bulk quantity $\gamma_{po}$ reflect
changes in the steady-state populations for the different physical
conditions. The timescale for conversion is $1/(\gamma_{po}n_{\rm H})$.
We found that these population-weighted values were dominated by the
contributions from the lower \vj\ states; furthermore, there is only a
small change as the populations change over a wide range. Because
recombination affects primarily the populations of the upper states, we
expect the exclusion of this process from our model not to affect this
bulk rate coefficient for ortho-para conversion.

%F Figure 5

In the low density limit, any collisional excitation (or exchange) is
followed by radiative emission which does not contribute (further) to
the exchange reaction. The conversion rate coefficients $\gamma_{op}$
and $\gamma_{po}$ are simply the sums of the state-to-state exchange
rates from (0,0) and (0,1), respectively (see equation~7).  Classically
a total energy of 0.424~eV (0.155~eV relative to the zero point) is the
minimum required for exchange, although quantum mechanically tunneling
though this barrier can occur. This energetic barrier is reflected in
the steep temperature dependence at lower temperatures.  The temperature
dependence weakens at higher temperatures as transitions to higher
energy states become accessible.

At higher densities excited states can be significantly populated
(depending on $T$) and so come into play in the conversion rate
coefficient. The rate coefficients from the excited states tend to be
higher with increasing energy of the state, accounting for the increase
in $\gamma_{po}$ with density (at a given $T$).  This is most pronounced
for higher temperatures where the excited-state populations become most
significant.  However, the higher density $\gamma_{po}$'s are not
monotonic with temperature above 14000~K because of the increased
importance of dissociation from the excited states near the dissociation
limit which depletes their populations.  For the same reason, these high
density ($n_{\rm H} = 10^9$~\cm3) ``limiting'' rate coefficients are
depressed relative to those obtained assuming Boltzmann populations.

\subsubsection{The density dependence of the ortho to para ratio}

%F Figure 6

In Figure~6 we consider the ratio of the population of ortho states to
that of para states.  If thermal equilibrium were reached (at least for
the states most relevant to the conversion), then at high temperatures
(many excited states accessible) a ratio near 3 would be expected on the
basis of nuclear spin considerations.  For Boltzmann populations at or
above 450~K, this is indeed the case and this ratio is approached in our
high density limiting cases too.  At lower temperatures, the ratio
decreases because of the diminishing Boltzmann factor for state (0,1)
with respect to (0,0) in the partition function (the energy difference
corresponds to 170.5~K).

As the density decreases the populations deviate further from the
Boltzmann distribution and in the limit only states (0,0) and (0,1) are
populated. This limit is reached by $n_{\rm H} = 10^{-3}$~\cm3\ except
for the highest temperatures in our grid; at 45000~K $p_{(0,2)} = 1
\times 10^{-2}$ and $p_{(0,3)} = 2 \times 10^{-3}$.  Without going
through the full master equation numerical solution (and hence as a
check), the fractional populations of the two ground states in this
limit are determined by solving
%E7
\begin{equation}
\label{eq:gpol}
\sum_{k:\ ortho} \gamma_{k,(0,0)}p_{(0,0)} \, = 
\sum_{k:\ para} \gamma_{k,(0,1)}p_{(0,1)} 
{\hskip 25 pt} {\rm and} {\hskip 25 pt}
p_{(0,0)} + p_{(0,1)} = 1 \, .
\end{equation}
%
The left-hand equation evaluates the steady-state condition
$\gamma_{op}=\gamma_{po}$ in this low density limit.  The low-density
ortho to para ratio is now simply $p_{(0,1)}/p_{(0,0)}$.

In the {\it additional} limit of very low temperature, it is often
assumed that only the transitions between $k \equiv (0,1)$ and (0,0),
respectively, are significant in equation~7.  Then by detailed balance
for the state-to-state $\gamma_{(0,1),(0,0)}$ and its inverse
counterpart the ratio $p_{(0,1)}/p_{(0,0)}$ (whence the ortho to para
ratio) would be $9 \exp(-170.5/T)$ (Mandy \& Martin 1993).  Note that at
low temperature this low density approximation would give precisely the
same value as the thermal equilibrium ratio (the factor 9 is misprinted
as 3 in equation~2 of Flower \& Watt 1984.)\ At 450~K this approximation
would give a ratio 6.2.  But the energy of state (0,2) relative to (0,1)
is equivalent to 339.5~K and so this low temperature limit is not
reached on our grid (450~K and above). Because of the accessibility of
the (0,2) state (and others) the ratio is lower than the simple
low-temperature estimate, though still typically larger than 3. With
increasing temperature the low-density ortho to para ratio becomes
closer to 3 (and even slightly lower) as more states (in the above sums)
become accessible to exchange transitions from the ground states.

\subsection{Dissociation}

In our master equation solution, \h2\ dissociates as the result of
collisions to the continuum of classically unbound states (CID), or
through dissociative tunneling from quasibound states.  The {\it total}
dissociation rate coefficient for the dissociation of
\h2\, $\gamma_d$, is:
%E8
\begin{equation}
\label{eq:gdiss}
\gamma_{d} = \gamma_{cid} + \gamma_{dt}
= \sum_{k=1,348}  \gamma_{349,k}p_k 
+ \sum_{k:\ quasibound}r_{349,k}p_k/n_H \, .
\end{equation}
%
%F Figures 7 and 8
%
Figures~7 and 8 shows $\gamma_{cid}$ and $\gamma_{dt}$, respectively, as
functions of \nh\ for various $T$. The timescale for dissociation is
$1/(\gamma_d n_{\rm H})$.  

In the low density limit, $\gamma_{cid}$ corresponds to dissociation
directly from (0,0) and (0,1), the only populated states.  Therefore, an
independent way of calculating $\gamma_{cid}$ without integrating the
master equation (and to check it) is
%E9
\begin{equation}
\label{eq:lowcid}
\gamma_{cid}  = \gamma_{349,(0,0)}p_{(0,0)} + \gamma_{349,(0,1)}p_{(0,1)}  
\; ,
\end{equation}
%
where $p_{(0,0)}$ and $p_{(0,1)}$ are from equation~7.  Similarly since
dissociative tunneling arises from quasibound states excited only from
the ground states
%E10
\begin{equation}
\label{eq:lowdt}
\gamma_{dt}  =  \sum_{k:\ quasibound}(1 - b_k)\
[\gamma_{k,(0,0)}p_{(0,0)} + \gamma_{k,(0,1)}p_{(0,1)}]  
\; .
\end{equation}
where $b_k = \sum_{l:\ E_k > E_l} A_{lk} / ( r_{349,k}+\sum_{l:\ E_k >
E_l} A_{lk})$ is the fraction of quasibound molecules in state $k$ that
are stabilized by the quadrupole emission.  Except for the six states
mentioned in \S 2.2, dissociative tunneling is much more probable than
quadrupole emission (branching factor $b_k = 0$), but as a refinement to
the treatment in Dove \& Mandy (1986a) it is important to allow for this
branching to match the master equation results precisely.  Equation~10
shows directly that in this low density limit $\gamma_{dt}$ is
independent of \nh\ (see Fig.~8); this can be reconciled with the
equation~8 by noting that the quasibound populations $p_k$ scale as \nh\
in this limit.  At low temperatures, since slightly less energy is
required to reach a quasibound state (which then undergoes dissociative
tunneling) than for direct CID, $\gamma_{dt} > \gamma_{cid}$. This
inequality is reversed at grid temperature 4500~K and above.

With increasing density, a significant part of the dissociation takes
place via ``ladder climbing'', a sequence of inelastic collisions in
which a molecule is first excited, with interconversion of vibrational
and rotational energy, and then dissociated by CID or transfer to a
quasibound state followed by dissociative tunneling.  In the master
equation solution, this is represented by the increased, and now
significant, steady-state populations $p_k$ of the excited states from
which dissociation by either route can take place (equation~8).
%
%F Figure 9a
%
Figure~9a shows the coefficients $\gamma_{349,k}$ state by state for the
intermediate temperature 2000~K. Note the huge increase relative to the
rate coefficient from the ground state; because of this $\gamma_{cid}$
increases with density for fixed $T$.  Internal energy is more effective
in promoting dissociation than is the same amount of translational
energy (Dove \& Mandy 1986b), with a slight preference for the
vibrational form.  The latter can be seen by comparing the shapes of the
contours in Figure~9a with those for constant internal energy in
Figure~2: for a given internal energy, $\gamma_{349,k}$ is higher for
higher $v$ and lower $j$.
%
%F Figure 9b
%
The populations $p_k$ for $T=2000$~K and \nh\ = $10^{3}$~\cm3\ have been
shown in Figure~3b; they decrease rapidly toward the continuum, opposite
to the trend for $\gamma_{349,k}$ seen in Figure~9a. The
population-weighted contribution to $\gamma_{cid}$, $\gamma_{349,k}p_k$,
is contoured in Figure~9b.  For low density, all of the CID dissociation
is from the ground state(s).  The value $10^{3}$~\cm3\ is a transitional
density where the highly excited {\it classically-bound} states near the
dissociation limit begin to dominate; this dominance is more pronounced
with increasing density as the excited-state populations increase
further. The increase of $\gamma_{cid}$ with density can be seen in
Figure~7, with the characteristic ``critical density'' behavior.  The
increase levels off at $n_{\rm H} \sim 10^6$~\cm3\ as the
classically-bound populations become collisionally dominated and
stabilize.

Figure~8 shows that $\gamma_{dt}$ also increases with density through a
critical region and then also reaches a plateau.  This plateau can be
reconciled with equation~8 by recalling that the classically-bound
states from which the quasibound states are being populated have been
stabilized; thus the quasibound populations $p_k$ scale as \nh\ just as
in the low density limit.  Except in the limit of very high density
where the quasibound states can be thermalized (well beyond our grid
which ends at \nh$=10^9$~\cm3; see \S~4.2.2), dissociative tunneling
{\it from these states} is much more probable than CID.

Also included (for the quasibound states) in the quantity contoured in
Figure~9b is $r_{349,k}p_k/n_H$, the contribution to $\gamma_{dt}$. CID
is negligible from the quasibound states allowing an unambiguous display
of the dissociative tunneling which is seen to be very relevant; thus
$\gamma_{dt} > \gamma_{cid}$ for $T=2000$~K and \nh\ = $10^{3}$~\cm3\ as
already seen for all \nh\ from Figures~7 and 8.  Individual quasibound
states tend to make comparable contributions to $\gamma_{dt}$ despite
vastly different populations, since those quasibound states with high
tunneling probabilities have low steady-state populations and
conversely. It is more a matter of how quickly they are populated by
collision from the lower-lying bound states; as seen in Figure~9b, the
maximum contribution is from (0,33) and other low lying quasibound
states like (7,21) which are most easily populated at 2000~K.  With
increasing $T$, the maximum in $\gamma_{349,k}p_k$ moves to lower $v$
and higher $j$ and the relative importance of $r_{349,k}p_k/n_H$
decreases.

At high density the relative importance of $\gamma_{dt}$ and
$\gamma_{cid}$ is like that at low density, although $\gamma_{dt}$ is
less dominant at low temperature and the reversal to $\gamma_{dt} <
\gamma_{cid}$ is at grid temperature 6000~K.  Compared to the
temperature dependences at low density, those at higher density are
compressed because transitions involving excited states are important.

Our model has excluded the effects of recombination (and fluorescence)
which might significantly affect the populations of the upper states
(\S~3.1).  The contribution of the upper states to the value of
$\gamma_d$ can be significant as we have shown, and so our values should
be used in other scenarios with caution. It can also be noted that the
cases most likely to be affected are those with low upper-state
populations, and these coincide with low or negligible $\gamma_d$.

\subsection{Cooling Rates}

The translational kinetic temperature $T$ is usually a distinct
parameter of importance in characterizing the physical conditions, and
so in this study cooling is considered to be a lowering of $T$ or the
removal of translational energy, $E_{tr}$.  In general, the lowering of
$E_{tr}$ is balanced exactly by changes in the internal energy $E_{ie}$
of the H$_2$ molecules (the {\it net} effect of collisional
redistribution and dissociation, and of radiative deexcitation) and
changes in the radiation field $E_r$ through radiation by excited \h2\
molecules.  Since we are considering steady state populations the
translational cooling arises from the usual irreversible effects of
radiative transitions and collisions leading to dissociation.  The
cooling rate per unit volume (\ecms) is usually called $\Lambda$, which
we separate into components related to internal energy ($\Lambda_{ie}$)
and radiation ($\Lambda_r$):
%E11
\begin{equation}
\label{eq:ecool}
\Lambda = -dE_{tr}/dt = dE_{ie}/dt + dE_{r}/dt = \Lambda_{ie} +
\Lambda_r \, .
\end{equation}
%
We report the cooling rate using our steady-state populations of \h2\
which were obtained assuming that the translational temperature, $T$, is
constant, consistent with the coolant, \h2, being infinitely dilute in
the H gas. In actual applications one would want to investigate whether
the steady-state populations at that $T$ would be reached on the cooling
timescale (in this case $n_{\rm H}kT/\Lambda$ since the energy of the
trace amounts of \h2\ is negligible).

\subsubsection{Internal Energy}

For steady fractional populations it follows from equation~3 that
$\Lambda_{ie} = -Ddn_{{\rm H}_2}/dt$, where $D$ is the absolute value of
the average binding energy of the \h2\ molecules for these populations.
This result also makes intuitive sense. For steady-state relative
populations the internal energy of an average molecule is constant; only
the number of molecules changes.  Since we have not included
recombination (and the accompanying heating) in our calculation, every
dissociation of a hydrogen molecule permanently removes the required
energy from translation.  Thus {\it in a steady state} $\Lambda_{ie}$ is
identical to what is sometimes referred to as dissociative cooling.

It is useful to define a rate coefficient $\gamma_{ie} \equiv
\Lambda_{ie}/(n_{\rm H}n_{{\rm H}_2})$ (units: \ecps). In terms of
previously-defined coefficients for CID and dissociative tunneling
(equation~8)
%E12
\begin{equation}
\label{eq:lamdiss}
\gamma_{ie} = D \gamma_d \, . 
\end{equation}
%
In earlier work (e.g., Hollenbach \& McKee 1979; Lepp \& Shull 1983) a
less precise definition has often been adopted, using instead of $D$ the
value $D_o$, the dissociation energy from the ground state (4.48~eV),
and also ignoring the contribution of dissociative tunneling to
$\gamma_d$.  In steady-state conditions in which \h2\ is excited $D$ can
be appreciably less than $D_o$ (on our grid $D$ is as small as 2.6~eV at
high $T$ and \nh).  We have already discussed (\S~3.3) the importance of
dissociative tunneling relative to CID and their behavior with $T$ and
\nh. 
%
%F Figure 10
%
Figure~10 displays $\gamma_{ie}$ in the familiar form of Figures~7 and
8. At 2000~K and below $\gamma_{ie}$ is small in comparison to
$\gamma_r$ to be discussed next. Also $\gamma_{ie}$ would be affected by
neglect of recombination in the same way as $\gamma_d$, but only where
$\gamma_{ie}$ is negligible.

\subsubsection{Radiation}

Whenever an \h2\ molecule in an excited \vj\ state emits a photon, the
energy carried away by the photon is irreversibly lost from the gas
system (which is assumed to be optically thin).  As described in
equation~13, the state-specific energy loss rate, shown in Figure~11,
depends on the Einstein $A$-values (the sum of which is shown in
Fig.~4a) and the energies of the photons emitted.
%
%F Figure 11
%
The bulk radiative cooling rate then arises from the products of these
state-specific rates with the steady-state populations (e.g., Fig.~3b).
For example, using the populations for $n_{\rm H}=10^3$~\cm3\ and $T =
2000$~K, we find that despite the bias seen in Figure~11 toward somewhat
excited states most of the radiation comes from the lowest states
because of the low excitation; there is a local maximum in energy loss
rate at state (0,7) for ortho states.  As the excitation increases with
density, the peak energy loss moves up to higher states, for example to
(1,3) and (1,5) for the ortho populations at 2000~K and $n_{\rm
H}=10^6$~\cm3\ and larger.  Increased temperature of course also raises
the excitation at constant density.  Because of this population
weighting, neglect of the potential importance of recombination for
high-lying states (\S~3.1) should have little effect on $\Lambda_r$.

A natural function to compute from the populations is the cooling per
\h2\ molecule (units: erg~s$^{-1}$):
%E13
\begin{equation}
\label{eq:rcool}
\Lambda_r/n_{{\rm H}_2} = \ \sum_{k=1,348} \
\sum_{l:\ E_l<E_k} (E_k - E_l){ A_{lk}p_k} \; .
\end{equation}
%
But for direct comparison with $\gamma_{ie}$ we plot in Figure~12 the
rate coefficient (units: \ecps) $\gamma_{r} \equiv \Lambda_r/(n_{\rm
H}n_{{\rm H}_2})$.

%F Figure 12

At low densities $\gamma_{r}$ is independent of \nh\ since the $p_k$
scale linearly with \nh (\S~3.1).  In this limit virtually all
collisional excitations are followed by a radiative cascade to (0,0) or
(0,1) (for the highest temperatures considered the quasibound states are
reached and so branching to dissociative tunneling has some effect).
Therefore, in this case it is possible to calculate $\gamma_{r}$
independently:
%E14
\begin{eqnarray}
\label{eq:gcool}
\gamma_{r} &=& 
\ \displaystyle\sum_{k:\ para} \
(E_k - E_{(0,0)})b_k(\gamma_{k,(0,0)}p_{(0,0)} + \gamma_{k,(0,1)}p_{(0,1)})  
\cr\cr
%
&+& 
\ \displaystyle\sum_{k:\ ortho} \
(E_k - E_{(0,1)})b_k(\gamma_{k,(0,0)}p_{(0,0)} + \gamma_{k,(0,1)}p_{(0,1)})
\; .
\end{eqnarray}
%
Here $p_{(0,0)}$ and $p_{(0,1)}$ are from equation~7 and $b_k$ is the
state-specific branching ratio accounting for tunneling from the
quasibound states (see equation~10; $b_k = 1$ for classically-bound
states).  Equation~14 also makes explicit that $\gamma_{r}$ is
independent of \nh\ in this low density limit and that its temperature
dependence arises from that of the state-to-state excitation rate
coefficients.

Above the critical density collisional energy transfer becomes more
influential than radiative energy transfer and steady-state relative
populations of the \vj\ states are reached which no longer change with
further increases in density.  Therefore, at high density
$\Lambda_r/n_{{\rm H}_2}$ does not depend on \nh\ and $\gamma_r \propto
n_{\rm H}^{-1}$.  Given this decrease of $\gamma_r$ with increasing \nh\
and the simultaneous growth of $\gamma_{ie}$, the latter comes to
dominate at high densities.

\section{DISCUSSION}

In this section we examine the various rate coefficients derived from
the master equation results in the context of previous estimates
(including the low density limit) and laboratory experiments
(approximating but not identical to our high density limit).  A
comparative discussion of the state-to-state collisional rate
coefficients on which this study is based was presented by Mandy \&
Martin (1993).

\subsection{Ortho-para Interconversion}

In many regions of interstellar space the ortho to para ratio might be
modified by collisions between \h2\ and protons (Dalgarno, Black, \&
Weisheit 1973; Flower \& Watt 1984).  The rate coefficient for the (0,1)
to (0,0) transition by this process is estimated to be relatively
independent of temperature and of order $3 \times 10^{-10}$~\cms.  This
can be compared to some state-to-state values given in Mandy \& Martin
(1993) and the values of $\gamma_{po}$ given in Figure~5 and
equation~A1. For example, at 1000~K $\gamma_{po} \simeq 2 \times
10^{-13}$~\cms\ for a wide range of densities.  Thus at 1000~K
collisions with H will come to dominate the conversion compared to
collisions with H$^+$ when $n_{\rm H^+}/n_{\rm H} < 10^{-3}$; at 2000~K
the ratio would be $10^{-2}$.  Thus collisions with H are important to
understanding measurements of the relative populations in the $j$ states
in warm regions (Takayanagi, Sakimoto, \& Onda 1987; Tanaka et al.\
1989).

%\subsubsection{Comparison with experiment}

Schofield (1967) has summarized some experimental measurements of the
rate of conversion of para \h2\ to the ortho form at high density and in
the temperature range 300 to 450~K.  Adopting Boltzmann values for the
para population distribution we have used our state-to-state rate
coefficients to simulate this situation; recall however that we regard
450~K as the lower bound on the range of validity of our coefficients.
Our value is (coincidentally) within 10 percent of the experimental
value at 450~K and a factor of 3 low at 300~K.  Schofield also
summarizes some theoretical calculations based on an earlier potential
energy surface and extending from 300~K to 1000~K. Compared to the
analytic function adopted by Schofield our simulations are a factor 2
low at 450~K, equal at 1400~K, and higher by a factor 4 at 4500~K. Given
the uncertainties, we consider the level of agreement acceptable.

\subsection{Dissociation}

\subsubsection{Comparison with other work}

Dalgarno \& Roberge (1979) and Roberge \& Dalgarno (1982) showed clearly
that radiative depopulation of the higher energy states (radiative
stabilization) could be significant in the interstellar medium where the
density and consequently the collisional probability was sufficiently
low to be comparable to the probability of quadrupole emission; their
particular emphasis was on the accompanying depression of the rate of
collision induced dissociation (CID).  Lepp \& Shull (1983) refined the
model for a warm gas ($T>1000$~K) by including vibrationally-excited
states and dissociation from all vibrational levels.  Not surprisingly,
different sets of rate coefficients and different levels of
approximation in the solution of the master equation have led to
discordant estimates of the rates of CID (and cooling; \S~4.2).

%F Figure 13

Figure~13 displays our results for $\gamma_{d}$ as a function of $T$ for
\nh\ = $10^{-3}$~\cm3, $10^{3}$~\cm3, and $10^{9}$~\cm3 (low density,
critical density, and high density, respectively: see Fig.~7).  This
form of display is that used by Roberge \& Dalgarno (1982) and is useful
in the context of the functional fits described in Appendix~A.  Also
shown for comparison are values from the study by Lepp \& Shull (1983)
for the same densities.  As pointed out earlier by Dove et al.\ (1987),
the temperature dependence at low density is incorrect. Furthermore, the
critical density that arises from their adopted excitation rate
coefficients is too high, so that their results for \nh\ =
$10^{3}$~\cm3\ are too close to those for the low density limit.
Mac~Low \& Shull (1986) published a modification to the CID rate
coefficient for the low density limit which is in closer agreement with
our results (Fig.~13). However, the discrepancies at intermediate
density are still present when this is combined with the critical
density formulation from Lepp \& Shull (as one infers should be
done). If instead we lower their critical density by one order of
magnitude, the behavior with density become much closer to what we find.
A lower critical density is qualitatively consistent with the
differences in the treatment of collisional transitions; we included all
rotational and vibrational-rotational transitions, some involving
smaller amounts of energy than do the vibrational-only transitions
accounted for in the calculations of Lepp \& Shull.

The uppermost curve in Figure~13 is the high-density limit from the
experimentally-determined values of Breshears \& Bird (1973), these
actually being constrained only in the limited range 3500~K to 7500~K.
This limit has been adopted as a firm constraint in these earlier
studies; the appropriateness of this is discussed next.

\subsubsection{The high density limit}

In the ``high density'' limit of this study, namely $n_{\rm H} =
10^9$~\cm3, collisional processes dominate relative to radiative
depopulation and the populations of the bound states {\it approach} the
Boltzmann distribution (\S~3.1).  In this regime, the dissociation rate
coefficient $\gamma_{d}$ is enhanced with increasing temperature because
of the increase in the populations of higher energy states.  However,
care must be exercised in comparing $\gamma_{d}$ for our high density
limit with experimentally measured values of the dissociation rate
coefficient (Breshears \& Bird).  As seen in \S~3.3 dissociative
tunneling can be an important process relative to CID even at $n_{\rm H}
= 10^9$~\cm3.  As the density increases many orders of magnitude further
to the range characteristic of laboratory experiments, $\gamma_{dt}$
fades in importance relative to $\gamma_{cid}$ (see equation~8); we have
confirmed this with some simulations with \nh\ as large as
$10^{20}$~\cm3.

Over the limited temperature range available for comparison, 3500~K to
7500~K, our ``high density'' values of $\gamma_{cid}$ are consistently
lower than those of Breshears \& Bird.  For example, at 4500~K we find
$\gamma_{cid} = 0.97 \times 10^{-13}$~\cms\ compared to $2.0 \times
10^{-13}$~\cms.  (Values of $\gamma_{d}$, which include dissociative
tunneling, are coincidentally very close: $2.0 \times 10^{-13}$~\cms\ at
4500~K, for example; see also Fig.~13).  Nevertheless, the agreement
between experiment and our purely theoretical calculations based on {\it
ab initio} CID rate coefficients $\gamma_{349,k}$ is quite remarkable,
and the slight systematic discrepancies can be understood at least
qualitatively. At $n_{\rm H} = 10^9$~\cm3, our values of $\gamma_{cid}$
are low because of an underpopulation of the upper states, which have
the largest individual CID rate coefficients (e.g., Fig.~9a). In higher
density simulations, this is gradually redressed but the values for the
CID rate coefficient are still slightly low.  If on the other hand we
simply take Boltzmann populations in combination with our individual CID
rate coefficients $\gamma_{349,k}$ we tend to overestimate the
laboratory results slightly since we have now overcompensated for the
population of the upper states.

Another reason that even our highest density simulations (including some
with radiative and dissociative tunneling rates set to zero) have
slightly depressed populations for the important high energy (including
quasibound) states is that in the simulations these states are not
repopulated by recombination, a process that is important at
experimental densities.  Schwenke (1990) has included recombination in
some higher density master equation calculations based on some similar
state-to-state rate coefficients.  One of his models (his Fig.~3) has
$n_{\rm H} \simeq 2 \times 10^{16}$~\cm3\ and $n_{\rm H_2} < n_{\rm H}$
consistent with our assumptions.  For this model at 5000~K his
time-dependent results show a ``plateau region'' at which a quasi-steady
state has been established (the timescale agrees with our time-dependent
integrations), before recombination has become important to the state
populations.  Interpolating from simulations at 4500~K and 6000~K and
$n_{\rm H} = 10^{16}$~\cm3\ we have been able to reproduce the
dissociation rate coefficient deduced from his results for the plateau
region.  Thereafter, his dissociation rate coefficient increases with
time by a factor 3 to a value close to the value we derive using simply
Boltzmann populations.

From these considerations it can be seen that it is not appropriate to
force the ``high density'' limit of $\gamma_{cid}$ (at $n_{\rm H} =
10^9$~\cm3) to match the measurements of Breshears \& Bird, as has been
done by Roberge \& Dalgarno (1982) and Lepp \& Shull (1983).
Furthermore, as discussed in \S~3.3, over the whole range one also has
to include dissociative tunneling which has been omitted in previous
studies.

\subsection{Cooling}

%F Figure 14

Early assessments of the cooling rate focussed on the role of rotational
excitation (Takayanagi \& Nishimura 1960, Nishimura 1968, Elitzur \&
Watson 1978).  For warmer gas this has to be extended by inclusion of
vibrational transitions.  Figure~14 displays our results for the
radiative cooling rate per \h2\ molecule, $\Lambda_r/n_{{\rm H}_2}$, as
a function of $T$ for fixed densities and compares them to the analytic
fits found by Lepp \& Shull (1983; see also their Fig.~3) and Hollenbach
\& McKee (1979; see also their Fig.~7 which shows $\gamma_r$).  The
uppermost curves in Figure~14 are for the ``high-density'' limit for
which the populations approach their Boltzmann values.  Our results
offer some improvement at high $T$ owing to our complete treatment of
the \vj\ states. The plateau beginning at $T = 10^4$~K does not appear
for Boltzmann populations, but occurs in our master equation results
because the high energy states near the dissociation limit remain
somewhat underpopulated as discussed in \S~3.1.  Larger differences with
respect to the earlier results occur for lower densities.  At the lowest
densities, where $\gamma_r$ becomes independent of \nh\ (Fig.~12;
equation~14), the differences can be attributed to use of different
excitation rate coefficients from the ground states and detailed
inclusion of a full suite of excited states.  Our results show that at
intermediate densities, a full master equation solution is required to
determine the populations of all states.

\subsection{Conclusions}

Over the past two decades computational speeds have improved to the
extent that it should be possible to include a full treatment of the
\vj\ populations of \h2\ in models of astrophysically interesting
environments, like shocks and photodissociation regions.  Our models
illustrate this possibility for time-dependent calculations (which in
our case lead to steady-state conditions) and show that a full
simulation is a desirable goal compared to approximate treatments.  Over
the same time there have been significant advances in the computation of
the potential energy surface for \hh2\ and in the calculation of
state-to-state (and CID) cross sections, to the point that a full set of
state-to-state rate coefficients is now available for the \hh2\
interaction (Mandy \& Martin 1995; Martin \& Mandy 1995) at intermediate
temperatures.  These rate coefficients need to be extended by quantum
calculations to lower temperatures, which is not insurmountable since
there only the lowest \vj\ states are important. Among the remaining
challenges is development of a set of state-to-state rate coefficients
for the \molh2\ interaction.

\acknowledgements 
%
This work was supported by the Natural Sciences and Engineering Research
Council of Canada. We thank R.~J.~LeRoy for an electronic version of the
dissociative tunneling probabilities and the referee for insightful
comments.

\newpage

\begin{appendix}
\title{APPENDIX: FUNCTIONAL APPROXIMATIONS}

The following fits have been made and are available in electronic form
(Appendix~C).

\section{Collisional Rate Coefficients}

Figure~5 ($\gamma_{po}$), Figure~7 ($\gamma_{cid}$), Figure~8
($\gamma_{dt}$), and Figure~9 ($\gamma_{c}$) show that the respective
rate coefficients all have qualitatively similar behavior.  For fixed
temperature there is a flat portion at low density, say $\gamma_l$, and
another higher flat portion at high density, say $\gamma_h$. The
transition occurs near some critical density $n_c$. This motivates the
following form, generalized from Lepp \& Shull (1983):
%EA1
\begin{eqnarray} 
\label{eq:fitct}
\log \gamma_{h,1} &=& 
\alpha_1+\alpha_2\log T+\alpha_3(\log T)^2+\alpha_4(\log T)^3
+\alpha_5\log(1+\alpha_6/T) \cr\cr
%
\log \gamma_{h,2} &=& \alpha_7/T \cr\cr
%
\log \gamma_{l,1} &=&
\alpha_8+\alpha_9\log T+\alpha_{10}(\log T)^2
+\alpha_{11}\log(1+\alpha_{12}/T) \cr\cr
%
\log \gamma_{l,2} &=& 
\alpha_{13}/T \cr\cr
%
\log n_{c,1} &=&
\alpha_{14}+\alpha_{15}\log T+\alpha_{16}(\log T)^2+\alpha_{17}/T \cr\cr
%
\log n_{c,2} &=&
\alpha_{18}+\log n_{c,1} \cr\cr
%
p &=& 
\alpha_{19}+\alpha_{20}\exp(-T/1850)+\alpha_{21}\exp(-T/440) \cr\cr
%
\log \gamma &=&
\log \gamma_{h,1} 
- (\log \gamma_{h,1} -\log \gamma_{l,1})/[1+(n_{\rm H}/n_{c,1})^p] \; + \cr\cr
&\phantom{+}& \log \gamma_{h,2} 
- (\log \gamma_{h,2} -\log \gamma_{l,2})/[1+(n_{\rm H}/n_{c,2})^p]
\; .
\end{eqnarray}
%
At density $n_c$, $\gamma$ is the geometric mean of $\gamma_l$ and
$\gamma_h$. 

%T Table I

Table~I contains the $\alpha$'s. Units for these parameters are not
given but are obvious from the context (c.g.s.\ units; see axis labels
in Figures and footnotes to the Table).  These parameters were obtained
by non-linear least squares fits to the computed grid of values for each
rate coefficient (low and high density limits were not optimized
separately). Grid values at 450~K were given low weight except in the
case of $\gamma_{po}$, and so the fits should be used with caution below
600~K.  To indicate the quality of the fit, the {\it rms} error along
with the maximum and minimum (signed) deviation (within the valid $T$
range) are also tabulated (all in percent).

\section{Radiative Cooling Rate}
 
Figure~11 for $\log \gamma_{r}$ shows a qualitatively different behavior
with respect to $n_{\rm H}$ and $T$. At low density there is a flat
portion which turns downward at some critical density $n_c$ to a high
density linear portion with slope $-1$. This motivates use of a
hyperbola to describe the data at fixed $T$, whence:
%EB1
\begin{eqnarray} 
\label{eq:fitra}
\log \gamma_h &=& - \log n_{\rm H} +
\alpha_1+\alpha_2\log T+\alpha_3(\log T)^2+\alpha_4(\log T)^3
+\alpha_5\log(1+\alpha_6/T) + \alpha_7/T \cr\cr
%
\log \gamma_l &=& \log \gamma_h + \log n_{\rm H} - \log n_c \cr\cr
%
f &=&  \log \gamma_l/\gamma_h \; \equiv \; \log n_{\rm H}/n_c \cr\cr
%
\log n_c &=&
\alpha_8+\alpha_9\log T+\alpha_{10}(\log T)^2+\alpha_{11}(\log T)^3
+\alpha_{12}\log(1+\alpha_{13}/T) + \alpha_{14}/T \cr\cr
%
w &=&
\alpha_{15}+\alpha_{16}\log T+\alpha_{17}(\log T)^2+\alpha_{18}(\log T)^3
+\alpha_{19}/T \cr\cr
%
\log \gamma_r &=&
\log \gamma_h + 0.5[f - \sqrt(f^2+2w^2)]
\; .
\end{eqnarray}
%
The fit is summarized in Table~I, again for 600~K and above.  The low
and high density limits, $\gamma_l$ and $\gamma_h$, were not optimized
separately.

\section{Availability in Electronic Form}

The following is a description of several files and associated test
programs that may be requested by e-mail.

A file from which the entries in Table~I were made is available, along
with a test program (in FORTRAN) that reads this file and evaluates the
various coefficients as a function of temperature and density using
equations A1 and B2.

The data assembled for solution of the master equation (equations~1 and
2) were described in \S~2.  The population $n_l$ on the RHS of the
master equation is an element of a column vector of length 348~+~1. The
elements were ordered as follows.  First come the 177 para-\h2\ \vj\
states, with an inner loop on $v$ (vertically in Figure~1):
%
(0,0), (1,0), (2,0), ..., (14,0), (0,2), (1,2), ..., (0,34), (1,34),
(0,36), (0,38).
%
Then come the 171 ortho-\h2\ \vj\ states:
%
(0,1), (1,1), (2,1), ..., (14,1), (0,3), (1,3), ..., (0,35), (1,35),
(0,37).  
%
Element 349 corresponds to a dissociated molecule.  Subroutines to
translate from \vj\ to $l$ and vice versa are available.  These can be
exercised with a test program that reads a file of the energies of the
\vj\ states (as plotted in Figure~1) to generate (elements of) the
corresponding energy column vector.

With the above ordering the matrix of elements $A_{kl}$ (radiative
transition probabilities from state $l$) contains two square blocks on
the diagonal of dimensions 177 and 171 since there are no
parity-changing quadrupole radiative transitions.  Each block is sparse
because of selection rules on $\Delta j$.  Row and column 349 are
empty. A test program to produce this matrix is available.
%
The dissociative tunneling transition probabilities $r_{349,l}$ can be
regarded as row 349 of an otherwise empty matrix.  Even this row is zero
except for the quasibound states.  In the test program it is added to
the matrix $A_{kl}$, which had row 349 empty.
%
In the test program a diagonal element $A_{kk}$ corresponding to the sum
of all radiative and dissociative tunneling transitions out of each of
the 348 states $k$ is generated as the negative of the sum along that
column.

Electronic files of the matrix $\gamma_{kl}$ with this same ordering
have been made available separately (Martin \& Mandy 1995).

\end{appendix}

\newpage

\begin{table}
\begin{center}
\caption[]{Analytic Fits of Steady-State Rate Coefficients for \hh2}
\vspace{1em}
\begin{tabular}{llllll}
\hline
\multicolumn{1}{c}{}&
\multicolumn{1}{c}{$\gamma_{po}$ \tablenotemark {a}} &
\multicolumn{1}{c}{$\gamma_{cid}$ \tablenotemark {a}} &
\multicolumn{1}{c}{$\gamma_{dt}$ \tablenotemark {a}} &
\multicolumn{1}{c}{$\gamma_{ie}$ \tablenotemark {b}} &
\multicolumn{1}{c}{$\gamma_{r}$ \tablenotemark {b}} \\
\hline
  \arm &  \phaa  1.6 &  \phaa  3.9 &  \phss  8.7 &  \phss  6.9 &  \phss  7.5 \sla
  \afx &  \phaa  3.6 &  \phaa  9.3 &  \phaa 23.8 &  \phaa 31.3 &  \phaa 19.1 \sla
  \afn &  \phbb  4.0 &  \phbb  8.8 &  \phbb 20.4 &  \phbb 17.5 &  \phbb 19.8 \sla
 \phcc   1 \tablenotemark {d} &  \pha 1.045999\phc 2 &  \phb 1.784239\phc 2 &  \phb 1.427664\phc 2 &  \phb 1.846781\phc 2 &  \phb 1.058656\phc 3 \sla
 \phcc   2  &  \phb 5.180703\phc 1 &  \phb 6.842243\phc 1 &  \pha 4.270741\phc 1 &  \pha 3.793892\phc 1 &  \pha 6.412920\phc 2 \sla
 \phcc   3  &  \pha 7.374957\phc 0 &  \pha 4.320243\phc 1 &  \phb 2.027365\phc 0 &  \pha 2.996993\phc 0 &  \phb 1.330667\phc 2 \sla
 \phcc   4  &  \phb 3.123496\phd 1 &  \phb 4.633167\phc 0 &  \phb 2.582097\phd 1 &  \phb 8.286897\phd 1 &  \pha 9.285717\phc 0 \sla
 \phcc   5  &  \phb 1.621112\phc 1 &  \pha 6.970086\phc 1 &  \pha 2.136094\phc 1 &  \pha 3.155289\phc 1 &  \pha 9.160106\phc 1 \sla
 \phcc   6  &  \pha 1.029437\phc 4 &  \pha 4.087038\phc 4 &  \pha 2.753531\phc 4 &  \pha 2.861577\phc 4 &  \pha 2.680075\phc 3 \sla
 \phcc   7  &  \phb 2.606107\phc 3 &  \phb 2.370570\phc 4 &  \phb 2.146779\phc 4 &  \phb 2.141734\phc 4 &  \pha 9.500433\phc 3 \sla
 \phcc   8  &  \phb 1.733190\phc 0 &  \pha 1.288953\phc 2 &  \pha 6.034928\phc 1 &  \pha 6.265847\phc 1 &  \phb 1.253746\phc 3 \sla
 \phcc   9  &  \phb 2.417648\phc 0 &  \phb 5.391334\phc 1 &  \phb 2.743096\phc 1 &  \phb 3.282533\phc 1 &  \pha 7.792906\phc 2 \sla
  10  &  \pha 1.532729\phd 1 &  \pha 5.315517\phc 0 &  \pha 2.676150\phc 0 &  \pha 3.275816\phc 0 &  \phb 1.628687\phc 2 \sla
  11  &  \phb 2.409915\phc 0 &  \phb 1.973427\phc 1 &  \phb 1.128215\phc 1 &  \phb 1.351127\phc 1 &  \pha 1.145088\phc 1 \sla
  12  &  \pha 2.284705\phc 4 &  \pha 1.678095\phc 4 &  \pha 1.425455\phc 4 &  \pha 1.430836\phc 4 &  \pha 1.057438\phc 2 \sla
  13  &  \phb 1.626040\phc 3 &  \phb 2.578611\phc 4 &  \phb 2.312520\phc 4 &  \phb 2.309884\phc 4 &  \pha 2.889762\phc 3 \sla
  14  &  \pha 3.061327\phc 1 &  \pha 1.482123\phc 1 &  \pha 9.305564\phc 0 &  \pha 5.903307\phc 0 &  \pha 1.382858\phc 4 \sla
  15  &  \phb 1.154976\phc 1 &  \phb 4.890915\phc 0 &  \phb 2.464009\phc 0 &  \phb 7.783954\phd 1 &  \phb 1.175497\phc 2 \sla
  16  &  \pha 1.177836\phc 0 &  \pha 4.749030\phd 1 &  \pha 1.985955\phd 1 &  \phb 6.723144\phd 3 &  \pha 7.886144\phc 1 \sla
  17  &  \phb 3.924553\phc 3 &  \phb 1.338283\phc 2 &  \pha 7.430600\phc 2 &  \pha 8.954974\phc 2 &  \phb 1.777082\phc 1 \sla
  18  &  \phb 1.684643\phd 1 &  \phb 1.164408\phc 0 &  \phb 1.174242\phc 0 &  \phb 1.156336\phc 0 &  \pha 1.338843\phc 0 \sla
  19  &  \pha 7.329367\phd 1 &  \pha 8.227443\phd 1 &  \pha 7.502286\phd 1 &  \pha 8.040080\phd 1 &  \pha 3.406424\phc 3 \sla
  20  &  \phb 7.165353\phd 1 &  \pha 5.864073\phd 1 &  \pha 2.358848\phd 1 &  \pha 1.464416\phd 1 &  \emp \sla
  21  &  \phb 4.441274\phd 1 &  \phb 2.056313\phc 0 &  \pha 2.937507\phc 0 &  \pha 2.615630\phc 0 &  \emp \sla
\hline
\end{tabular}
\end{center}
%\tablenotetext {a}{units: \cms}
%\tablenotetext {b}{units: \ecps}
%\tablenotetext {c}{rms, and maximum and minimum of signed deviations in \%}
%\tablenotetext {d}{coefficients $\alpha_i$ for analytic fits}
\noindent
$^a$units: \cms; $^b$units: \ecps; $^c$rms and maximum and minimum of
signed deviations in \%; $^d$coefficients $\alpha_i$ for analytic fits
\end {table}

\clearpage

\newpage

%REF
\begin{references}
\bigskip
%definitions for journals not defined is aasms.sty
\def\cpl{Chem.\ Phys.\ Letters}
\def\ijck{Int.\ J.\ Chem.\ Kinetics}
\def\jcp{J.\ Chem.\ Phys.}
\def\jms{J.\ Mol.\ Spectrosc.}
\def\jpc{J.\ Phys.\ Chem.}
\def\jpcrd{J.\ Phys.\ Chem.\ Ref.\ Data}

\reference Black, J. H., \& Dalgarno, A. 1976, \apj, 203, 132

\reference Black, J. H., Porter, A., \& Dalgarno, A. 1981, \apj, 249, 138

\reference Blais, N.\ C., \& Truhlar, D.\ G.\ 1979, \jcp, 70, 2962
  
\reference Blais, N.\ C., Truhlar, D.\ G., \& Mead, C.\ A.\ 1988, \jcp,
89, 6204

\reference Boothroyd, A.\ I., Keogh, W.\ J., Martin, P.\ G., \&
Peterson, M.\ R.\ 1991, \jcp, 95, 4343,

\reference Boothroyd, A.\ I., Keogh, W.\ J., Martin, P.\ G., \&
Peterson, M.\ R.\ 1995, \jcp, in preparation

\reference Breshears, W. D., \& Bird, P. F. 1973, in 14th Symposium
(Int.) on Combustion (Pittsburgh: Combustion Institute), 211

\reference Burton, M. G., Hollenbach, D., \& Tielens, A.\ G. G. M. 1990,
\apj, 365, 620

\reference Chang, C.\ A., \& Martin, P.\ G.\ 1991, \apj, 378, 202

\reference Dalgarno, A., Black, J. H., \& Weisheit, J. C. 1973,
Astrophysical Lett., 14, 77

\reference Dalgarno, A., \& Roberge, W.\ G.\ 1979, \apj, 233, L25

\reference Dove, J.\ E., \& Mandy, M.\ E.\ 1986a, \apj, 311, L93

\reference Dove, J.\ E., \& Mandy, M.\ E.\ 1986b, \ijck, 18, 993

\reference Dove, J. E., Rusk, A. C. M., Cribb, P. H., \& Martin,
P. G. 1987, \apj, 318, 379

\reference Duley, W. W., \& Williams, D. A. 1986, \mnras, 223, 177

\reference Elitzur, M., \& Watson, W.\ D.\ 1978, \aap, 70, 443

\reference Flower, D., \& Watt, G.\ D.\ 1984, \mnras, 209, 25

\reference Fox, L., \& Wilkinson, J. H. 1990, The NAG Fortran Library
Manual, Mark 14 (Oxford: NAG Ltd.)

\reference Hollenbach, D., \& McKee, C. F. 1979, \apjs, 41, 555

\reference Hollenbach, D., \& McKee, C. F. 1989, \apj, 342, 306

\reference Keogh, W.\ J., Boothroyd, A.\ I., Martin, P.\ G., Mielke, S.\
L., Truhlar, D.\ G., \& Schwenke D.\ W.\ 1992, \cpl, 195, 144

\reference Lepp, S., Buch, V.  \& Dalgarno, A. 1995, \apjs, in press
%not j,f,march

\reference Lepp, S., \& Shull, J.\ M.\ 1983, \apj, 270, 578

\reference LeRoy, R.\ J., \& Bernstein, R.\ B.\ 1971, \jcp, 54, 5114

\reference Liu, B.\ 1973, \jcp, 58, 1925

\reference Mac~Low, M., \& Shull, J. M. 1986, \apj, 302, 585

\reference Mandy, M.\ E., \& Martin, P.\ G.\ 1991, \jpc, 95, 8726

\reference Mandy, M.\ E., \& Martin, P.\ G.\ 1992, \jcp, 97, 265

\reference Mandy, M.\ E., \& Martin, P.\ G.\ 1993, \apjsupp, 86, 199

\reference Mandy, M.\ E., \& Martin, P.\ G.\ 1995, \jcp, in preparation

\reference Martin, P.\ G., \&  Mandy, M.\ E.1995, \apjl, xx, Lyy (submitted) %letters

\reference Nishimura, S. 1968, Ann. Tokyo Astron. Obs. Ser. 2, 11, 33

\reference Roberge, W., \& Dalgarno, A.\ 1982, \apj, 255, 176

\reference Schofield, K. 1967, Planet.\ Space Sci., 15, 643

\reference Schwartz, C., \& LeRoy, R.\ J.\ 1987, \jms, 121, 420

\reference Schwenke, D.\ W. 1990, \jcp, 92, 7267

\reference Siegbahn, P., \& Liu, B.\ 1978, \jcp, 68, 2457

\reference Shapiro, P. 1992, in Astrochemistry of Cosmic Phenomena, ed.\
P. D. Singh (Dordrecht: Kluwer), 73

\reference Sternberg, A., \& Dalgarno, A.\ 1989, \apj, 338, 197

\reference Sternberg, A., \& Dalgarno, A.\ 1995, \apjsupp, 99, 565

\reference Sun, Y., \& Dalgarno, A.\ 1994, \apj, 427, 1053

\reference Takayanagi, K., \& Nishimura, S.\ 1960, \pasj, 12, 77

\reference Takayanagi, K., Sakimoto, K., \& Onda, K. 1987, \apjl, 318,
L81

\reference Tanaka, M., Hasegawa, T., Hayashi, S. S., Brand, P. W. J. L.,
\& Gatley, I.  1989, \apj, 336, 207

\reference Tielens, A.\ G. G. M., \& Hollenbach, D.\ M.\ 1985, \apj,
291, 722 

\reference Truhlar, D.\ G., \& Horowitz, C.\ J.,\ 1978, \jcp, 68, 2466
and 71, 1514E

\reference Turner, J., Kirby-Docken, K., \& Dalgarno, A.\ 1977,
\apjsupp, 35, 281

\reference Varandas, A.\ J.\ C., Brown, F.\ B., Mead, C.\ A., Truhlar,
D.\ G., \& Blais, N.\ C.,\ 1987, \jcp, 86, 6258

\reference Wagenblast, R., \& Hartquist, T. W. 1989, \mnras, 237, 1019

\reference Wu, Y.\ M., \& Kuppermann, A.\ 1993, \cpl, 201, 178

\end{references}

\newpage

\centerline{\bf Figure Captions}

\noindent Fig.~1. -- The \vj\ energy levels of the \h2\ molecule in its
ground electronic state.  The dissociation energy is indicated by the
dashed line.

\noindent Fig.~2. -- The \vj\ energy levels of the \h2\ molecule
contoured as functions of $v$ and $j$.  Negative energies (bound states)
are represented by the dashed lines; the 14 contours correspond to the
energies of 14 excited $(v,0)$ states. The first solid contour
corresponds to the dissociation energy.  The other solid contours are
for positive energies for a few of the many quasibound states: (0,32),
(2,30), (1,32), (0,34), (1,34).

\noindent Fig.~3a. -- (a) Contour plot of the steady-state fractional
populations of \vj\ states of \h2\ at $T=2000$~K and $n_{\rm
H}=10^9$~\cm3.  Contours are at 1~dex intervals, decreasing from a value
$10^{-1}$ in the lower left.  Steep falloff in populations near the
dissociation limit arises from collision induced dissociation and, for
quasibound states, from dissociative tunneling as well.
%
(b) Same as (a), but for $n_{\rm H}=10^3$~\cm3. The effect of the cap on
the tunneling rate (\S~2.4) can be seen in the overpopulations of the
highest energy quasibound states.

\noindent Fig.~4 -- (a) Contour plot of the probabilities of \vj\ states
depopulating by quadrupole radiative emission.  The first set of 13
contours are at 0.1~dex intervals, decreasing from a value $5 \times
10^{-6}$~s$^{-1}$ for the contour encircling the maximum at (0,27) and
the local maximum at (9,3) [the contour passing near (5,15) dividing
these two maxima is at a value $2.5 \times 10^{-6}$~s$^{-1}$].
Thereafter, in the lower left as the probabilities continue to decrease,
at contour level $-6.5$~dex the contour interval increases to 0.5~dex.
%
(b) Contour plot of the rate coefficient for \vj\ states depopulating by
collisional energy transfer and dissociation at $T=2000$~K.  Contours
are at 0.1~dex intervals increasing from a value $5 \times
10^{-11}$~\cms\ at the lower left to $3.2 \times 10^{-9}$~~\cms\ at the
upper left.
%
(c) Contour plot of critical density of \vj\ states at $T=2000$~K,
formed from the ratio of the values in Fig.~4a and 4b (equation~5).  The
first 14 contours are at 0.1~dex intervals decreasing away from a value
$6.3 \times 10^{3}$~\cm3\ for the closed contour around (2,5).
Subsequently in the lower left at contour level 2.5~dex the contour
interval increases to 0.5~dex.  Along $v=0$ the critical density rises
rapidly from a value $6.2 \times 10^{-1}$~\cm3\ for (0,2) to a maximum
of $1.2 \times 10^{4}$~\cm3\ for (0,17).

\noindent Fig.~5. -- Steady-state rate coefficient $\gamma_{po}$ for
conversion of ortho to para \h2\ (and vice versa) as a function of
collision partner density $n_{\rm H}$. Values for constant $T$ are
joined by curves whose analytic form is given in the Appendix. Beginning
at the lower part of the diagram, the 15 values in the $T$ grid are (in
K): 450, 600, 780, 1000, 1400, 2000, 3000, 4500, 6000, 7800, 10000,
14000, 20000, 30000, 45000; the systematic order of different line types
(and symbols) helps in deciphering $T$ in the upper part of the diagram
where lines overlap.

\noindent Fig.~6. -- Like Fig.~5 but for ratio of ortho to para \h2\ (no
analytic fit has been performed).  Note that $T$ grid begins with 450~K
at the upper part of the diagram.

\noindent Fig.~7. -- Like Fig.~5 but for rate coefficient $\gamma_{cid}$
for collision induced dissociation.

\noindent Fig.~8. -- Like Fig.~5 but for rate coefficient $\gamma_{dt}$
for dissociative tunneling.

\noindent Fig.~9. -- (a) Contour plot of the CID rate coefficient
$\gamma_{349,k}$ (for ortho states of \h2) at $T=2000$~K.  Contours are
at 0.5~dex intervals, decreasing from a value $10^{-9}$~\cms\ at the
upper left to $10^{-23}$~\cms\ at the lower left.
%
(b) Contour plot of the population-weighted contributions
$\gamma_{349,k}p_k + r_{349,k}p_k/n_H$, for \nh\ = $10^{3}$~\cm3\ (see
Fig.~3b for the fractional population $p_k$).  Contours are at 0.5~dex
intervals. The highest CID contributions are from the states (14,3) and
(13,5); the contour level just below these states is $10^{-23.5}$~\cms.
Contours at $10^{-24}$~\cms\ are adjacent to this and also above the
local maximum at (0,5).  At $10^{-24.5}$~\cms\ contours loop down
through (9,9) and rise from (0,15); contours at $10^{-25}$~\cms\ loop
down through (9,3) and rise from (0,23). Although the CID contribution
from the sparsely populated quasibound states is very low, the
dissociative tunneling is important. The highest contours,
$10^{-23}$~\cms, surround the states with highest dissociative tunneling
contribution, (0,33) and (6,23).

\noindent Fig.~10. -- Rate coefficient $\gamma_{ie}$ for the cooling
contribution related to decreasing the internal energy (often called
dissociative cooling for steady-state conditions; see equation~12 and
also Fig.~7 for $\gamma_{cid}$ and Fig.~8 for $\gamma_{dt}$).  Lines
types and symbols like Fig.~5, except that results are not shown for the
lowest temperatures in our grid ($T < 2000$~K).  On our grid
$\gamma_{ie}$ is important relative to $\gamma_r$ only for temperatures
3000~K and above, and even then only for \nh\ above the critical density
(see Fig.~12).

\noindent Fig.~11 -- Contour plot of the radiative energy loss rate from
individual \vj\ states.  The first set of 29 contours are at 0.1~dex
intervals, decreasing from a value $6.3 \times 10^{-18}$~erg~s$^{-1}$
for the contour encircling the maximum at (9,7); the contour passing
near (4,19) leading to a local maximum at (2,29) is at a value $2.5
\times 10^{-18}$~s$^{-1}$].  Thereafter, in the lower left, the contour
interval increases to 0.5~dex. The minimum loss rate is $2.1 \times
10^{-24}$~erg~s$^{-1}$ at (0,3).

\noindent Fig.~12. -- Rate coefficient $\gamma_r$ for radiative cooling
in same format as Fig.~10 for $\gamma_{ie}$. Lines types and symbols
like Fig.~5, except that lowest locus is for 600~K. For comparison with
Fig.~10, 3000~K is the second lowest dotted locus; it sinks below that
for $\gamma_{ie}$ at log(\nh)$\sim 7.5$.

\noindent Fig.~13. -- Rate coefficient $\gamma_{cid}$ for collision
induced dissociation as a function of $T$ for fixed densities.  Circles:
our values for \nh\ = $10^{-3}$~\cm3\ (lowest curve), $10^{3}$~\cm3\
(near the critical density), and $10^{9}$~\cm3.  Solid curves: our
analytic fit (equation A1). Dashed curves: analytic form from Lepp \&
Shull (1983) for same densities.  Dash-dotted curves: analytic form for
two lowest densities using low density limit from Mac~Low \& Shull
(1986), experimental high density limit from Breshears \& Bird (1973)
(same as upper dashed curve from Lepp \& Shull) and interpolation
formula of Lepp \& Shull.  Dotted curve: like dash-dotted case for \nh\
= $10^{3}$~\cm3, but with critical density lowered by factor 10.

\noindent Fig.~14. -- Radiative cooling rate per \h2\ molecule,
$\Lambda_r/n_{{\rm H}_2}$ as a function of $T$ for fixed densities: \nh\
= $10^{0}$~\cm3\ (lowest curve), $10^{2}$~\cm3\, $10^{2}$~\cm3\, and
$10^{8}$~\cm3\ (``high density'' limit).  Solid curves: our analytic fit
(equation B2). Dashed curves: analytic form from Lepp \& Shull (1983).
Dash-dotted curves: analytic form from Hollenbach \& McKee (1979).

\end{document}
