\documentclass[11pt]{article}       
\usepackage{a4wide,amsmath,epsf}
\usepackage{graphicx}

\setlength{\textheight}{25cm}
\setlength{\textwidth}{15cm}
\oddsidemargin=0mm
\evensidemargin=0mm
%\setlength{\parindent}{0mm}
\setlength{\itemsep}{-5pt plus-5pt minus-5pt}
\setlength{\parskip}{4pt plus1pt minus1pt}
\setlength{\parsep}{2pt plus1pt minus1pt}
\voffset -1.5cm

    \def\abl#1#2{\frac {d #1}{d #2}}
    \def\pabl#1#2{\frac {\partial #1}{\partial #2}}
    \def\eg{e.\,g.\ }
    \def\ie{i.\,e.\ }
    \def\Tg{T_{\rm g}}
    \def\Trad{T_{\rm rad}}

\begin{document}

\sloppy

\centerline{\huge\bf Recombination Cooling in ProDiMo}
{\ }\\[-1.7ex]
\centerline{\large Peter Woitke, May 2009}

\section{Photo-processes and related heating and cooling rates}

\noindent Consider a pair of photo-process as
\begin{equation}
  {\rm AB} \,+\, h\nu \;\;\begin{array}{c} {\textstyle k_f} \\[-0.6ex]
                                           \rightleftharpoons \\[-0.6ex]
                                           {\textstyle k_r}  
                          \end{array}\; {\rm A \,+\, B}
  \label{reac}
\end{equation}
As main application in mind, we will consider the photo-ionisation
and direct recombination of H$^-$ (AB=H$^-$, A=H, B=e$^-$), but the
following method may also be applicable to a pair of
photo-dissociation and radiative association processes (?). The
forward rate coefficient $k_{\rm f}$ is given by
\begin{equation}
  k_f = 4\pi\!\int\limits_{\nu_{\rm thr}}^\infty
        \!\frac{J_{\nu}}{h\nu}\,\sigma^{\rm f}(\nu)\,d\nu
  \label{kf} \ ,
\end{equation}
such that $n_{\rm AB}\,k_f\;\rm[cm^{-3} s^{-1}]$ is the number of
photoionisations per volume and time. $\sigma^{\rm f}(\nu)$ is the
photo cross section [$\rm cm^2$] of AB for the forward reaction (only
continuous cross sections considered here).  Energy conservation
implies that
\begin{equation}
  h\nu = E_b + E_{\rm th}
\end{equation}  
where $E_b=h\nu_{\rm thr}$ is the binding energy of AB, $\nu_{\rm
thr}$ is the threshold frequency, and $E_{\rm th}$ is the excess thermal
energy $\frac{1}{2}m_e v^2$ liberated by the photo-reaction. 

Hence, the {\sl ionisation heating rate} $\Gamma
\rm[erg\,s^{-1}cm^{-3}]$ is given by
\begin{equation}
  \Gamma = n_{\rm AB}\;4\pi\!\int\limits_{\nu_{\rm thr}}^\infty \,J_{\nu}\, 
           \frac{\nu-\nu_{\rm thr}}{\nu}\,\sigma^{\rm f}(\nu)\,d\nu \ .
  \label{heat}
\end{equation}
The reverse process, the recombination cooling rate, is slightly
more difficult to understand.  A direct recombination destroys
one thermal electron, so the cooling rate should be of order
$n_{\rm A} n_{\rm B} k_r(\Tg) \langle E_{\rm th} \rangle$, where
the recombination rate coefficient is given by an Arrhenius
law with fit coefficients $\alpha$, $\beta$, $\gamma$ as
\begin{equation}
  k_r(\Tg) = \alpha\,\Big(\frac{\Tg}{300\rm \,K}\Big)^{\beta} 
             \exp\left(-\frac{\gamma}{\Tg}\right)
  \label{kr} \ .
\end{equation}
In order to investigate the details of $\langle E_{\rm th} \rangle$,
we consider the pair of photo-processes (\ref{reac}) in the case of
thermodynamical equilibrium, where $J_\nu=B_\nu(\Tg)$ and where every
reaction is in detailed balance with its direct reverse, \ie
\begin{equation}
  n^\star_{\rm AB} k_f\big\vert_{J_\nu=B_\nu(\Tg)} 
  = n^\star_{\rm A} n^\star_{\rm B} k_r(\Tg) \ .
  \label{equil}
