text:
\section{Introduction}

High-order harmonic generation (HHG) is one of the most studied
nonlinear phenomena in intense laser-matter interaction
\cite{brabec00}. HHG is most easily understood using the
three-step model \cite{kulander,corkum,lewenstein} -- first the
electron is released by tunnel ionization; second, it is
accelerated by the oscillating electric field of the laser and
later driven back to the target ion; and third, the electron
recombines with the ion to emit a high energy photon. A
semiclassical formulation of the three-step model based on the
strong-field approximation (SFA) is given by Lewenstein {\it et
al} \cite{lewenstein}. In this model (often called the Lewenstein
model), the liberated continuum electron experiences the full
effect from the laser field, but the effect of the target
potential is neglected. Despite this limitation, the Lewenstein
model has been widely used for studying HHG from atoms and
molecules, and for characterizing the nature of the harmonics
generated. Since the continuum electron needs to come back to
revisit the parent ion in order to emit radiation, the neglect of
the electron-ion interaction in the SFA model is rather
questionable. In the past years, various efforts have been made to
improve upon the SFA model, i.e., by including the Coulomb
distortion \cite{ivanov,kaminski,ciappina,smirnova}, or by using
Ehrenfest theorem \cite{gordon}. All of these improvements can
still be considered insufficient since they fail to include
accurately the scattering of the continuum electrons with the
parent ions.

With moderate effort, direct numerical solution of the
time-dependent Schr\"odinger equation (TDSE) for atomic targets
can be accurately carried out, at least within the single active
electron (SAE) approximation. The TDSE methods are computationally
more demanding for molecular targets. Even for simplest diatomic
molecules, one has to solve a three-dimensional time-dependent
Schr\"odinger equation provided the molecular axis has an
arbitrary orientation with respect to the polarization of the
laser field. Such calculations require fast CPU and large memory
computers and have been accomplished only recently
\cite{lein03,kamta,telnov07}. In order to compare with
experimental results, thus most of the existing calculations for
HHG from molecules were carried out using the SFA model
\cite{zhou,atle06,hoang}. Based on the experience from high
harmonics generated by atomic targets, such SFA model is not
expected to offer accurate predictions. The lack of a reliable
theory for describing HHG from molecules has prevented the
possible exploitation of molecular structure using infrared laser
pulses. This is unfortunate since it has been well recognized that
infrared lasers have the potential for probing time-resolved
molecular dynamics, with temporal resolutions down to sub-ten
femtoseconds.

 In view of the difficulty in solving the TDSE as a practical theoretical tool  for
 predicting the nonlinear interactions between molecules with
 intense lasers, recently Le {\it et al} \cite{atle08} and Morishita {\it et
 al} \cite{toru08} have developed an alternative approach for
 calculating HHG, and have tested the model for atomic targets. The theory is based on the
 three-step model but without the approximations that lead to the SFA (or the Lewenstein model).
 The HHG is viewed as resulting from the photo-recombination of
 the returning electrons generated by the laser pulses. The HHG yield is
 expressed as the product of a returning wave packet with the
 photo-recombination cross section (PRCS). The shape of the wave packet has been
 shown to be largely independent of the targets for the same laser
 pulse and can thus be obtained from a reference atom or from the SFA. Thus there are
 two ways to obtain HHG spectra without
 actually solving the TDSE. The first method is called SW-SFA
 model where the HHG spectra are obtained by replacing the dipole
 matrix element using scattering waves for the continuum electrons
 instead of the plane waves used in the SFA, but retaining the
 electron wave packet from the SFA. The other method is to compare
 the HHG spectra with another reference atom where   the HHG
 spectra have been obtained by accurate calculations or by
 experiment. If the accurate PRCS for both the reference atom and
 the target are known, then the HHG spectra can be
 obtained from the HHG spectra of the other atom. The validity of
 the two methods have been established in Le {\it et al} \cite{atle08} and Morishita {\it et
 al} \cite{toru08} for atomic systems where the results from the
 models are compared to accurate results from the TDSE. While the
 model has been tested only for atoms so far, there is a general
 belief that the same models would  apply to molecular targets.
 In this paper our goal is to apply the models to the H$_2^+$ molecules
 where accurate HHG spectra have been calculated by solving the
 TDSE \cite{lein03,kamta,telnov07}. The results from this comparison confirm
 that the models indeed work well and the method can be further extended
 as a general tool for calculating HHG from molecular targets.

  First let us be specific about the two methods for obtaining HHG
  spectra. The first method is called the scattering wave based strong-field approximation
