\documentclass[12pt]{iopart}
\usepackage{iopams}
\usepackage{graphicx}% Include figure files
\renewcommand{\vec}[1]{\bi{#1}}% Vectors
% \bfeta=\vec{\eta} !
\renewcommand{\eref}[1]{equation~(\ref{#1})}% For some reason "equation" was not printed before
\newcommand{\eeref}[2]{equations~(\ref{#1}) and (\ref{#2})}% equations (M) and (N)
\newcommand{\frmijk}{v}% former \mathrm{f}
\newcommand{\fijk}{h}% former f
\newcommand{\betaijk}{w}% former \beta
\newcommand{\ps}{phase space}
\newcommand{\ha}{harmonic approximation}
\newcommand{\ho}{harmonic oscillator}
\newcommand{\wf}{Wigner function}
\newcommand{\jp}{jumping point}
\newcommand{\bos}{Born - Oppenheimer surface}
\newcommand{\dr}{Duschinsky rotation}
\newcommand{\se}{Schr\"odinger equation}
\begin{document}
\title[Most probable path in \ps]
{Most probable path in \ps\ for\\
a radiationless transition in a polyatomic molecule}
\author{Alexey V. Sergeev and Bilha Segev\footnote[3] {To whom correspondence should be addressed (bsegev@bgumail.bgu.ac.il)}}
\address{Department of Chemistry, Ben-Gurion University of the Negev,\\
POB 653, Beer-Sheva, 84105 Israel}
\begin{abstract}
We study a radiationless transition in a polyatomic molecule when
the electronic energy of an excited electronic state is
transferred to the vibrational degrees of freedom of the nuclei,
and when some nuclear coordinates change abruptly. This jump
between the donor energy surface to the acceptor one gives the
initial conditions for the subsequent dynamics on the acceptor
surface, and the partition of energy between competing accepting
modes. In the Wigner representation, the physical problem of
recognizing the accepting modes for a radiationless vibronic
relaxation reduces to the mathematical problem of finding the
maximum of a function of many variables under a constraint. The
function is the initial \wf\ of the nuclei and the constraint is
energy conservation. In a \ha\ for the potential surfaces, the
problem is equivalent to finding the distance from a given point
to a multidimensional ellipsoid. This geometrical problem is
solved in closed form. For nonharmonic potentials, the
optimization problem is solved perturbatively.
\end{abstract}
\pacs{33.50.Hv, 82.20.Rp, 03.65.Sq, 2.30.Wd}
%% http://www.physica.org/pacs_index.html
%% 33.50.Hv Radiationless transitions, quenching
%% 82.20.Rp Energy distribution and transfer; relaxation (see also 31.70.H Time-dependent phenomena—in atomic and molecular physics)
%% 03.65.Sq Semiclassical theories and applications
%% 02.30.Wd Calculus of variations and optimal control
%% 31.15.-p Calculations and mathematical techniques in atomic and molecular physics (excluding electron correlation calculations) (for computational techniques, see 02.70)
%% 33.70.Ca Oscillator and band strengths, lifetimes, transition moments, and Franck–Condon factors
\submitto{\JPA}
% Comment out if separate title page not required
\maketitle
\section{\label{int}Introduction}
Molecules are made of heavy particles, the nuclei, and light
particles, the electrons. In the Born - Oppenheimer approximation
one uses this separation of scales in mass and hence in
velocities to obtain electronic terms. Each electronic term
corresponds to some wavefunction of the electrons: $S_0$ is the
ground electronic state, $S_1$ and $S_2$ the first and second
singlet excited states, etc. For a given electronic term, the
nuclei are moving on a potential surface created by the
electrostatic repulsion between them and the attractive force of
the chemical bounds that the electrons create. As long as the
molecule stays at the same electronic state, the motion of the
nuclei can be analyzed by propagating wave packets or by a
semi-classical solution for the \se\ or a classical solution of
the equation of motion with a given potential. What happens to the
nuclei during an electronic transition, say a radiative one, in
which a photon is emitted or absorbed, or in a radiationless
(forbidden) transition? Usually, nothing. Since the nuclei are so
heavy and the electronic transition is so fast, the initial
conditions for the motion of the nuclei on the accepting
potential surface are then given by the position and the velocity
of the nuclei on the donor potential surface. This is called a
vertical transition. Sometimes, however a vertical transition can
not take place, for example when it violates energy conservation.
In this case, the nuclei must ``jump" during the electronic
transition. The jump can be a change in position or momentum of
one or many nuclei on the fast time scale of the electronic
transition. Predicting this jump for each molecule and transition
is a difficult problem and an old one in molecular physics and
photochemistry. We have developed a simple method to predict this
jump. In a number of cases our results indicate that the
dimension of the problem reduces considerably: if only one
nucleus moves out of, say 14, in a given molecule, we need not
consider the other 13 nuclei with the same careful detailed
analysis when calculating and studying properties of the
transition. Usually one form of motion of the nuclei, which is
called the accepting mode (and need not be a normal mode) does
all the ``jumping", while the rest of the degrees of freedom
undergo a practically vertical transition. The theory of surface
jumping, as we call this process, was developed in
Ref.~\cite{Segev:2000:PSD}, where it was also applied to a simple
model. Recently, we have applied it to a 30-dimensional model of
the benzene molecule.
The mechanism of surface jumping is important for applications in
molecular physics and photochemistry \cite{Siebrand:1967:RTP,
Berry:1970:DCP}. The subject of this paper however is not the
application of this mechanism to specific molecules but the
general mathematical problem one is facing when analyzing these
jumps. We aim here at a rigorous formulation of the problem and
its general formal solution. The physical system is presented in
\sref{physical}, the mathematical formulation of the problem is
given in \sref{math}. The complete solution in the \ha, examples
of simple cases, and corrections due to anharmonicity are
respectively studied in sections \ref{harm}, \ref{simple}, and
\ref{anharmsec}. \Sref{sum} concludes.
\section{\label{physical}The physical system}
Let us consider a molecule described by a potential
$V(\vec{R},\vec{r})$ where $\vec{R}$ and $\vec{r}$ are
coordinates of the nuclei and electrons respectively. In the
adiabatic approximation, the quantum states of the molecule are
determined in two steps. First, the \se\ is solved with respect to
the coordinates of the electrons for an arbitrary frozen
configuration of the nuclei described by coordinates $\vec{R}$
that are treated as parameters of the \se. The result is an
eigenfunction $\varphi _{n}(\vec{r}\,;\vec{R})$ and an eigenvalue
$U_{n}( \vec{R})$ where $n$ is a set of electronic quantum
numbers. The second step is to solve the \se\ with respect to the
coordinates of the nuclei, $\vec{R}$, in an effective potential
$U_{n}(\vec{R})$ (which was found in the first step), and to find
eigenfunctions $\chi _{N,n}(\vec{R})$ and eigenvalues $E_{N,n}$
where $N$ is a set of vibrational and rotational quantum numbers
of the nuclear motion (in this equation, $n$ are parameters
because the potential $U_{n}(\vec{R})$ depends on $n$). The
result of the adiabatic approximation is the wavefunctions in the
form
\begin{equation}
\Psi _{N,n}(\vec{R},\vec{r})=\chi _{N,n}(\vec{R})\varphi
_{n}(\vec{r}\,;\vec{ R}) \label{adi}
\end{equation}
and the corresponding energies $E_{N,n}$. Even when the adiabatic
approximation does not apply, one can use these wavefunction as a
basis set.
We study the following problem. Let some initial state $\Psi
_{(N_{\mathrm{I}}, n_{\mathrm{I}})}$ of a molecule have quantum
numbers $(N_{ \mathrm{I}}, n_{\mathrm{I}})$. Now, suppose that
there exist several possible final states with another set of
electronic quantum numbers $n_{\mathrm{F}}$ and different sets of
nuclear quantum numbers $N_{\mathrm{F}}$, e.g.
\begin{equation}
(N_{\mathrm{F}}^{(1)}, n_{\mathrm{F}}), (N_{\mathrm{F}}^{(2)},
n_{\mathrm{F} }), (N_{\mathrm{F}}^{(3)}, n_{\mathrm{F}}),...
\label{fin}
\end{equation}
Since the state $\Psi _{(N_{\mathrm{I}},n_{\mathrm{I}})}$ is
actually a mixed quantum state due to small nonadiabatic effects,
in the process of its evolution in time there is nonzero
probability of finding the molecule in one of states (\ref{fin})
even if it was initially prepared in the state $ (N_{\mathrm{I}},
n_{\mathrm{I}})$. According to the theory of radiationless
transitions \cite{Jortner:1969:RTP}, this probability is maximal
for states having the same energy as the initial state
$(N_{\mathrm{I}}, n_{\mathrm{I}})$, and is proportional to the
density of final states multiplied by the square of the Frank -
Condon integral -- an overlap integral between the nuclear
components of the wavefunctions,
\begin{equation}
\int \chi _{N_{\mathrm{I}}, n_{\mathrm{I}}}(\vec{R})\chi
_{N_{\mathrm{F}}, n_{\mathrm{F}}}(\vec{R})d\vec{R}, \label{fci}
\end{equation}
where here and in the following all integrals are from $-\infty$
to $+\infty$.
The purpose of this work is to develop a method of choosing a
state $(N_{ \mathrm{F}}^{\ast }, n_{\mathrm{F}})$ or a
superposition of such states among all possible final states
(\ref{fin}) satisfying the energy conservation condition
\begin{equation}
E\left( N_{\mathrm{F}}^{\ast }, n_{\mathrm{F}}\right) =E\left(
N_{\mathrm{I} }, n_{\mathrm{I}}\right) \label{ecc}
\end{equation}
for which the square of the integral (\ref{fci}) reaches its
maximum. This state is the most preferable accepting mode for a
radiationless transition.
The idea of our approach is to use the Wigner transformations of
the wavefunctions. The Wigner transform of a given wavefunction
$\psi (\vec{R})$ is defined as
\begin{equation}
\rho (\vec{R},\vec{P})=\left( \frac{1}{2\pi }\right)^{N}\int
d\bfeta \rme^{-\rmi\vec{P}\cdot \bfeta}
\psi^\star(\vec{R}+\bfeta/2)\psi (\vec{R}-\bfeta/2), \label{wig}
\end{equation}
where $N$ is the number of independent coordinates. Here and
henceforward we use units where $\hbar=1$. In the Wigner
representation an overlap integral squared can be rewritten as an
integral over \ps\ coordinates,
\begin{equation}
\left| \int\psi _{1}^{\ast }\psi _{2}d \vec{R}\right| ^{2}=\left(
{2\pi }\right)^{N} \int \int d\vec{R} d\vec{P}\rho _{1}\rho _{2},
\end{equation}
where $\rho _{1}=\rho _{1}(\vec{R},\vec{P})$ and
$\rho_{2}=\rho_{2}(\vec{R},\vec{P})$ are Wigner transforms of the
functions $\psi _{1}=\psi _{1}(\vec{R})$ and $\psi _{2}=\psi
_{2}(\vec{R })$, respectively.
%If the function $\psi_{{\rm f}}(x)$ is a wavefunction of a Hamiltonian $%
%H(x,p)$ with energy $E$, then up to physical intuition the \ps\
%distribution $\rho_{{\rm f}}(x,p)$ peaks at points where $H(x,p)=E$.
The total rate of transition from a state $(N_{\mathrm{I}},
n_{\mathrm{I}})$ to a manifold of states
$(N_{\mathrm{F}},n_{\mathrm{F}})$ with a definite
$n_{\mathrm{F}}$ and all possible $N_{\mathrm{F}}$ is proportional
to a sum
\begin{equation}
\sum_{E\left( N_{\mathrm{F}},n_{\mathrm{F}}\right) =E}\left( \int
\chi _{N_{\mathrm{I}},n_{\mathrm{I}}}( \vec{R})\chi
_{N_{\mathrm{F}},n_{\mathrm{F}}}(\vec{R})d\vec{R}\right) ^{2}.
\label{rat}
\end{equation}
where both the Frank-Condon factor and the density of final
states are included in the expression, $E=E\left( N_{\mathrm{I}
},n_{\mathrm{I}}\right)$, and $\chi_{N_{\mathrm{I}},
n_{\mathrm{I}}}(\vec{R})$ is assumed to be real. In terms of \wf
s (\ref{rat}) is proportional to
\begin{equation}
{\left( 2\pi \right) ^{N}}\int \int d\vec{R}d\vec{P}\rho _{N_{
\mathrm{I}},n_{\mathrm{I}}}\sum_{E\left( N_{\mathrm{F}},n_{
\mathrm{F}}\right) =E}\rho _{N_{\mathrm{F}}, n_{\mathrm{F}}},
\label{wigrat}
\end{equation}
where $\rho_{N,n}=\rho_{N,n}(\vec{R},\vec{P})$. Here we study the
expression in (\ref{wigrat}) to be integrated. We are especially
interested in finding a maximum of this integrand. Importance of
the point of maximum of the \ps\ integrand was stressed in the
paper \cite{Segev:2000:PSD} where the \ps\ derivation of
propensity rules for energy transfer processes between \bos s was
presented.
We use the following approximation for the second factor in the
integrand of (\ref{wigrat}):
\begin{equation}
(2\pi)^N \sum_{E\left( N_{\mathrm{F}}, n_{\mathrm{F}}\right)
=E}\rho _{N_{\mathrm{F} },n_{\mathrm{F}}}(\vec{R},\vec{P})=
\delta \left( H_{\mathrm{F}}(\vec{R},\vec{ P})-E\right)
\label{cla}
\end{equation}
which is equivalent to replacing this function by the zero-order
classical term of its semiclassical expansion in powers of $\hbar
^{2}$ \cite {Hupper:1998:USE,Heller:1978:QCC}. The approximation
(\ref{cla}) reduces the integral (\ref {wigrat}) to
\begin{equation}
\int_{H_{\mathrm{F}}(\vec{R},\vec{P} )=E}
|\vec{\nabla}H_{\mathrm{F}}|^{-1} \rho _{N_{\mathrm{I}},
n_{\mathrm{I}}}(\vec{R},\vec{P})d\vec{R}d\vec{P}, \label{rat1}
\end{equation}
where $|\vec{\nabla}H_{\mathrm{F}}|=
\left[(\vec{\nabla_{\vec{R}}}H_{\mathrm{F}})^2+
(\vec{\nabla_{\vec{P}}}H_{\mathrm{F}})^2\right]^{1/2}$ is the
gradient of the function $H_{\mathrm{F}}$ in \ps. The rest of the
paper is devoted to finding a maximum of the \wf\ $\rho
_{N_{\mathrm{I}}, n_{\mathrm{I}}}(\vec{R},\vec{P})$ on an
equipotential surface defined through the equation
$H_{\mathrm{F}}(\vec{R},\vec{P})=E$ for both harmonic and
anharmonic potentials. In doing so we set the ground for the
future analysis of radiationless transitions of specific large
polyatomic molecules. In addition, we formulate and prove some
general yet simple thumb rules for predicting the accepting mode
of a given radiationless transition.
\section{\label{math}Formulation of the problem}
The Hamiltonian of the acceptor is approximated by a \ho\ plus
third order anharmonic terms,
\begin{equation}
H_{\mathrm{F}}=\frac12\sum_{i=1}^{N}\left( p_{i}^{2}+\omega
_{i}^{2}q_{i}^{2}\right) +\frac{1}{6}
\sum_{i,j,k=1}^{N}\frmijk_{ijk}q_{i}q_{j}q_{k}, \label{finham}
\end{equation}
where $p_{i}$ and $q_{i}$ are mass weighted normal momenta and
coordinates, $ q_{i}=R_{i}\sqrt{m_{i}}$ and
$p_{i}=P_{i}/\sqrt{m_{i}}$. Similarly, the Hamiltonian of the
donor surface is
\begin{equation}
H_{\mathrm{I}}= \frac12\sum_{i=1}^{N}\left( {p_{i}^{\prime
}}^{2}+{\omega _{i}^{\prime }}^{2}{q_{i}^{\prime }}^{2}\right)
+\frac{1}{6} \sum_{i,j,k=1}^{N}\frmijk_{ijk}^{\prime
}q_{i}^{\prime }q_{j}^{\prime }q_{k}^{\prime }. \label{iniham}
\end{equation}
The mass weighted normal coordinates $p_{i}^{\prime
}=P_{i}^{\prime }/\sqrt{ m_{i}^{\prime }}$ and $q_{i}^{\prime
}=R_{i}^{\prime }\sqrt{m_{i}^{\prime }}$ are generally some
linear combinations of $p_{i}$ and $q_{i}$,
\begin{eqnarray}
\vec{q}^{\prime } &=&\mathbf{S}\left(
\vec{q}-\vec{q}^{(0)}\right), \nonumber
\\
\vec{p}^{\prime } &=&\mathbf{S}\vec{p} \label{doncor}
\end{eqnarray}
where $\mathbf{S}$ is an orthogonal $N\times N$ matrix
($\mathbf{S}^{\mathrm{ T}}=\mathbf{S}^{-1}$) and the vector
$\vec{q}^{(0)}$ corresponds to the change of the equilibrium
structure of the molecule relative to the donor state. An element
$S_{ij}\neq \delta _{ij}$ only when the $i$th and $j$th normal
coordinates have the same symmetry (so called \dr). The same
matrix transforms both $\vec{q}$ and $\vec{p}$ since the
transformation preserves the commutation relations $
[q_{i}^{\prime },p_{j}^{\prime }]=[q_{i},p_{j}]=i\hbar \delta
_{ij}$ and since the Hamiltonians (\ref{finham}) and
(\ref{iniham}) have the same kinetic energy term,
$\sum_{i=1}^{N}p_{i}^{2}=\sum_{i=1}^{N}{ p_{i}^{\prime }}^{2}$.
%
We restrict ourselves to the ground state in the donor potential,
\begin{equation}
\chi _{N_{\mathrm{I}},n_{\mathrm{I}}}(\vec{q})=C\exp \left(
-\frac12 \sum_{i=1}^{N}\omega _{i}^{\prime }{q_{i}^{\prime
}}^{2}\right) +\chi _{1},
\end{equation}
where $C$ is a normalization factor, and $\chi _{1}$ is the first
anharmonic correction (a linear function of the coefficients
$\frmijk_{ijk}$) derived in \sref{anharmsec} below. The Wigner
transform of $\chi _{N_{\mathrm{I}},n_{\mathrm{I}}}( \vec{q})$ is
$C^{\prime }\exp \left( -2W\right) $ where $C^{\prime }$ is a
constant pre-factor,
\begin{equation}
W=\frac12\sum_{i=1}^{N}\frac{{1}}{\omega _{i}^{\prime }}\left( {
p_{i}^{\prime }}^{2}+{\omega _{i}^{\prime }}^{2}{q_{i}^{\prime
}}^{2}\right) +W_{1} \label{wig0}
\end{equation}
and $W_{1}$ is the first anharmonic correction derived in
\sref{anharmsec} below.
The jumping between the donor and acceptor states occurs at a
point of minimum of $W$ subject to a constraint
$H_{\mathrm{F}}=E$. There are several approaches to solve a
problem of constraint minimum \cite {Beveridge:1970:OTP}. One
could use a method of direct substitution by eliminating one of
the variables from the function $W$. This method is not
symmetrical with respect to the treatment of the variables
$\{x_{i}\}$ ($i=1,2,\ldots,M$, $M=2N$) that are arguments of the
functions $W$ and $H$, i.e. $\{q_i, p_i\}$ for $i=1, 2, \ldots,
N$ (later these collective variables will be redefined). To avoid
distinction between the variables, we use a method of Lagrange
multiplier by introducing an undetermined constant $\lambda $ and
forming a function $F(\vec{x},\lambda )=W(\vec{x})-\lambda
H(\vec{x})$. This function is to be made stationary with respect
to all variables $\{x_{i}\}$, so that
\begin{equation}
\frac{\partial F}{\partial x_{i}}\left( \vec{x}^{\ast },\lambda
^{\ast }\right) =0 \label{opt}
\end{equation}
for $i=1,2,...,M$, and the constant $\lambda ^{\ast }$ is to be
selected so that
\begin{equation}
H\left( \vec{x}^{\ast }\right) =E \label{opte}
\end{equation}
Conditions (\ref{opt}) and (\ref{opte}) provide a system of $M+1$
equations for $M+1$ unknowns, $x_{1}^{\ast },x_{2}^{\ast
},...,x_{M}^{\ast }$, and $ \lambda ^{\ast }$ which can be
briefly summarized as:
\begin{equation}
\vec{\nabla}W=\lambda \vec{\nabla}H,\quad H=E \label{lagr}
\end{equation}
The Lagrange multiplier $\lambda $ has concrete physical meaning.
Since
\begin{equation}\fl
\frac{\rmd}{\rmd E}W\left( \vec{x}^{\ast }\right)
=\sum_{i=1}^{M}{\frac{\partial W}{
\partial x_{i}}\left( \vec{x}^{\ast }\right) \frac{dx_{i}^{\ast
}}{{dE}}} =\lambda ^{\ast }\sum_{i=1}^{M}{\frac{\partial
H}{\partial x_{i}}\left( \vec{ x}^{\ast }\right)
\frac{dx_{i}^{\ast }}{{dE}}}=\lambda ^{\ast }\frac{d}{dE} H\left(
\vec{x}^{\ast }\right) =\lambda ^{\ast },
\end{equation}
the parameter $\lambda ^{\ast }$ is the sensitivity of the
minimum value of $W$ to the energy gap.
After finding all the stationary points $\vec{x}^{\ast }$ it is
necessary to determine for each point whether it is a minimum of
the function $W$ under restriction (\ref{opte}), a saddle point or
a maximum, and which of all the local minima gives the smallest
value for $W$. The global minimum found in this way is a true
solution of the optimization problem, see \fref{contours}.
\begin{figure}[bp]
\begin{center}
\includegraphics[scale=0.7]{contours.eps}
\end{center} \caption{\label{contours}Finding the minimum of the
function $W$ under the energy constraint $H=E $. For this example,
$W(x_1,x_2)= 0.4(x_1-0.1)^2+ 0.6(x_2-0.2)^2+ 0.05[(x_1-0.1)^3+
(x_1-0.1)^2 (x_2-0.2)+ (x_1-0.1) (x_2-0.2)^2- (x_2-0.2)^3]$,
$H(x_{1},x_{2})=\frac12x_1^2+\frac12 x_2^2+ 0.1[-x_1^3-3 x_1^2
x_2+2 x_1 x_2^2]$, and $E=1$. Dashed lines represent stationary
points of the function $F=W-\protect\lambda H$, \eref{opt}.
Energy-constraint points satisfying \eref{opte} lie on the border
of the dark area, $H
\alpha_{\mathrm{min}}$ for $j>L$, where $\alpha_{\mathrm{min}}$
is the minimal number among $\alpha_i$ $(i=1,2...,M)$, and $L$ is
the number of entries of $\alpha_{\mathrm{min}}$ in the set
$\{\alpha_i\}$ $(i=1,2...,M)$. When $\lambda$ increases from
$-\infty $ to $\alpha _{\mathrm{min}}$, the function $h(\lambda
)$ monotonously increases from $0$ to
\begin{equation}
E_1= \frac12\sum_{i>L}\left( \frac{\alpha_i}{\alpha _i-\alpha
_{\mathrm{min}}}\right) ^{2}X_i^2
\end{equation}
if $X_1= X_2= \ldots= X_L= 0$, otherwise it increases from $0$ to
$E_1= \infty$.
There are two possible cases. In the first case, when $X_1= X_2=
\ldots= X_L= 0$ and $E> E_{1}$, the minimal root of (\ref{opte})
is $\lambda ^{\ast }= \alpha_{\mathrm{min}}$, $x_{j}^{\ast }$ for
$j>L$ are expressed through $\lambda ^{\ast }$ by (\ref{coord}),
and from (\ref{restr}) we get a set of possible $(x_1^\ast,
x_2^\ast, \ldots, x_L^\ast)$. In the second case, when $E\leq
E_{1}$, there is a unique root $ \lambda ^{\ast }$ of
\eref{restr1} on the interval $(-\infty, \alpha_{\mathrm{min}})$,
the coordinates of this minimum are expressed through $\lambda
^{\ast }$ by (\ref{coord}), and the minimum of $W$ is given by
(\ref{wmin}).
%It can be shown that there are no other cases...
\subsection{Results}
Let us summarize the solution in the \ha. Given an initial \wf\
and an accepting Hamiltonian, applying a \ha\ and a change of
variables, re-enumerating the variables so that $\alpha_1=
\alpha_2= \ldots= \alpha_L $ equal to the smallest of all
$\alpha_i$, and explicitly solving \eeref{wig3}{restr}, we get
the \jp\ for the radiationless transition. There are two cases:
\underline{\bf Case I}
In this case (when $X_1= X_2= \ldots= X_L= 0$ and $E>E_{1}$) there
exist several points of minimum with the same value of $W$. If
$L=1$, then two \jp s differ in sign of the first coordinate
\begin{equation}
x_1^{\ast }= \pm \left[ 2\left( E-E_{1}\right) \right] ^{1/2}.
\label{x1c}
\end{equation}
For $i\neq 1$ coordinates are
\begin{equation}
x_i^{\ast }= \frac{\alpha _{i}}{\alpha _{i}-\alpha
_{\mathrm{min}}}X_{i}. \label{xic}
\end{equation}
If $L>1$, then $(x_1^{\ast }, x_2^{\ast }, \ldots, x_L^{\ast })$
fill a $L-1$-dimensional sphere of radius $\left[ 2\left(
E-E_{1}\right) \right] ^{1/2}$,
\begin{equation}
\frac12(x_1^2+ x_2^2+ \ldots+ x_L^2)= E-E_1,
\end{equation}
which follows from \eref{restr}. The rest of coordinates, for
$i>L$, are defined unambiguously by \eref{xic}.
The minimum of $W$ is
\begin{equation}
W^{\ast }=\alpha _{\mathrm{min}}E-\frac{\alpha
_{\mathrm{min}}}{2}\sum_{i>L}\frac{\alpha _{i}X_{i}^{2}}{\alpha
_{i}-\alpha _{\mathrm{min}}}. \label{wmin2}
\end{equation}
Below, for example in \eref{xine1}, we assume that $L=1$, or that
$\alpha_i> \alpha_1$ for $i=2,3, \ldots,M$.
\underline{\bf Case II}
This case applies when at least one of $X_i$ with $i\leq L$ is
nonzero, or when $X_1= X_2= \ldots= X_L= 0$ and $E\leq E_{1}$. The
coordinates at the \jp\ are given by \eref{coord}:
\begin{equation}
x_{i}^{\ast }=\frac{\alpha _{i}}{\alpha _{i}-\lambda _{\ast
}}X_{i}, \nonumber
%\label{x0b}
\end{equation}
where $\lambda _{\ast }\leq \alpha_{\mathrm{min}}$ is the unique
root of the equation:
\begin{equation}
\frac12 \sum_{i}\left( \frac{\alpha _{i}}{\alpha _{i}-\lambda
}\right) ^{2}X_{i}^{2}=E. \label{minroot}
\end{equation}
\subsection{Discussion}
In a radiative vertical transition, only displaced modes are
excited. The initial conditions for dynamics on the accepting
potential energy surface, which we call the \jp\ are then given
by:
\begin{equation}
x_{i}^{\ast }=X_{i}, \label{vertical}
\end{equation}
The energy that goes into vibration and the value of the
logarithm of the \wf\ at the \jp\ are then given respectively by:
\begin{equation}
E_{0}=\frac12\sum_{i=1}^{M} \alpha _{i}^{2}X_{i}^{2}, \quad
W_{0}=\frac12\sum_{i=1}^{M} \alpha _{i}X_{i}^{2}.
\end{equation}
Energy is conserved because the photon takes the rest of the
energy
\begin{equation}
E_{\rm photon}=E-E_0, \label{photon}
\end{equation}
where E is the energy gap between minima of the donor and
accepting surfaces.
In a radiationless transition there is no photon. The released
electronic energy must become vibrational energy. The two cases I
and II differ in how this energy is distributed between the
different vibrations.
In case I, despite the fact that $X_1=0$, i.e. there is no
displacement along the $x_1$ direction in \ps\ (be it a coordinate
or a momentum), $x_1$ is an accepting mode for this transition
{\it because the initial \ps\ quasidistribution comes closest to
the final energy hypersurface in this direction}. We shall refer
to $x_1$ as the major accepting mode.
%This major accepting mode replaces the photon as the acceptor of
%the access energy relative to vertical transition.
This mode that absorbs the excess energy $E- E_1$ plays the same
role as electro-magnetic field modes during a radiative
transition, when an emitting photon accepts the excess energy
relative to the vertical transition.
In contrast, case II does not look essentially different from a
vertical transition. Namely, only displaced modes are involved in
the transition, and the jump for each coordinate involved is
proportional to this mode's displacement. No momentum jumps
exist, since the Hamiltonians are never displaces along a
momentum \ps\ coordinate. The smaller is $|\lambda _{\ast }|$, the
closest is the transition to a vertical one, as ${\alpha
_{i}}/({\alpha _{i}-\lambda _{\ast })}$ is closer to 1.
Physically a small $|\lambda _{\ast }|$ corresponds to the
special case when $H_F(\vec{X})=E$.
\subsection{Dependence on the energy gap}
In case I the dependence on the energy gap is trivial. All the
\ps\ coordinates at the \jp\ but one do not depend on the energy
gap but only on their respective displacements and the relative
difference of the parameters $\alpha_i$
%widths
between each one of them and the major accepting mode. $x_1^\ast$,
the jumping coordinate of the major accepting mode grows with the
energy gap. The larger is the energy gap - the more important is
this accepting mode.
In order to consider the properties of the \jp\ as a function of
the energy gap $E$ in case II, for both limits of small and large
$E$ \eref{minroot} is rewritten here in a simplified form
\begin{equation}
\sum_{i=1}^{M}\left( \frac{\beta _{i}}{\alpha _{i}-\lambda
}\right) ^{2}\ =1, \label{restr1a}
\end{equation}
where $\beta _{i}=\alpha _{i}\left| X_{i}\right| \left( 2E\right)
^{-1/2}$.
For a small energy gap, \eref{restr1a} has two real roots
$\lambda =\pm \left( \sum \beta _{i}^{2}\right) ^{1/2}$. The
solution corresponding to the minimal (negative) root
asymptotically behaves as
\begin{eqnarray}
\lambda ^{\ast } &=-\left( E_{0}/E\right) ^{1/2}, \label{small1} \\
x_{i}^{\ast } &=\alpha _{i}X_{i}\left( E/E_{0}\right) ^{1/2},
\label{small2}
\\
W^{\ast } &=W_{0}-2\left( EE_{0}\right) ^{1/2}. \label{small3}
\end{eqnarray}
In the limit of a large energy gap, when $\beta _{i}\rightarrow
0$, \eref{restr1a} has $2M$ roots $\lambda =\alpha _{i}\pm \beta
_{i}$, ($i=1,2,...,M$). The minimal root corresponding to the
minimum of $W$ is $\lambda ^{\ast }=\alpha _{1}-\beta _{1}$. It
is clear that $X_{1}=0$ belongs to the case I for sufficiently
large $E$. A perturbative solution, with $\varepsilon
=\mathrm{sign}X_{1}(2E)^{-1/2}$ as a small parameter shows that
for $ X_{1}\neq 0$ as well, although $\lambda ^{\ast }$ depends
on the energy, this dependence approaches zero for a large enough
energy gap, for which:
\begin{eqnarray}
x_{1}^{\ast } &\approx \sqrt{2E}\ \mathrm{sign}(X_{1}), \\
x_{i\neq1}^{\ast } &= \frac{\alpha _{i}}{\alpha _{i}-\alpha
_{1}}X_{i}
%\nonumber \\
-\frac{\alpha _{i}\alpha _{1} }{\left( \alpha _{i}-\alpha
_{1}\right) ^{2}} \frac{|X_{1}|X_{i}}{\sqrt{2E}} +O(\varepsilon
^{2}) \label{xine1}.
\end{eqnarray}
Here, again, in the limit of a large energy gap $x_1$ is the
major accepting mode regardless of its being displaced or not. If
the displacement $X_1=0$ there are two \jp s with opposite signs,
while if $X_1\neq 0$ the sign of the jump is determined by the
sign of the displacement.
\section{\label{simple}Simple cases}
In the previous section a complete solution in the \ha\ was
derived. For any given radiationless transition the accepting
mode(s) can be found by applying this procedure. To gain some
intuition, we apply it here to some simple examples. We
separately check the influence of the frequencies, \dr s, and the
displacements, on the results. We also check the predictive power
of the results on the energy distribution between the modes, for
an example where this energy distribution is well defined.
\subsection{Frequencies}
In the simplest case, when the normal coordinates of the initial
and the final states are the same ($\vec{q}=\vec{q}^{\prime}$,
$\vec{p}=\vec{p}^{\prime}$) and $m_{i}=m^{\prime}_{i}$, the set of
variables $ \{x_{1},x_{2},...,x_{M}\}$ consists of
$\sqrt{m_{i}}\omega _{i}R_{i}$ and $1/\sqrt{m_{i}}P_{i}$, and the
set $\{\alpha _{1},\alpha _{2},...,\alpha _{M}\}$ consists of
$\omega^{\prime}_{i}/\omega _{i}^{2}$ and $1/\omega^{\prime}_{i}$
sorted in ascending order. Since $E_{1}=0$, this system belongs
to the case I considered in the previous section, with
$x_{1}^{\ast }=\pm \sqrt{2E}$ and $ x_{2}^{\ast }=x_{3}^{\ast
}=...=x_{M}^{\ast }=0$. There is only one accepting mode. In
terms of normal mode coordinates it means that if the minimum
number in the set $ \{\omega^{\prime}_{i}/\omega
_{i}^{2},1/\omega^{\prime}_{i}\}$ is $\alpha
_{1}=\omega^{\prime}_{i_{0}}/\omega _{i_{0}}^{2}$, then the
launching point for the transition is at $R_{i_{0}}=\pm \left(
2E/m_{i_{0}}\right) ^{1/2}/\omega _{i_{0}}$ and the other
coordinates and momenta are zero. If $\alpha
_{1}=1/\omega^{\prime}_{i_{0}} $, then $ P_{i_{0}}=\pm \left(
2Em_{i_{0}}\right) ^{1/2}$ and the other coordinates and momenta
are zero.
\subsection{\dr}
\label{dusch}
Now, suppose that $\vec{q}_{0}=0$ and $m_{i}=m^{\prime}_{i}=m$,
$N=2$, and the matrix $\mathbf{S}$ is a unitary matrix of a
general form \[\mathbf{S}= \left(
\begin{array}{cc}
\cos \varphi & \sin \varphi \\
-\sin \varphi & \cos \varphi
\end{array}
\right), \] where $\varphi $ is a rotation angle. Then, $\{\alpha
_{1},\alpha _{2},\alpha _{3},\alpha _{4}\}$ consists of the
following four numbers: $ \{A\pm \left[ A^{2}-a_{1}a_{2}\right]
^{1/2},1/\omega^{\prime}_{1},1/ \omega^{\prime}_{2}\}$, where
$A=\frac12(a_{1}+a_{2})\cos ^{2}\varphi +
\frac12(b_{1}+b_{2})\sin ^{2}\varphi $,
$a_{1}=\omega^{\prime}_{1}/ \omega _{1}^{2}$,
$a_{2}=\omega^{\prime}_{2}/\omega _{2}^{2} $,
$b_{1}=\omega^{\prime}_{2}/\omega _{1}^{2}$,
$b_{2}=\omega^{\prime}_{1}/ \omega _{2}^{2}$. If $\varphi =0$,
then it is $\{a_{1},a_{2},1/\omega^{
\prime}_{1},1/\omega^{\prime}_{2}\}$. If $\varphi =\pi /2$, then
it is $
\{b_{1},b_{2},1/\omega^{\prime}_{1},1/\omega^{\prime}_{2}\}$.
The next example is numerical. It demonstrates that a \dr\ can
influence the \jp\ and the value of the \wf\ at that point.
Suppose that $\vec{q}_{0}=0$ and $m_{i}= m^{\prime}_{i}=m$, $N=3$,
\[\mathbf{S}=\left(
\begin{array}{ccc}
\cos \varphi & \sin \varphi & 0 \\
-\sin \varphi & \cos \varphi & 0 \\
0 & 0 & 1
\end{array}
\right),\] $\omega _{1}=0.6$, $\omega _{2}=0.3$, $\omega
_{3}=0.603$, $ \omega^{\prime}_{1}=0.595$, $\omega^{\prime}_{2}=
0.298$, $\omega^{\prime}_{3}= 0.6$ ($\omega^{\prime}_{i}$ are
taken slightly smaller than $\omega _{i}$ as it usually happens in
molecules). The dependence of $ \{\alpha _{i}\}$ on $\varphi $ is
shown in \fref{exmp}.
\begin{figure}[bp]
\begin{center}
\includegraphics[scale=0.6]{alpha.eps}
\end{center} \caption{\label{exmp}The dependence of the smallest four
eigenvalues $\protect\alpha_{i}$, $ i=1, 2, 3, 4$, on the
rotation angle $\protect\varphi$ for the example of
\protect\sref{dusch}. The rest of the eigenvalues
$\{1/\protect\omega^{\prime}_{2}, A+ \left[ A^{2}-a_{1}a_{2}
\right] ^{1/2}\}$ are larger than $3$. The smallest eigenvalue
$\protect \alpha_\mathrm{min}$ determines the minimum of $W$
($W_{\mathrm{min}}=\protect\alpha_\mathrm{min} E $). }
\end{figure}
The \jp\ for $\varphi >0.033$ differs from the \jp\ at $\varphi
=0$. Here, since $\alpha _{i}$ are independent of the energy, and
because of the fact that
$dW_{\mathrm{min}}/dE=\alpha_\mathrm{min}$, $W_{
\mathrm{min}}=\alpha_\mathrm{min} E$, and the value of the \wf\ is
proportional to $\exp (-2\alpha_\mathrm{min} E)$.
\subsection{Displacements}
Consider a simple case with non-zero displacement, with $N=2$,
$m_{i}=m^{\prime}_{i}$, and $\vec{q}_{0}=(Q,0)$, i.e.\
$q^{\prime}_{1}=q_{1}-Q$, $q^{ \prime}_{2}=q_{2} $. We define:
\begin{eqnarray}
\lambda _{\ast }&=& \frac{\omega^{\prime}_{1}}{\omega
_{1}^{2}}\left[ 1-\left( \frac{m_{1}\omega
_{1}^{2}Q^{2}}{2E}\right) ^{1/2}\right],
\label{lambdaq1} \\
\sigma&=&\min \left(\lambda _{\ast
},\frac{\omega^{\prime}_{1}}{\omega_{1}^{2}},
\frac{\omega^{\prime}_{2}}{\omega_{2}^{2}},
\frac{1}{\omega^{\prime}_{1}},
\frac{1}{\omega^{\prime}_{2}}\right).
\end{eqnarray}
There are four sub-cases. (1) If $\sigma=1/\omega^{\prime}_{1}$ or
$\sigma=\lambda _{\ast }$, there is only one non-zero jumping
coordinate
\begin{equation}
q_1^\ast =Q/(1-\lambda _{\ast }\omega _1^2/\omega^\prime_1).
\label{q1}
\end{equation}
Otherwise, there are two non-zero jumping coordinates
$q_{1}^\ast$ and (2) $q_{2}^{\ast}$ if
$\sigma=\omega^{\prime}_{2}/\omega_{2}^{2}$, (3) $p_{1}^{\ast }$
if $\sigma=1/\omega^{\prime}_{1}$,(4) $p_{2}^{\ast }$ if
$\sigma=1/\omega^{\prime}_{2}$.
The larger is the displacement and the smaller is the energy gap,
the smaller is $\lambda _{\ast }$, and $q_1$ becomes the only
accepting mode. In contrast, in the limit of a very large energy
gap, the displacement no longer plays a role in the minimization
problem predicting the jump. In this limit, the frequencies alone
determine the jump as in the previous case of zero displacement.
\subsection{Predictive power of the \jp}
The latter case with additional simplifications $m_1=m_2=1$ and $
\omega_1=\omega_2=1$ was considered in \cite{Segev:2000:PSD}.
This paper has plots of the initial wave function,
\begin{equation}
\Psi_{\mathrm{I}}(q_1,q_2) = \psi_0(q_1-Q) \psi_0(q_2) \label{wfi}
\end{equation}
vs. the final wave function
\begin{equation}
\Psi_{\mathrm{F}}(q_1,q_2) = \sum_{j=0}^n C_j \psi_j(q_1)
\psi_{n-j}(q_2) \label{wff}
\end{equation}
where $\psi_i(q) $ is a \ho\ wave function, $ E=n+1$, and $C_j$ is
an overlap integral between the ground state $\Psi_{\mathrm{I}}$
and the excited wave function $\psi_j(q_1) \psi_{n-j}(q_2)$.
It was demonstrated that the pattern of the final wave function
depends on the position of the \ps\ jump. Here, we reconsider six
numerical examples from the paper \cite{Segev:2000:PSD} by a
\textit{quantitative} comparison with the \ps\ results. We
calculate partial energies of excitations along two different
modes,
\begin{equation}
E_1=P_E^{-1}\sum^n_{j=0}C_j^2\left(j+\frac12 \right), \quad
E_2=P_E^{-1}\sum^n_{j=0}C_j^2\left(n-j+\frac12\right), \nonumber
\end{equation}
where $P_E=\sum^n_{j=0}C_j^2$ is the total probability of a
transition to $E_F=E$. $E_1$ and $E_2$ are well defined physical
observables because the two dimensional \ho\ here considered is
separable along $q_1$ and $q_2$. They can be calculated exactly
and compared to their \ps\ counterparts
\begin{equation}
E_1^\ast=\frac12\left({p_1^\ast}^2+{q_1^\ast}^2\right) +\frac12,
\quad E_2^\ast= \frac12\left({p_2^\ast}^2+{q_2^\ast}^2\right)
+\frac12, \label{estar}
\end{equation}
where $q_1^\ast$, $p_1^\ast$, $q_2^\ast$, $p_2^\ast$ are the \ps\
coordinates of the jump, \[E= E_1^\ast+ E_2^\ast= n+1.\] Notice
that the classical energies in \eref{estar} are corrected by
incorporating the quantum energy of zero vibrations, $1/2$.
Without the zero vibrational energy, we found that results are
less satisfactory.
We compare the percentage of energy going into the first mode,
exact vs \ps\ result,
\begin{equation}
R_1=E_1/E, \quad R_1^\ast=E_1^\ast/E \label{perc}
\end{equation}
\Tref{r1} shows that $R_1$ and $R_1^\ast$ agree within 6\% for
all examples with $N\geq20$. The examples are chosen in pairs
(a,b), (c,d), and (e,f) to demonstrate that very different
initial states can give similar results for $R_1$ and $R_1^\ast$.
We have also compared the exact wave functions with the
trajectory of a classical particle moving in the potential
$H_{\mathrm{F}}$ with the initial conditions of the \ps\
coordinates of the jump: at the position $(q_1^\ast, q_2^\ast)$,
with momentum $(p_1^\ast, p_2^\ast)$. The results are depicted in
\fref{trajectories}. The agreement is remarkable.
\begin{figure}[bp]
\begin{center}
\includegraphics[scale=0.4]{traject.eps}
\end{center} \caption{\label{trajectories}Density plot of the final
wave function. The white dot marks the \jp\ $(q_1^{\ast},
q_2^{\ast})$ that with $(p_1^{\ast}, p_2^{\ast})$ defines the
initial conditions for the classical trajectory shown by light
line or ellipse. Here, $n=10$. The parameters $\omega_1^{\prime}$,
$\omega_2^{\prime}$ and $Q$ are listed in \protect\tref{r1}. For
the example (e), we show only one of two symmetric \jp s of equal
significance.
}
\end{figure}
\Table{\label{r1}Accuracy of the prediction by the \ps\ method of
the partition of energy between different modes for the model of
two coupled \ho s, $H_ \mathrm{I}
=\frac12({\protect\omega^{\prime}_1}^2 p_1^2 +{\protect\omega
^{\prime}_2}^2 p_2^2+(q_1-Q)^2+q_2^2)$, $H_\mathrm{F} =
\frac12(p_1^2+p_2^2+q_1^2+q_2^2)$ for examples studied earlier in
\protect\cite{Segev:2000:PSD}. Percentage of energy going to the
first mode is given by \eref{perc}.} \br
&\centre{3}{Parameters}&\centre{5}{$R_1$ (\%)}\\
\ns&\crule{3}& \crule{5}\\
Label&$\omega'_1$&$\omega'_2$&$Q$&$n=2$
&$n=6$&$n=12$&$n=20$&$n=30$\\ \mr
a&0.02&0.18&0&60.4&74.0&82.5&91.6&94.5 \\
& & & &83.3$^{\rm a}$ &92.9$^{\rm a}$ &96.2$^{\rm a}$ &97.6$^{\rm a}$ &98.4$^{\rm a}$
\\
b&10&2.2&0&71.8&87.8&93.8&96.3&97.5 \\
& & & &83.3$^{\rm a}$ &92.9$^{\rm a}$ &96.2$^{\rm a}$ &97.6$^{\rm
a}$ &98.4$^{\rm a}$
\\
c&0.45&0.01&0&25.3&10.4&5.4&3.3&2.2 \\
& & & &16.7$^{\rm a}$ &7.1$^{\rm a}$ &3.8$^{\rm a}$ &2.4$^{\rm
a}$ &1.6$^{\rm a}$
\\
d&2&18&0&24.8&10.1&5.2&3.2&2.2 \\
& & & &16.7$^{\rm a}$ &7.1$^{\rm a}$ &3.8$^{\rm a}$ &2.4$^{\rm a}$ &1.6$^{\rm a}$
\\
e&2&0.1&3&82.6&82.0&44.9&27.2&18.3 \\
& & & &83.3$^{\rm a}$ &78.4$^{\rm a}$ &42.2$^{\rm a}$ &26.1$^{\rm a}$ &17.7$^{\rm a}$
\\
f&2&10&3&82.6&82.0&44.9&27.2&18.3 \\
& & & &83.3$^{\rm a}$ &78.4$^{\rm a}$ &42.2$^{\rm a}$ &26.1$^{\rm a}$ &17.7$^{\rm a}$
\\ \br
\end{tabular}
\item[] $^{\rm a}$ The \ps\ result $R_1^\ast$.
\end{indented}
\end{table}
\section{\label{anharmsec}Anharmonicity}
In this section we study the effect of anharmonicities on the
jump. We consider anharmonic potential surfaces for the donor's
and acceptor's Hamiltonians, focusing here on Hamiltonians of \ho
s perturbed by cubic anharmonic terms as in
\eeref{finham}{iniham}. Generalization to any polynomial
anharmonicity is straightforward.
\subsection{The ground-state \wf\ for an an\ho}
The Hamiltonian of the donor, \eref{iniham}, is rewritten in this
section as
\begin{eqnarray}
H(\vec{q}, \vec{p})= \frac12\sum_{i=1}^N p_i^2+
V(\vec{q}),\nonumber\\ V(\vec{q})= \frac12\sum_{i=1}^N\omega_i^2
q_i^2 +\xi V_1(\vec{q}), \quad
V_1(\vec{q})= \frac16 \sum_{i,j,k=1}^N\frmijk_{ijk} q_i
q_j q_k, \label{iniham1}
\end{eqnarray}
where we omit the primes, although we have in mind the excited
donor surface which is marked by primed variables in the previous
and subsequent sections. We omit the subscript ``I'' in
$H_\mathrm{I}$ and $V_\mathrm{I}$ too. In (\ref{iniham1}), we
introduce a dummy expansion parameter $\xi=1$.
Rewriting the \se\ for the ground state wavefunction
$\Psi(\vec{q})$ in terms of the function $S(\vec{q})\equiv -\ln
\Psi (\vec{q})$ \cite{Bender:1978:AMM} gives the equation:
\begin{equation}
-\frac12\sum_{i=1}^N\left( \frac{\partial S}{\partial q_i}\right)
^2+\frac12\sum_{i=1}^N\frac{\partial^2 S}{\partial q_i^2}+
V(\vec{q})- E= 0. \label{ric}
\end{equation}
Without the second sum, \eref{ric} reduces to the Hamilton -
Jacobi equation for the action of a classical particle moving in
the potential $E-V(\vec{q })$. In such a quasiclassical limit, a
perturbation theory for $S(\vec{q})$ is easily developed
\cite{Popov:1994:LOE}. The more general quantum case which is
considered here is still solvable analytically, but the
corrections obtained have more monomial terms.
The function $S(\vec{q})$ and the energy $E$ are expanded in
powers of $\xi$,
\begin{eqnarray}
S(\vec{q})&=& S_{0}(\vec{q})+\xi S_{1}(\vec{q})+O(\xi ^{2}),
\label{expans} \\
E &=&E_{0}+\xi E_{1}+O(\xi ^{2}),
\end{eqnarray}
where the zero order terms are
\begin{equation}
S_{0}(\vec{q}) =\frac12\sum_{i=1}^{N}\omega _{i}q_{i}^{2}+
\mathrm{const}, \label{s0} \quad E_{0}
=\frac12\sum_{i=1}^{N}\omega _{i},
\end{equation}
and the constant in $S_{0}(\vec{q})$ is responsible for the
normalization of the wavefunction.
Let us find the first order anharmonic corrections. To the first
order in $\xi$, \eref{ric} reduces to a linear equation with
respect to the number $E_{1}$ and the function $S_{1}(\vec{q})$:
\begin{equation}
-\sum_{i=1}^{N}\frac{\partial S_{0}}{\partial
q_{i}}\frac{\partial S_{1}}{
\partial q_{i}}+\frac12\sum_{i=1}^{N}\frac{\partial ^{2}S_{1}}{\partial
q_{i}^{2}}+V_{1}(\vec{q})-E_{1}=0. \label{ric1}
\end{equation}
Let us suppose that $S_{1}(\vec{q})$ is a polynomial, then, as
soon as the inhomogeneous part of \eref{ric1},
$V_{1}(\vec{q})-E_{1}$, is a third degree polynomial, it may be
shown that the solution is at most a third degree polynomial too,
\begin{equation}
S_{1}(\vec{q})=\sum_{i=1}^{N}{A}_{i}q_{i}
+\frac12\sum_{i,j}^{N}{B}_{ij}q_{i}q_{j}
+\frac{1}{6}\sum_{i,j,k=1}^{N}{C}_{ijk}q_{i}q_{j}q_{k}\ .
\label{s1a1}
\end{equation}
We substitute (\ref{s1a1}) and (\ref{iniham1}) in \eref{ric1} and
solve to obtain:
\begin{eqnarray}
E_1=0, \quad {A}_{i} = \frac{1}{2\omega
_{i}}\sum_{j=1}^{N}\frac{\frmijk_{ijj}}{\omega _{i}+2\omega _{j}},
\quad {B}_{ij} = 0, \quad {C}_{ijk} =
\frac{\frmijk_{ijk}}{\omega _{i}+\omega _{j}+\omega _{k}}.
\label{c1}
\end{eqnarray}
Having calculated $S_1$, we would like to calculate the logarithm
of the \wf\ expanded in powers of $\xi$:
\begin{eqnarray}
\rho(\vec{q}, \vec{p})&= C^\prime \exp\left(-2 W(\vec{q},
\vec{p})\right),
\nonumber\\
W(\vec{q},\vec{p})&=W_0(\vec{q},\vec{p})+\xi W_1(\vec{q},
\vec{p})+O(\xi ^2), \label{pwig}
\end{eqnarray}
where
\begin{equation}
W_0(\vec{q},\vec{p})=\sum_{i=1}^N \exp \left(
-\frac{p_i^2}{\omega_i}-\omega_i q_i^2\right), \label{w0}
\end{equation}
$W_1(\vec{q},\vec{p})$ is the first anharmonic correction to be
determined here. Substituting the perturbed wavefunction
\begin{equation}
\Psi (\vec{q})=\left[ 1-\xi S_{1}(\vec{q})\right] \exp \left(
-S_{0}(\vec{q})\right) +O(\xi ^{2}),
\end{equation}
in the definition of the \wf, we get
\begin{eqnarray}
\rho (\vec{q},\vec{p})&\approx\left( \frac{1}{2\pi
}\right)^{N}\int d\bfeta \rme^{-\rmi\vec{p}\cdot \bfeta} \left[
1-\xi S_{1}(\vec{q}+\bfeta/2)-\xi S_{1}(\vec{q}-\bfeta/2)\right]
\nonumber \\ &\times \exp \left[
-S_{0}(\vec{q}+\bfeta/2)-S_{0}(\vec{q}-\bfeta/2) \right],
\label{wigpert}
\end{eqnarray}
and, using \eref{c1},
\begin{equation}
W_1(\vec{q},\vec{p})= \sum_{i,j}
\frac{\frmijk_{ijj}}{2\omega_i\omega_j}q_i +\sum_{i,j,k}
\frac{\frmijk_{ijk}}{\omega_i+\omega_j+\omega_k}
\left(\frac{q_{i}q_{j}q_{k}}{3}-
\frac{q_{i}p_{j}p_{k}}{\omega_j\omega_k}\right).
\label{w1explicit}
\end{equation}
The complete expression for $W=W_0+ \xi W_1$ is obtained by
substituting (\ref{w1explicit}) in \eref{wig0} and changing
variables as in \sref{math} above:
\begin{equation}\fl
W= \frac12\sum_i\alpha_i(x_i-X_i)^2+ \xi\sum_i\gamma_i(x_i-X_i)+
\frac\xi 6\sum_{ijk}\betaijk_{ijk}(x_i-X_i)(x_j-X_j)(x_k-X_k),
\label{wanh}
\end{equation}
where $x_i$ are variables collecting coordinates and momenta, and
$X_i$ are displacements. We eliminate linear terms in equation
(\ref{wanh}) by including them into effective displacements
$\tilde{X}_i= X_i+ \xi \gamma/\alpha$:
\begin{equation}
W= \frac12\sum_i\alpha _{i}\bar{x}_i^2+ \frac\xi
6\sum_{ijk}\betaijk_{ijk}\bar{x}_i\bar{x}_j\bar{x}_k+ O(\xi^2),
\label{wanh1}
\end{equation}
where $\bar{x}_i=x_i-\tilde{X}_i$. When $\frmijk_{ijj}=0$,
$\tilde{X}_{i}=X_{i}$.
%Finally, we omit terms $\sim O(\xi^2)$ and set the dummy
%perturbation parameter $\xi=1$ in equation (\ref{wanh1}).
\subsection{Anharmonic effects on the jump}
There are two effects of anharmonicities on the minimization
problem that we solve using equations~(\ref{lagr}). First, the
initial \wf\ is
shifted, redefining the effective displacement between the
two potential surfaces as $\tilde{X}_{i}$ instead of $X_{i}$.
Second, third order terms are added to the harmonic terms in both
$W$ and $H$. We treat these third-order terms by perturbation
theory, for the functions $ H=H^{(0)}+H^{(1)}\xi $,
$W=W^{(0)}+W^{(1)}\xi$, where
\begin{eqnarray} H^{(0)} &=\frac12\sum_{i}x_{i}^{2},\qquad H^{(1)}&=\frac{1}{6}
\sum_{i,j,k}\fijk_{ijk}x_{i}x_{j}x_{k}, \label{h} \\
W^{(0)} &=\frac12\sum_{i}\alpha _{i}\bar{x}_{i}^{2},\qquad
W^{(1)}&=\frac{1}{6}\sum_{i,j,k}\betaijk
_{ijk}\bar{x}_{i}\bar{x}_{j}\bar{x}_{k}. \label{w}
\end{eqnarray}
$\alpha_i,x_i,\bar{x}_i,\betaijk_{ijk}$ were defined above and
$\fijk_{ijk}$ are linear combinations of $\frmijk_{ijk}$.
\Eref{lagr} is equivalent to
\begin{eqnarray}
\alpha_{i}\bar{x}_{i}+\frac{\xi }{2}\sum_{j,k}\betaijk
_{ijk}\bar{x}_{j}\bar{ x}_{k} =\lambda \left( x_{i}+\frac{\xi
}{2}\sum_{j,k}\fijk_{ijk}x_{j}x_{k}
\right), \label{lagr1} \\
\frac12\sum_{i}x_{i}^{2}+\frac{\xi}{6}
\sum_{i,j,k}\fijk_{ijk}x_{i}x_{j}x_{k} =E. \label{lagr2}
\end{eqnarray}
The unknown variables $x_{i}$ ($i=1,...,M$) and the Lagrange
multiplier $\lambda $ are searched in the form
\begin{eqnarray}
x_{i} &=x_{i}^{(0)}+x_{i}^{(1)}\xi +o(\xi ), \label{x0} \\
\lambda &=\lambda ^{(0)}+\lambda ^{(1)}\xi +o(\xi ).
\label{lambda0}
\end{eqnarray}
In the zero order approximation ($\xi =0$), \eeref{lagr1}{lagr2})
are
\begin{eqnarray}
\alpha _{i}\bar{x}_{i}^{(0)} =\lambda ^{(0)}x_{i}^{(0)}, \nonumber \\
\frac12\sum_{i}x_{i}^{(0)2} =E, \nonumber
\end{eqnarray}
where here $\bar{x}_{i}^{(0)}=x_{i}^{(0)}-\tilde{X}_{i}$. In the
first order in $\xi $, \eeref{lagr1}{lagr2} are
\begin{equation}\fl
\alpha _{i}x_{i}^{(1)}+\frac12\sum_{j,k}\betaijk
_{ijk}\bar{x}_{j}^{(0)} \bar{x}_{k}^{(0)}= \lambda ^{(0)}\left(
x_{i}^{(1)}+\frac12
\sum_{j,k}\fijk_{ijk}x_{j}^{(0)}x_{k}^{(0)}\right) +\lambda
^{(1)}x_{i}^{(0)}, \label{lagr11}
\end{equation}
\begin{equation}
\sum_{i}x_{i}^{(0)}x_{i}^{(1)}+\frac{1}{6}
\sum_{i,j,k}\fijk_{ijk}x_{i}^{(0)}x_{j}^{(0)}x_{k}^{(0)} =0.
\label{lagr12}
\end{equation}
Let us find the first correction to the \ha\ for the two cases
discussed above using these formulas.
\underline{\bf Case (1)}
The unperturbed coordinates are given in this case by
\eeref{x1c}{xic} while the unperturbed Lagrange multiplier is
$\lambda^{(0)}=\alpha_1$. It then follows from (\ref{lagr11}) for
$i=1$ that
\begin{equation}
\lambda ^{(1)}=\frac{1}{2x_{1}^{(0)}}\sum_{j,k}\left( \betaijk
_{1jk}\bar{x} _{j}^{(0)}\bar{x}_{k}^{(0)}-\alpha
_{1}\fijk_{1jk}x_{j}^{(0)}x_{k}^{(0)}\right), \label{lambda1c2}
\end{equation}
and from (\ref{lagr11}) for $i\neq 1$ that
\begin{equation}\fl
x_{i}^{(1)}=\frac{1}{\alpha _{i}-\alpha _{1}}\biggl[ \frac12
\sum_{j,k}\left( \alpha _{1}\fijk_{ijk}x_{j}^{(0)}x_{k}^{(0)} -
\betaijk_{ijk}\bar{x }_{j}^{(0)}\bar{x}_{k}^{(0)}\right) +\lambda
^{(1)}x_{i}^{(0)}\biggr], \quad i\neq 1. \label{x1ic2}
\end{equation}
Finally, the remaining unknown variable $x_{1}^{(1)}$ can be found by
substituting (\ref{x1ic2}) into (\ref{lagr12}),
\begin{equation}
x_{1}^{(1)}= -\frac{1}{x_{1}^{(0)}}\left[ \sum_{i\neq
1}x_{i}^{(0)}x_{i}^{(1)}+ \frac16
\sum_{i,j,k}\fijk_{ijk}x_{i}^{(0)}x_{j}^{(0)}x_{k}^{(0)}\right] .
\label{x11c2}
\end{equation}
In zero order (\ha), there are two points of minimum differing by
a sign of $x_{1}$ with the same $W_{\min }$ given by
(\ref{wmin2}). In the first order approximation, $W_{\min }$ is
given by (\ref{wmin1a}) below, and it is no longer the same for
the two points, corresponding to different signs in \eref{x1c}.
So, a true minimum is the one for which (\ref{wmin1a}) is smaller.
\underline{\bf Case (2)}
In this case the unperturbed coordinates and Lagrange multiplier
are given by \eeref{coord}{minroot}. It follows from
(\ref{lagr11}) that
\begin{equation}
x_{i}^{(1)}=\frac{1}{\alpha _{i}-\lambda^{(0)}}\biggl[\frac12
\sum_{j,k}\left( \lambda^{(0)}
\fijk_{ijk}x_{j}^{(0)}x_{k}^{(0)}-\betaijk _{ijk}
\bar{x}_{j}^{(0)}\bar{x}_{k}^{(0)}\right) + \lambda
^{(1)}x_{i}^{(0)} \biggr], \label{x1a}
\end{equation}
Inserting (\ref{x1a}) into (\ref{lagr12}), we find
\begin{equation}\fl
\lambda ^{(1)}=\frac{1}{6}\left(
\sum_{i}\frac{x_{i}^{(0)2}}{\alpha _{i}-\lambda^{(0)}}\right)
^{-1}\sum_{i,j,k}\frac{x_{i}^{(0)}}{\alpha _{i}-\lambda^{(0)}
}\biggl[ 3\betaijk _{ijk}\bar{x}_{j}^{(0)}\bar{x}_{k}^{(0)}
-(2\lambda^{(0)} +\alpha _{i})\fijk_{ijk}x_{j}^{(0)}x_{k}^{(0)}
\biggr]. \label{lambda1a}
\end{equation}
In this case, we determine first $\lambda^{(1)}$ by
\eref{lambda1a}, and then $x_i^{(1)}$ by substituting
$\lambda^{(1)}$ into \eref{x1a}.
\underline{\bf $W$ in both cases}
In both cases expanding the minimum value of the constrained
function $W$ into power series $W_{\min }=W_{\min }^{(0)}+W_{\min
}^{(1)}\xi +O(\xi ^{2}) $, we find that $W_{\min }^{(0)}$ is
given by \eref{wmin}, and
\begin{equation}
W_{\min }^{(1)}=\sum_{i}\alpha
_{i}\bar{x}_{i}^{(0)}x_{i}^{(1)}+\frac{1}{6} \sum_{i,j,k}\betaijk
_{ijk}\bar{x}_{i}^{(0)}\bar{x}_{j}^{(0)}\bar{x}_{k}^{(0)}.
\label{wmin1a}
\end{equation}
\subsection{Anharmonic potentials with no effective displacements}
As a simple example, consider the case of no effective
displacements and no \dr, with $\bar{x}_{i}={x}_{i}$,
${x}_{i}^{(0)}=0$ for all $i\neq 1$ and
${x}_{1}^{(0)}=\pm\sqrt{2E}$. The unperturbed Lagrange multiplier
is $\lambda^{(0)}=\alpha_1$. It then follows that
\begin{eqnarray}
\lambda ^{(1)}&=&\frac12 {x}_{1}^{(0)} \left( \betaijk
_{111}-\alpha _{1}\fijk_{111}\right),
\\
x_{i\neq 1}^{(1)}&=&-\frac12 \frac{{x}_{1}^{(0)2}}{\alpha
_{i}-\alpha _{1}} \left( \betaijk _{i11}-\alpha
_{1}\fijk_{i11}\right),
\\
x_{1}^{(1)}&=&-\frac{1}{6}{x_{1}^{(0)2}}\fijk_{111}.
\\
W_{\min }^{(1)}&=& \frac{1}{6} {x}_{1}^{(0)3} \left( \betaijk
_{111}-\alpha _{1}\fijk_{111}\right).
\end{eqnarray}
For $\xi\ne0$ the two minima $W_{\min }^{(0)}+ \xi W_{\min
}^{(1)}$ corresponding to two different signs of $x_1^{(0)}$ are
not the same, the global minimum is reached when ${x}_{1}^{(0)}$
has the opposite sign to $\xi(\betaijk_{111}-
\alpha_1\fijk_{111})$.
\subsection{A numerical example}
The example of \fref{anhpots} shows a typical behavior of the
minimum with change of the anharmonicity.
\begin{figure}[bp]
\begin{center}
\includegraphics[scale=0.55]{anhpots.eps}
\end{center} \caption{\label{anhpots}Evolution of the minimum of
the function $W$ under the energy constraint $H=E$ with
strengthening of the anharmonicity. In $W$ and $H$, we introduced
an overall coupling parameter $\xi$: $W(x_1,x_2)= 0.4(x_1-0.1)^2+
0.6(x_2-0.2)^2+ 0.05\xi [2(x_1-0.1)^2 (x_2-0.2)- 2(x_1-0.1)
(x_2-0.2)^2+ (x_2-0.2)^3]$, $H(x_{1},x_{2})=\frac12x_1^2+\frac12
x_2^2+ 0.1\xi [-x_1^3- 3 x_1^2 x_2- 3 x_1 x_2^2+ x_2^3]$, $E=1$.
The upper left corner section refers to the \ha\ ($\xi=0$).
Position of the global minimum marked by the largest circle is a
discontinuous function of the anharmonicity since between
$\xi=0.52$ and $\xi=0.53$ as well as between $\xi=0.75$ and
$\xi=0.76$ the global minimum swaps with one of the secondary
local minima. This example shows that for strong anharmonicities
the dependence $\vec{x}^{\ast}(\xi)$ cannot be approximated by an
analytic function.}
\end{figure}
Evolution of the \jp, or the location of the minimum, is shown in
\fref{x1x2}.
\begin{figure}[bp]
\begin{center}
\includegraphics[scale=0.55]{x1x2.eps}
\end{center} \caption{\label{x1x2}
The location of the minimum $\vec{x}^{\ast}$ as a function of the
strength of the anharmonicity $\xi$ (solid line). Here, $W$, $H$,
and $E$ are the same as in \fref{anhpots}. The discontinuity is
the result of competition for the global minimum between several
local minima, see \fref{anhpots}. The dashed line is the result
of first order perturbation theory (\sref{anharmsec} of the
paper). After the first discontinuity occurs at $\xi=0.53$, this
linear approximation completely fails. There are appreciable
errors of this approximation even at $\xi>0.2$ due to presence of
quadratic terms.}
\end{figure}
For $\xi<0.1$ first order perturbation theory gives good results,
but for larger $\xi$ it breaks both because of presence of
nonlinear terms and because of an abrupt change of position of the
global minimum.
\subsection{Discussion}
The following observation can be made: when anharmonicities are
small enough to allow for a perturbative treatment their
influence on the \jp\ grows with an increased energy gap, linearly
in the energy, and increases for degrees of freedom $i$ for which
$\alpha_i\approx\alpha_1$. There is a compensating effect of the
anharmonicities on the two surfaces: their effects occur with
opposite signs. Finally, we note that not all anharmonic
potentials could be treated perturbatively, and that for
completely general potential surfaces an analytic treatment
becomes impossible. In these cases a numerical approach is needed.
\section{\label{sum}Summary and conclusions}
{Intramolecular energy transfer} is an important part of many
chemical processes. Very often, energy transfer processes in
molecules involve degrees of freedom with a separation of
timescales, for example electronic and vibrational-rotational
motion. Sometimes, the electrons transfer much energy to the
vibrations and rotations in a sudden process. This exchange is
usually followed by further intramolecular vibrational energy
transfer. The nuclei have to make a leap in coordinate or
momentum to reach the other \bos\ in such nonvertical
transitions, and the direction of the jump from one surface to
another can be very specific.
A procedure for recognizing the \jp s in \ps\ in the noncrossing
regime was introduced in Refs.~\cite{Segev:2000:PSD,
Heller:1983:RTN}. Here we have presented a closed-form complete
solution for finding the \jp\ for any transition between harmonic
potentials and a perturbative treatment of the nonharmonic
effects.
The ingredients needed to predict the quantum jump are the
classical Hamiltonian for the nuclei at the final (accepting)
electronic state and the \wf\ of the initial state $\exp{-2W}$. In
a special coordinate system, described in detail in \sref{math},
both assume the particularly simple form of \eeref{sim1}{sim2}.
In the \ha\ \eeref{wig3}{restr} replace \eeref{sim1}{sim2}
allowing for a closed form, exact solution.
In the \ha, the parameters characterizing the transition are the
energy gap $E$, the normalized displacements $\left\{
X_{i}\right\} $ between the acceptor and donor potential surfaces
and the parameters $\left\{ \alpha _{i}\right\} $.
%inversely
%proportional to the square root of the widths of the initial
%state in \ps.
For a small energy gap the transition is almost
vertical and the accepting modes are the displaced ones. For a
large energy gap the major accepting mode is the one with the
smallest $\alpha$, in general some mixed coordinate and momentum
for which the initial distribution comes closest to the final
energy surface.
The richness of possible phenomena and excitations comes about
from the transformation of coordinates from the mathematically
accessible ones to the physical coordinates of the atoms in the
molecule. Some simple examples were presented.
The harmonic case is completely solved here. Future applications
to harmonic models of specific molecules is now a straightforward
procedure. When nonharmonic potential surfaces are known, it is
possible to use the harmonic model and perturbatively correct for
the anharmonicities, only if the anharmonicities are small. More
complicated surfaces require a numerical approach.
The predictive power of the \jp\ for energy distribution between
different accepting modes was also demonstrated. Predictions for
the physical features of the radiationless transition based on
location of the \jp s deserve farther study.
\ack
This research was supported by Grant No. 2000118 from the United
States - Israel Binational Science Foundation (BSF), Jerusalem,
Israel. We gratefully acknowledge many useful discussions with
E.~J.~Heller and Shimshon Kallush. A.~S. thanks Ben Gurion
University of the Negev for hospitality.
\Bibliography{8}
\bibitem{Segev:2000:PSD} Segev B and Heller E J 2000 \JCP {\bf
112} 4004
\bibitem{Siebrand:1967:RTP} Siebrand W 1967 \JCP {\bf
46} 440
\bibitem{Berry:1970:DCP} Berry R S and Nielsen S E 1970 \PR A {\bf
1} 383
\bibitem{Jortner:1969:RTP} Jortner J, Rice S A, and Hochstrasser R M 1969 {\it Advances in Photochemistry} {\bf 7}
149
\bibitem{Hupper:1998:USE} H{\"u}pper B and Eckardt B 1998 \PR A {\bf
57} 1536
\bibitem{Heller:1978:QCC} Heller E J 1978 \JCP {\bf
68} 2066
\bibitem{Beveridge:1970:OTP} Beveridge G S G and Schechter R S, eds. 1970 {\it Optimization: Theory and
Practice} (New York: McGraw-Hill)
\bibitem{Bender:1978:AMM} Bender
C M and Orszag S A 1978 {\it Advanced Mathematical Methods for
Scientists and Engineers} (New York: McGraw-Hill)
\bibitem{Popov:1994:LOE} Popov V S and Sergeev A V 1994 {\it Phys. Lett.} A {\bf 193}
165
\bibitem{Heller:1983:RTN} Heller E J and Brown R C 1983 \JCP {\bf
79} 3336
%\bibitem{Heller:1993:QCC} Heller E J and Beck D 1993 {\it Chem. Phys. Lett.} {\bf 202} 350
\endbib
\end{document}