\end{equation}
Thus, we find the recombination rate coefficient to be
\begin{equation}
  k_r(\Tg) = \frac{n^\star_{\rm AB}}{n^\star_{\rm A} n^\star_{\rm B}}\;
             4\pi\!\int\limits_{\nu_{\rm thr}}^\infty
             \!\frac{B_{\nu}(\Tg)}{h\nu}\,\sigma^{\rm f}(\nu)\,d\nu
  \label{kr_Milne} \ ,
\end{equation}
which is known as the bound-free Milne relation. Equation
(\ref{kr_Milne}) allows us to calculate the {\sl recombination cooling rate}
$\Lambda \rm[erg\,s^{-1}cm^{-3}]$ analogous to Eq.\,(\ref{heat})
\begin{equation}
  \Lambda = n_{\rm A} n_{\rm B}\;
            \frac{n^\star_{\rm AB}}{n^\star_{\rm A} n^\star_{\rm B}}\;
            4\pi\!\int\limits_{\nu_{\rm thr}}^\infty
            B_{\nu}(\Tg)\,\frac{\nu-\nu_{\rm thr}}{\nu}
            \,\sigma^{\rm f}(\nu)\,d\nu
  \label{cool} \ ,
\end{equation}
Obviously, $\Gamma\!=\!\Lambda$ in thermodynamical equilibrium, as
required from first principles. The fraction $n^\star_{\rm
AB}/(n^\star_{\rm A} n^\star_{\rm B})$ can, in principal, be determined
from thermochemical data.
\begin{equation}
  p^\star_{\rm AB} = p^\star_{\rm A}\,p^\star_{\rm B}\,k_p(\Tg)
\end{equation}
which relates the partial pressures in thermodynamical equilibrium.
$k_p(\Tg)$ is the equilibrium constant of H$^-$. Figure~\ref{fig:kr_Milne}
shows an attempt to use Eq.\,(\ref{kr_Milne}) directly to calculate
the recombination rate coefficient $k_r$ which clearly shows some
numerical problems at low temperatures where the integral goes to
zero and the fraction $\frac{n^\star_{\rm AB}}{n^\star_{\rm A}
  n^\star_{\rm B}}$ goes to infinity. Also, the range of validity of 
the polynomial fit of $k_p$ is clearly too narrow.

\begin{figure}
  \centering
  \includegraphics[height=8.5cm]{kr_Milne.eps}
  \caption{Comparison of the recombination rate coefficient $k_r(\Tg)$ of
  H$^-$ according to Eq.\,(\ref{kr}) with data from UMIST 2006 (full
  line) and according to Eq.\,(\ref{kr_Milne}) with photo
  cross-sections data from Leiden's photoionisation database and
  thermo-chemical data from the JANAF tables (1985), involving
  a polynomial fit of $k_p(\Tg)$ between 1000\,K and 6000\,K.} 
   \label{fig:kr_Milne}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[height=8.5cm]{Eth_over_kTg.eps}
  \caption{Mean thermal energy $\langle E_{\rm th} \rangle$ destroyed
    by a direct recombination of $\rm H + e^- \rightarrow H^- + h\nu$ in
    units of $k\Tg$, according to Eq.~(\ref{cool2}).}
\end{figure}

Therefore, we will rely on Eq.\,(\ref{kr}) rather than on
Eq.\,(\ref{kr_Milne}) to compute the recombination rate coefficient
$k_r$, also because these rates are deeply imbedded in ProDiMo's treatment
of chemistry. However, we will indirectly use Eq.\,(\ref{kr_Milne}),
in form of Eq.\,(\ref{cool}), to get a handle on $\langle E_{\rm th}
\rangle$.  Inserting Eq.\,(\ref{equil}) into Eq.\,(\ref{cool})
results in 
\begin{equation}
  \Lambda = \frac{n_{\rm A}\,n_{\rm B}\,k_r(\Tg)}
                 {k_f\big\vert_{J_\nu=B_\nu(\Tg)}}\;
            4\pi\!\int\limits_{\nu_{\rm thr}}^\infty
            B_{\nu}(\Tg)\,\frac{\nu-\nu_{\rm thr}}{\nu}
            \,\sigma^{\rm f}(\nu)\,d\nu
  \label{cool2} \ ,
\end{equation}
where $k_r$ is computed from Eq.\,(\ref{kr}).  ProDiMo uses directly
Eq.~(\ref{cool2}) to calculate recombination cooling rates. Figure 2
shows the results for $\rm H + e^- \rightarrow H^- + h\nu$, with the
cross section data from Fig.~3, in units of $n_{\rm A} n_{\rm B}
k_r(\Tg) k\Tg$, showing that the mean thermal energy destroyed by a
direct recombination is roughly given by $\langle E_{\rm th}
\rangle\approx (1.5-2.5)k\Tg$.