(SW-SFA) \cite{atle08}. The  HHG yield $S(\theta,\omega)$ for a
molecule whose internuclear axis makes an alignment angle $\theta$
with respect to the laser polarization direction can be obtained
from the SFA result $S^{SFA}(\theta,\omega)$ by
\begin{equation}
S(\theta,\omega)=\left|\frac{T(\theta,E)}{T^{PWA}(\theta,E)}\right|^2
S^{SFA}(\theta,\omega) \label{SW-SFA}
\end{equation}
where $T$ and $T^{PWA}$ are the {\em exact} transition dipole
matrix element and its plane-wave approximation (PWA),
respectively. Here the electron energy $E$ is related to the
emitted photon energy   $\omega$ by $E=k^2/2=\omega-I_p$, with
$I_p$ being the ionization potential of the target (atomic units
are used throughout the paper unless otherwise indicated). As only
the ratio of the transition dipoles enters Eq.~(1), one can
formulate the SW-SFA model by using the above equation with the
photoionization or its time-reversed process, i.e., the
photo-recombination cross sections. In this paper we refer to the
photoionization process, as it is more widely available
theoretically and experimentally. Note that, one can identify, up
to some factor, $S^{SFA}/|T^{PWA}|^2$  as the flux of the
returning electron, which we will call a ``wave-packet''. The
returning electron that contributes most to the HHG is the one
that propagates along the laser polarization direction. Therefore,
the relevant differential cross section is for ${\bi k}$ that is
parallel to the laser polarization axis, that is, for $\theta_{\bi
k}=0$ and $\pi$. Note that for asymmetric molecules, the two
contributions from $\theta_{\bi k}=0$ and $\pi$ differ, in
general.

Clearly, the formulation of the SW-SFA given in Eq.~(1) is not unique. In fact, one
can also use the ``wave-packet'' from a reference atom
(i.e., atom with the same $I_p$ as for the target). In other words, the HHG yield
can also be written as
\begin{equation}
S(\theta,\omega)=\left|\frac{T(\theta,E)}{T^{ref}(E)}\right|^2
S^{ref}(\omega) \label{scaledH}
\end{equation}
where the superscript ``{\it ref}'' refers to the reference atom.
In the following, we will check the validity of both models by
comparing their predictions against the benchmark data. For the
reference atom, we will use a scaled hydrogen with the effective
nuclear charge chosen such that it has the same $1s$ binding
energy as H$_2^+$. The advantage of using  the scaled hydrogen as
a reference atom is that one can accurately and efficiently
calculate both the HHG yield by solving the TDSE numerically, and
the exact photo-recombination cross section analytically.
Experimentally one can also replace the scaled atomic hydrogen
with an atomic target with comparable ionization energy.

In this paper we will focus only on the HHG component in the
parallel polarization direction. The theoretical methods used to
generate the data for the present study have all been given
previously. The ``benchmark'' HHG spectra of H$_2^+$ are calculated
by solving the TDSE using the methods given in
\cite{lein03,telnov07}. The calculation of the transition
dipole matrix elements are described in \cite{picca06,picca07}
while the HHG calculated within the SFA is given in
\cite{zhou,atle06,hoang}. In terms of computational effort and
thus the accuracy, the solution of the TDSE is the most demanding.

\begin{figure}
\centering
\mbox{\rotatebox{0}{\myscaleboxb{
\includegraphics{Fig1a.eps}}}}
\mbox{\rotatebox{0}{\myscaleboxb{
\includegraphics{Fig1b.eps}}}}
\caption{Comparison of the photoionization differential cross sections from the
exact (i.e.,with scattering wave) and the plane-wave approximation
calculations, shown in (a) and (b), respectively, for alignment
angles $\theta=0^{\circ}$, $30^{\circ}$, $45^{\circ}$ and
$60^{\circ}$.} \label{fig1}
\end{figure}

In Fig.~1(a) we show the  photoionization differential cross
sections (DCS) (electrons moving along the polarization axis), as
functions of emitted electron energy, for H$_2^+$ at the
equilibrium distance $R=2.0$ a.u., alignment angles
$\theta=0^{\circ}$, $30^{\circ}$, $45^{\circ}$ and $60^{\circ}$.
The data are shown for the energy range up to 200 eV, relevant to
harmonic generation from typical infrared lasers.   The noticeable
features of these curves are the pronounced minima which move to
the higher energies, as the alignment angle increases. The minimum
positions are $12$ eV, $36$ eV and $117$ eV, for
$\theta=30^{\circ}$, $45^{\circ}$, and $60^{\circ}$, respectively.
As we will show later, these minima are responsible for the
interference minima, seen in the HHG spectra. The corresponding
data obtained from the PWA are plotted in Fig.~1(b). Apart from
the clear differences in the shape of these curves in comparison
to the exact data, there are strong discrepancies in the positions
of the minima. For $\theta=30^{\circ}$, the minimum is shifted to
higher energy as much as about $100$ eV. For the two larger
angles, the minima are not even seen within the energy range of
the plot. That is the reason that the minima in the HHG spectra
from SFA are shifted to much higher harmonic orders in comparison
to the exact TDSE calculations, as observed by Kamta and Bandrauk
\cite{kamta} and by Chirila and Lein \cite{chirila}.

\begin{figure}
\centering \mbox{\rotatebox{0}{\myscaleboxc{
\includegraphics{Fig2a.eps}}}}
\mbox{\rotatebox{0}{\myscaleboxc{
\includegraphics{Fig2b.eps}}}}
\mbox{\rotatebox{0}{\myscaleboxc{
\includegraphics{Fig2c.eps}}}}
\mbox{\rotatebox{0}{\myscaleboxc{
\includegraphics{Fig2d.eps}}}}
\caption{Comparison of the HHG spectra obtained from the TDSE (solid red lines),
SFA (dashed lines), SW-SFA (solid blue lines), and with the use of scaled hydrogen
atom (solid black lines) for the alignment angle $\theta=0^{\circ}$, $30^{\circ}$,
$45^{\circ}$ and $60^{\circ}$. A 20-cycle laser with peak intensity of
$3\times 10^{14}$ W/cm$^2$, wavelength of 800 nm is used.} \label{fig2}
\end{figure}

Next we compare the HHG spectra calculated using the different
models with those obtained from the solution of the TDSE.  In
Fig.~2(a-d) we show the TDSE (solid red lines) and  SFA (dashed
black lines) results against those from using SW-SFA (solid blue
lines) and with the wave-packet from scaled hydrogen (solid black
lines), for $R=2.0$ a.u. and $\theta=0^{\circ}$, $30^{\circ}$,
$45^{\circ}$, and $60^{\circ}$. We use laser pulse with peak
intensity of $3\times 10^{14}$ W/cm$^2$, wavelength of 800 nm, and
20-cycle duration with a sine-square envelope. The TDSE results
for H$_2^+$ were obtained by using the method described in
Ref.~\cite{telnov07}. Here we have normalized the data near the
cutoff and only the odd harmonics are shown.

For $\theta=0^{\circ}$ and $30^{\circ}$, one can see that the SFA
results disagree greatly with the TDSE results. The improvements by
the SW-SFA and with the wave-packet from the scaled hydrogen are
quite striking. Indeed, results from both models reproduce quite
well the overall shape of the spectra over the broad range down to
H21, i.e., just above the ionization threshold ($I_p=30$ eV). For
$\theta=30^{\circ}$ the minimum in both model calculations near H27
comes from the dip near $E=12$ eV in the exact DCS [see Fig.~1(a)].
This minimum is consistent with the TDSE data  \cite{telnov07},
although not clearly seen in Fig.~1(b). In fact, this is also
consistent with the earlier results by Lein {\it et al}
\cite{lein03} and Kamta and Bandrauk \cite{kamta}. For
$\theta=45^{\circ}$, the minimum in both of our model calculation is
near H43, also agrees well with the prediction from all the
available TDSE results reported in Refs.~\cite{kamta} (see their
Table II), and \cite{telnov07}. Again, this minimum comes from the
dip near $36$ eV in the exact DCS. We note that the minima from our
model calculations here tend to be more pronounced, compared to the
TDSE results. The remaining discrepancies could be due to the
limitation of our model, or because of the possible lack of
convergence in the TDSE results. Finally, for $\theta=60^{\circ}$,
we notice that all the data, including the SFA, agree quite well
with the TDSE results. This is not entirely surprising since the
shape of the photoionization DCS for this angle from the PWA is
fairly close to that of the exact DCS for the relevant electron
energy range below $70$ eV [see Fig.~1(a) and (b)].