\begin{figure}
  \centering
  \includegraphics[height=8.5cm]{Hmin-cross.eps}
  \caption{Bound-free ionisation cross-section $\sigma^f(\lambda)$ of
  the negative hydrogen ion H$^-$ from Leiden's photoionisation
  database. Binding energy of H$^-$ is $E_b=0.755\,eV$, which
  corresponds to threshold wavelength $\lambda_{\rm thr}=1.642\mu$m. }
\end{figure}

\section{How does recombination cooling work?}

Equations (\ref{heat}) and (\ref{cool2}) are sufficient to calculate
ionisation heating and recombination cooling rates, and ProDiMo uses
them directly. However, these equations need some interpretation to
understand what's going on.  In kinetic chemical equilibrium, free
electrons are destroyed and created at the same rate, and whether
H$^-$ bound-free and free-bound transitions cause radiative heating or
cooling is given by the {\sl difference between the thermal energies
of created and destroyed free electrons}.

Let's assume that reaction (1) is in kinetical equilibrium, \ie
$n_{\rm AB} k_f = n_{\rm A} n_{\rm B} k_r$ (other reaction channels
like $\rm H^- + H \leftrightarrow H_2 + e^-$ are ignored in this
section). What's the net (heating$\,-\,$cooling) rate 
$Q_{\rm net}=\Gamma-\Lambda$?
\begin{equation}
  Q_{\rm net} = n_{\rm AB}\;4\pi\!\int\limits_{\nu_{\rm thr}}^\infty
                J_\nu \frac{\nu-\nu_{\rm thr}}{\nu}
                \,\sigma^{\rm f}(\nu)\,d\nu
             \;-\; n_{\rm A} n_{\rm B} k_r \langle E_{\rm th}\rangle
\end{equation}
Using the assumption of detailed kinetic equilibrium we find
\begin{equation}
  Q_{\rm net} = n_{\rm AB} k_f 
   \Big( \langle E_{\rm abs}\rangle - \langle E_{\rm th}\rangle \Big)
\end{equation}
where
\begin{equation}
  \langle E_{\rm abs}\rangle = \frac
   {\int\limits_{\nu_{\rm thr}}^\infty
    J_\nu \frac{h\nu-E_b}{h\nu} \,\sigma^{\rm f}(\nu)\,d\nu}
   {\int\limits_{\nu_{\rm thr}}^\infty
    J_\nu \frac{1}{h\nu} \,\sigma^{\rm f}(\nu)\,d\nu}
   = \frac
   {\int\limits_{\nu_{\rm thr}}^\infty
    J_\nu \,\sigma^{\rm f}(\nu)\,d\nu}
   {\int\limits_{\nu_{\rm thr}}^\infty
    J_\nu \frac{1}{h\nu} \,\sigma^{\rm f}(\nu)\,d\nu} - E_b
   \label{Erad}
\end{equation}
is the mean excess photon energy $\langle h\nu\rangle-E_b$ as weighted
by the photo-cross section in a given radiation field $J_\nu$.  It is
instructive to consider a dilute Planckian radiation field
$J_\nu=W\,B_\nu(\Trad)$ in the Wien limit ($h\nu\gg k\Trad$), \ie
$J_\nu\approx W\,\frac{2h\nu^3}{c^2}\exp(-h\nu/k\Trad)$. For this
case, we can derive
\begin{eqnarray}
  \pabl{k_f}{\Trad} &=& \pabl{}{\Trad} 4\pi 
   \int\limits_{\nu_{\rm thr}}^\infty
    W\,\frac{2h\nu^3}{c^2}\exp\Big(-\frac{h\nu}{k\Trad}\Big) 
    \frac{1}{h\nu} \,\sigma^{\rm f}(\nu)\,d\nu   \\
   &=& \frac{4\pi}{k\Trad^2} \int\limits_{\nu_{\rm thr}}^\infty
    W\,\frac{2h\nu^3}{c^2}\exp\Big(-\frac{h\nu}{k\Trad}\Big) 
    \,\sigma^{\rm f}(\nu)\,d\nu 
\end{eqnarray}
Identifying the terms in Eq.~(\ref{Erad}), we find
\begin{equation}
  \langle E_{\rm abs}\rangle 
   = \frac{k\Trad^2\pabl{k_f}{\Trad}}{k_f} - E_b
\end{equation}

\end{document}