\begin{figure}
\centering \mbox{\rotatebox{0}{\myscaleboxb{
\includegraphics{Fig3a.eps}}}}
\centering \mbox{\rotatebox{0}{\myscaleboxb{
\includegraphics{Fig3b.eps}}}}
\caption{Alignment dependence of some selected harmonics from the SW-SFA (a) and the
TDSE (b).} \label{fig3}
\end{figure}

It is important to emphasize that for each new alignment the
returning wave-packet from the SFA is calculated again, whereas
the wave-packet from scaled hydrogen is, of course, alignment
independent. It has been shown \cite{toru08,atle08} that the shape
of the wave-packet largely depends only on the laser field (for
systems with identical ionization potentials). Using scaled
hydrogen, one therefore can get only the {\em shape} of the
wave-packet, as function of energy. To get the absolute magnitude
of the HHG spectra from the latter, one needs to multiply by a
factor, which accounts for the alignment dependent ionization
yield. This can be done, for example, by using the MO-ADK theory
\cite{moadk}.

To better quantify the prediction of the SW-SFA model for
different molecular alignment against the TDSE results, we show in
Fig.~3(a) and (b) the selective harmonics, H25, H35, H45, H55 and
H65, respectively, against the alignment angles. For convenience,
the data were normalized to $1.0$ at $\theta=0^{\circ}$. Again,
the general shape, as well as the minimum position for each
harmonic, are in quite good agreement.

\begin{figure}
\centering \mbox{\rotatebox{0}{\myscaleboxa{
\includegraphics{Fig4.eps}}}}
\caption{Comparison of the HHG spectra obtained from the TDSE (red
line), SW-SFA (blue line), and with the use of the returning
wave-packet from scaled hydrogen (black line), for alignment angle
$\theta=50^{\circ}$. The data have been shifted vertically to show
the details. A laser pulse with peak intensity of $5\times
10^{14}$ W/cm$^2$, wavelength of 780 nm is used. The pulse is
turned on and off over three cycles and kept constant for four
cycles. The arrow indicates the position of the interference
minimum.} \label{fig4}
\end{figure}

We now show detailed comparison of the HHG spectra from different
calculations in Fig.~4 for the alignment angle
$\theta=50^{\circ}$. Here we use 10-cycle laser pulse with peak
intensity of $5\times 10^{14}$ W/cm$^2$ and wavelength of $780$
nm, the same as has been used by Lein {\it et al} \cite{lein03}.
The pulse is turned on and off over three cycles and kept constant
for four cycles. The TDSE results were obtained by the method
described in Ref.~\cite{lein03}. For clarity, we have shifted the
data vertically. Clearly, the SW-SFA data (blue line) are in quite
good agreements with the TDSE results (red line). The results from
our model using the wave-packet calculated by solving the TDSE
equation for scaled hydrogen (black line) agree even better with
the TDSE results. In particular, the minimum near H55, indicated
by an arrow in the figure, is well reproduced by both model
calculations.

\begin{figure}
\centering \mbox{\rotatebox{0}{\myscaleboxa{
\includegraphics{Fig5.eps}}}}
\caption{Projected internuclear separation vs ``effective''
electron wavelength at the minima in the {\em exact}
photoionization cross sections. Data were obtained with $R=1.8$,
$2.0$, $2.5$, $3.0$, and $4.0$ a.u., and for different alignment
angles as shown in the labels. Note that we have used the
``effective'' wavelength according to Lein {\it et al}
\cite{lein02}, see text.} \label{fig5}
\end{figure}

 In recent years, the ``interference minima'' in the HHG spectra
 from molecular targets have been widely discussed in the
 literature \cite{leinPRL,kanai,vozzi,atle06}. As shown by Lein {\it et al}
\cite{lein02}, the minimum positions in HHG spectra from the TDSE
for H$_2^+$ with $R=1.8$, $2.0$ and $2.5$ a.u. agree quite well
with the two-emitter model \cite{leinPRL,lein02}. Specifically,
the minima satisfy the relation $R\cos(\theta)=\lambda^{eff}/2$,
where $\lambda^{eff}$ is the ``effective'' wavelength of the
continuum electron defined such that the ``effective'' wave vector
is $k^{eff}=\sqrt{2\omega}$, with $\omega$ being the energy of the
emitted photon. In other words, the ``effective'' energy is
shifted by $I_p$ with respect to the usual relation, i.e.,
$k=\sqrt{2(\omega-I_p)}$.

 According to the two models presented here, the minima in HHG are
 attributed to the minima in the photoionization DCS. In Fig.~5
 the projected internuclear
separation $R\cos(\theta)$ is plotted vs the ``effective''
electron wavelength at the minima in the photoionization DCS. The
data were obtained  for $R=1.8$, $2.0$, $2.5$, $3.0$ and $4.0$
a.u. and for alignment angles $\theta=30^{\circ}$, $40^{\circ}$,
$50^{\circ}$ $60^{\circ}$ and $70^{\circ}$. The results are indeed
scattered nicely around $R\cos(\theta)=1/2\lambda^{eff}$, shown as
the dashed line, for small $\lambda^{eff}$ and tend to lie above
it for $\lambda^{eff}>3.5$ a.u. This is consistent with the
results from the numerical solution of the TDSE by Lein {\it et
al} \cite{lein02} [see their Fig.~(3)]. This result demonstrates
that the minima in the HHG spectra from H$_2^+$ are fully and
accurately reproducible by the minima in the photoionization cross
sections, and hence, by the SW-SFA model, for other internuclear
distances as well. We note in this connection that analysis of the
recombination cross sections in relation with the minima positions
in the HHG spectra from H$_2$ and N$_2$ has been reported recently
by Zimmermann {\it et al} \cite{zimmermann}.

In this paper we have shown that harmonic generation from H$_2^+$
molecule can be factored out into the product of a returning
electron wave packet with the photo-recombination cross section. The
wave packet can be extracted from the SFA model or from HHG
generated by atoms with similar ionization potential. Using either
wave packet, but replacing the PRCS of either case by the accurate
PRCS of the molecule, we show that the resulting HHG spectra
compared favorably with those from solving the TDSE. Similar
conclusions have been reached for atomic targets recently
\cite{atle08,toru08}. The models have been tested using different
laser parameters and for different internuclear distances, and we
expect the models to work well for molecules in general. Since
photoionization (or photo-recombination) is a linear process, it is
intrinsically much easier to treat theoretically than the nonlinear
laser-molecule interaction. While accurate calculations of molecular
photoionization DCS are by no means trivial, sophisticated packages
have been developed \cite{gianturco,lucchese,cherepkov} over the
years. Furthermore, for the purpose of obtaining HHG spectra for
molecules in the intense laser pulses, high-precision cross sections
in a narrow energy region such as those measured with synchrotron
radiation are not needed. Instead, moderately accurate results over
a broad range of energy such as those based on the one-electron
model developed by Tonzani \cite{tonzani} is likely adequate.

The significance of the present result is not limited to a
workable theory for HHG from molecules. Most importantly, the
results suggest that it is possible to extract photoionization
cross sections over a broad energy range from the measured HHG
spectra which may then be used to unravel the structure of the
target molecule. Since infrared lasers with sub-ten femtoseconds
are already widely available, high-order harmonics generated from
molecules undergoing transformation may then be used to extract
the structure of the transient molecules at different time delays.
The potential of using HHG from molecules for time-resolved
chemical imaging appears quite promising.

This work was supported in part by Chemical Sciences, Geosciences
and Biosciences Division, Office of Basic Energy Sciences, Office
of Science, U.S. Department of Energy. M.L. acknowledges support
from Deutsche Forschungsgemeinschaft.

\section*{References}


text:
\section{QUANTUM MECHANICAL APPROACH}

The quantum scattering theory of a linear rigid-rotor with a structureless atom
has been developed  by \citet{art63}.
In the current calculations, we treat CO as a rigid-rotor with  bond-length fixed at
its equilibrium distance.
Calculations were carried out using the quantum
CC method and the coupled-states (CS) approximation
\citep[see, for example,][]{flo07}.
The two-dimensional CO-H interaction potential is given by $V(R,\theta)$,
where $R$ is the separation between
 the CO center of mass and the H atom,
and  $\theta$ is the angle between $\vec{R}$ and the CO molecular axis.
 The PES was expressed as
\begin{equation}
V(R,\theta)=\sum_{\lambda}^{\lambda_{\rm max}}v_{\lambda}(R)P_{\lambda}(\textrm{cos} \theta),
\end{equation}
\label{eq_vr}
where $P_{\lambda}$ are Legendre polynomials of order $\lambda$.

For a rotational transition from an initial state  $j$ to
a final state  $j'$,
the degeneracy-averaged-and-summed integral cross section
can be given, within the CC method, as
\begin{equation}
\sigma_{j\rightarrow j'}(E_{j})
=\frac{\pi}{(2j+1)k_{j}^2}\sum_{J=0}(2J+1)\sum_{l=|J-j|}^{J+j}
\sum_{l'=|J-j'|}^{J+j'}|\delta_{jj'}\delta_{ll'}
-S_{jj'll'}^J(E_j)|^2,
\label{eq_cross}
\end{equation}
where $S_{jj'll'}^J$ is an element of the scattering matrix and $\vec{j}$ and $\vec{l}$ denote the rotational angular momentum of the CO molecule
and the orbital angular momentum of the collision complex, respectively.
$\vec{J}$ is the total angular momentum given by
 $\vec{J}=\vec{l}+\vec{j}$.
 $k_{j}=\sqrt{2\mu(E-\epsilon_{j})}/\hbar$  is the wave vector
for the initial channel,  $E$ the total energy, $\epsilon_j$ the rotational energy of
CO, $E_j$ the relative translational energy for the
initial channel, and $\mu$ the reduced mass of the CO-H system.
In the CS formulation, the integral cross section is
\begin{equation}
\sigma_{j\rightarrow j'}(E_{j})
=\frac{\pi}{(2j+1)k_{j}^2}\sum_{J=0}(2J+1)\sum_{\Omega=0}^{\Omega_{\bf max}}
(2-\delta_{\Omega0})|\delta_{jj'}
-S_{jj'}^{J\Omega}(E_j)|^2,
\label{cs_cross}
\end{equation}
where $\Omega$ is the projection of the angular momentum quantum number of the diatom
along the body-fixed axis.

The nonreactive scattering code MOLSCAT \citep{molscat} was used in all our calculations.
The modified
log-derivative Airy propagator of \citet{ale87}, with a variable step-size, was adopted to solve
the coupled-channel equations.
The propagation was carried out from $R=1$ $a_0$ to a maximum distance of
$R=50$ $a_0$, where $a_0$ is the atomic unit of length (the Bohr radius).
Cross sections were calculated for collision energies between
 $10^{-5}$  and 15,000 cm$^{-1}$.
The CC calculations have been carried out for initial kinetic energies
up to 11,000 cm$^{-1}$ on the MRCI potential.
Computations with the CS approximation were performed for the MRCI, CCSD,
and WKS potentials above 2000 cm$^{-1}$, but CC calculations were performed for the
latter two surfaces for collision energies less than 2000 cm$^{-1}$.
The angular dependence of the
potential expressed in Eq.~(1) was expanded to $\lambda_{\rm max}=20$
with a 22-point Gauss-Legendre quadrature.
In order to ensure convergence of the state-to-state cross sections,
at least five to ten closed channels in the basis and a sufficient
number of angular momentum partial waves were included in our calculations.
The CO rotational constant of 1.9225 cm$^{-1}$ \citep{lov79} was used to calculate the
rotational energy levels.

\section{RESULTS AND DISCUSSION}

\subsection{State-to-state deexcitation cross sections}

Calculations of
state-to-state rotational quenching cross sections have been carried out for initial
levels $j=1, \  2, \ \cdots, \ 5$ for the deexcitation of CO($v=0,j$) by
H.\footnote{All state-to-state
cross sections and rate coefficients
are available on the UGA Molecular Opacity Project website
(www.physast.uga.edu/ugamop/). The rate coefficients are also available in
the format of the Leiden Atomic and Molecular Database (LAMDA) \citep{sch05}
and in BASECOL \citep{dub06} format on our website.}
To illustrate the different behavior of the cross sections on the WKS, CCSD, and
MRCI potentials, as an example we present in Fig.~\ref{fig1}, a comparison of
the $j=1\rightarrow 0$ deexcitation cross section.
Due to the limited $R$ range (3-20 $a_0$) of the CCSD potential,
the quenching cross sections are given for collision energies greater than  0.1 cm$^{-1}$.
From Fig.~\ref{fig1},  one can see that the cross sections on the WKS PES are significantly
larger than that obtained using the CCSD and MRCI potentials,
especially at ultracold collision energies.
In the energy region between 0.1 and $\sim$20 cm$^{-1}$ ,
the cross sections computed on all three potentials exhibit
resonances caused by the van der Waals interaction,
but due to differences in the potential well depths \citep[see][]{she07}
very different structures are displayed, in particular the resonances on the WKS potential
are largely suppressed.
There exists generally good agreement between the cross sections on the CCSD and MRCI PESs
above $\sim$100 cm$^{-1}$, while the cross sections on all three surfaces appear to
converge above 10$^4$ cm$^{-1}$.
From Fig.~\ref{fig1}, it can also be seen that the CS approximation
is in excellent agreement with the CC method.

For initial rotational states $j$=2, 3, 4, and 5, the
state-to-state quenching cross sections were calculated
using the MRCI PES of \citet{she07} only.
Figure~\ref{fig2} presents the cross sections
from initial $j$=2, 3, and 4 into the individual final states $j'$.
Over the entire collision energy range considered,
similar structure is observed for the different initial rotational $j$ states.
At collision energies below $\sim$10$^{-3}$ cm$^{-1}$, only $s$-wave ($l=0$)
scattering in the incident channel contributes, with the cross section for
inelastic collisions varying inversely with the relative velocity,
indicative of Wigner's threshold Law  \citep{wig48}.
Similar to that shown in Fig.~\ref{fig1}, a number of resonances
at collision energies  between $\sim$1  and $\sim$30 cm$^{-1}$ are evident,
but the resonances become suppressed compared with the $j=1\rightarrow$~0 resonances
\citep[see also][]{yan06}.
In general, the cross sections increase
with increasing $j'$, with that for $j'=0$ being the smallest for initial levels $j$=3 and 4.
For each initial $j$ state, the  $|\Delta j|=|j'-j| = 2$  transition is seen to dominate the
quenching.
This finding is consistent with the fact that odd-$\Delta j$ transitions in homonuclear
molecules are forbidden and since CO is nearly homonuclear, it should follow an even-$\Delta j$ propensity \citep[e.g.,][]{lee87,she07}.

\subsection{State-to-state deexcitation rate coefficients}

By thermally averaging state-to-state quenching cross sections
over a Maxwellian distribution of collision energy, the state-to-state quenching rate
coefficients for deexcitation from specific initial rotational states  can be obtained
at a temperature $T$.  For the $j=1\rightarrow j'=0$ rate coefficients
in Fig.~\ref{fig3}, we compare the present
results on the CCSD and MRCI potentials
with previous calculations of \citet{chu75}, \citet{gre76},  and BYD.
Significant differences between the present results and those
of \citet{chu75} and BYD are evident.
In particular, the BYD rate coefficients obtained using the WKS potential
are over an order of magnitude larger than the present results.
However, the results of \citet{gre76} are in agreement with the present
rate coefficients on the MRCI PES, though their adopted semiempirical potential
is expected to be less accurate.
As one can see from Fig.~\ref{fig3}, the current rate coefficients on the
CCSD and MRCI potentials agree very well for temperatures above 150 K, though
differences occur for lower temperatures. The cross section resonances, as shown
in Fig.~\ref{fig1}, are responsible for the undulations in the rate coefficients.

For $j$=2, 3, 4, and 5,
the state-to-state deexcitation rate coefficients were computed
using the MRCI potential only
and are presented in Figures~\ref{fig4} and \ref{fig5}
for temperatures ranging from 1 to 3000 K.
For temperatures below  $\sim$20~K, the
rate coefficients display a less prominent temperature dependence
due to suppression of the  resonances in the cross sections.
For temperatures above $\sim$20~K, the rate coefficients show a general
increasing trend for all transitions.
Further, the rate coefficients increase with decreasing
$\Delta j$, except for $\Delta j=-2$ transitions which dominate.
For those dominant transitions out of each initial level,
$2 \rightarrow 0$, $3 \rightarrow 1$, $4 \rightarrow 2$, and $5 \rightarrow 3$,
 the rate coefficients
of \citet{chu75}, \citet{gre76}, and BYD are also
displayed in Figures~\ref{fig4} and \ref{fig5}.
Similar to Fig.~\ref{fig3} for $j=1\rightarrow 0$,
Figures~\ref{fig4} and \ref{fig5} show large differences between
the present rate coefficients and
those of BYD as well as those of \citet{gre76}. However,
the rate coefficients of \citet{chu75} approach the current results.

\section{ASTROPHYSICAL APPLICATIONS}

Most of the observed FIR/submm absorption and emission lines originate from rotational
transitions due to simple molecular species such as CO.
To accurately model
molecular spectral features, in addition to spectroscopic data, a knowledge of
collisional excitation rate coefficients due to the dominant impactors is required.
As discussed in the Introduction,
carbon monoxide is nearly ubiquitous, being detected in a vast number of objects with
varying physical, chemical, and radiative properties. As a consequence, CO plays an
important role in the thermal balance and chemistry and is a key target of many modeling efforts.
For example, the radiative and collisional rates of CO were used by Martin et al. (2004)
in estimating the gas density and kinetic temperature in a model of the Galactic Center molecular
region. In a study of the observed properties of Galactic diffuse gas,
Pety et al. (2008) used  the $j$=1~$\rightarrow$~0 and 2~$\rightarrow$~1 rotational
transitions of CO to perform non-LTE models of CO brightness, while
Srianand et al. (2008) detected CO absorption
in a damped Lyman-$\alpha$ system at high redshift.
From their modeling with the radiative transfer code RADEX \citep{tak07}, \citet{sri08} found
that the cosmic microwave background dominates the CO excitation and that
the rotational excitation temperatures of CO are higher than measured in the
Galactic ISM. More highly
excited CO transitions, from $j=14$ to 30, have been observed with {\it Herschel/PACS} in
the extragalactic source NGC 1068 \citep{hai12}.

To investigate the excitation mechanism of CO rotational and rovibrational emission seen in Herbig Ae
disks \citep[for example, the {\it Herschel/PACS} observations of][]{mee12},
Thi et al. (2013) have constructed a complete CO rovibrational model
which adopted existing rate coefficients of CO in collision with H, He, H$_2$, and electrons into
the radiative-photochemical code ProDiMo. Due to the incompleteness of  the
available rate coefficients,
the missing values were estimated using
various scaling laws. For the case of CO-H, the rate coefficients of BYD were used, though
they are not reliable as pointed by
Shepler et al. (2007) and indicated in Figures \ref{fig3}-\ref{fig5}. From  modeling of protoplanetary disks, \citet{thi12} found that the
abundance of atomic H
in the CO line-emitting region is significant so that CO-H collisions are important.

As CO collisional rate coefficients are needed for CO line intensity and
abundance profile modeling, the  current  CO-H rate coefficients will be very
useful for interpreting
observed CO emission lines. In particular, recent observations with
{\it Herschel/PACS} and {\it Herschel/HIFI} provide a large amount of CO FIR/submm data
for  warm gas around evolved stars.
Modeling the CO emission with accurate collisional rate coefficients will enhance
our ability to understand the
evolution of these objects. The availability of the new rotational CO-H rate coefficients provided
here, along with our earlier CO-H$_2$ \citep{yan10} results, should improve CO modeling for the
environments discussed above and many others.
Further, the excitation rates due to H and H$_2$ become comparable for abundance
 ratios $n_{\rm H} /n_{\rm H_2}$ of 2-3 at 10~K and 0.5-1 at 3000~K, highlighting the importance
 of H collisions in  regions with a significant atomic hydrogen fraction.

\section{CONCLUSIONS}

We have performed explicit quantum-mechanical calculations of
rotational quenching of CO due to H collisions
using the close-coupling approach
and the coupled-states approximation on two different interaction potential surfaces.
Deexcitation cross sections reveal even-$\Delta j$ propensity for the CO-H system
 over the collision energy range considered as expected on physical grounds.
State-to-state quenching rate coefficients
for initial CO rotational levels $j=$1, 2, $\cdots$, 5 to all lower $j'$ levels
were computed for temperatures ranging from 1 to 3000 K and are available in tables formatted
for astrophysical applications.
Discrepancies in the  cross sections and rate coefficients on the two new potentials
are due to the differences in the well depth and anisotropy of the surfaces.
Comparison of the present rate coefficients with previous theoretical results
show significant differences primarily due to the quality of the adopted potentials.
The present CO-H rate coefficients can be used to aid astrophysical
modeling, while quenching from higher rotational excited states as well as vibrational
transitions will be studied in the future.

\begin{acknowledgements}

 BHY, PCS, and JMB acknowledge support from NASA under Grant No.
 NNX12AF42G, while the work of NB and RCF was partially supported
by the NSF through grant PHY-1205838 and grant PHY-1203228,
respectively.

\end{acknowledgements}

redpajama_set_name: RedPajamaArXiv