# Diagonalizing the Born-Oppenheimer Hamiltonian via Moyal Perturbation
Theory, Nonadiabatic Corrections and Translational Degrees of Freedom

Robert Littlejohn<EMAIL_ADDRESS>Department of Physics,
University of California, Berkeley, California 94720 USA  Jonathan Rawlinson
<EMAIL_ADDRESS>School of Mathematics, University of
Manchester, Manchester UK  Joseph Subotnik<EMAIL_ADDRESS>Department
of Chemistry, University of Pennsylvania, Philadelphia, PA, USA

###### Abstract

This article describes a method for calculating higher order or nonadiabatic
corrections in Born-Oppenheimer theory and its interaction with the
translational degrees of freedom. The method uses the Wigner-Weyl
correspondence to map nuclear operators into functions on the classical phase
space and the Moyal star product to represent operator multiplication on those
functions. The result is a power series in $\kappa^{2}$, where
$\kappa=(m/M)^{1/4}$ is the usual Born-Oppenheimer parameter. The lowest order
term is the usual Born-Oppenheimer approximation while higher order terms are
nonadiabatic corrections. These are needed in calculations of electronic
currents, momenta and densities. The method was applied to Born-Oppenheimer
theory by Littlejohn and Weigert (1993), in a treatment that notably produced
the correction $K_{22}$ to the Born-Oppenheimer Hamiltonian (see infra).
Recently Mátyus and Teufel (2019) have applied an improved and more elegant
version of the method to Born-Oppenheimer theory, and have calculated the
Born-Oppenheimer Hamiltonian for multiple potential energy surfaces to order
$\kappa^{6}$. One of the shortcomings of earlier methods is that the
separation of nuclear and electronic degrees of freedom takes place in the
context of the exact symmetries (for an isolated molecule) of translations and
rotations, and these need to be a part of the discussion. This article
presents an independent derivation of the Moyal expansion in molecular Born-
Oppenheimer theory, with special attention to the translational degrees of
freedom. We show how electronic currents and momenta can be calculated within
the framework of Moyal perturbation theory; we derive the transformation laws
of the electronic Hamiltonian, the electronic eigenstates, and the derivative
couplings under translations; we discuss in detail the rectilinear motion of
the molecular center of mass in the Born-Oppenheimer representation; and we
show how the elimination of the translational components of the derivative
couplings leads to a unitary transformation that has the effect of exactly
separating the translational degrees of freedom. Our discussion is framed in
terms of dressing transformations, that is, unitary transformations that
change the physical meanings of the operators of the theory. We show, for
example, that in the Born-Oppenheimer representation, the operator that looks
like the momentum of the nuclei actually includes, from a physical standpoint,
the momentum of the electrons. Moreover, this result is exact. In regard to
hybrid classical-quantum models such as used in surface hopping calculations,
we note that the theory presented here will be useful in developing new
semiclassical approximations going beyond the Born-Oppenheimer approximation,
as well as taking into account a glaring error of existing schemes, namely,
the failure to conserve linear and angular momentum.

††preprint: AIP/123-QED

##  I Introduction

In the Born-Oppenheimer theory of molecules the Hamiltonian becomes an
infinite-dimensional matrix of nuclear operators once an electronic basis is
chosen. Moyal perturbation theory provides a method of block-diagonalizing
this matrix to obtain effective dynamics on one or a small number of potential
energy surfaces. This is an elegant way of deriving Born-Oppenheimer dynamics
on a small number of surfaces in a manner that is more rigorous than the usual
approaches, with other advantages enumerated below.

By “Moyal perturbation theory” we refer to the use of the Wigner-Weyl
correspondence and the Moyal star product to block-diagonalize matrices of
operators. It would be more descriptive to call it “Wigner-Weyl-Moyal-van
Vleck” perturbation theory, but we will use the shorter designation. The
Wigner-Weyl correspondence maps quantum operators into functions on a
classical phase space, and the Moyal star product is used to represent
operator multiplication on those functions. The Wigner-Weyl-Moyal formalism is
fairly well known, and for convenience is summarized in Appendix A, along with
literature citations. The application in this article of Moyal perturbation
theory to Born-Oppenheimer problems is based on the treatment of Weigert and
Littlejohn (1993), with many improvements.

Moyal perturbation theory is an efficient way of obtaining higher order or
nonadiabatic corrections in Born-Oppenheimer theory. Such corrections are
needed in cases where the lowest order contribution vanishes, notably in the
calculation of electronic currents, densities, linear and angular momenta, and
matrix elements needed for radiative transitions. There is an extensive
literature on such problems, including Mead and Moscowitz (1967); Nafie
(1983); Stephens (1985); Stephens and Lowe (1985); Freedman and Nafie (1986);
Buckingham, Fowler, and Galwas (1987); Nafie (1992, 1997); Barth  _et al._
(2009); Okuyama and Takatsuka (2009); Patchkovskii (2012); Diestler (2012a, b,
2013); Diestler  _et al._ (2013); Bredtmann  _et al._ (2015); Pohl and
Tremblay (2016); Albert, Hader, and Engel (2017); Schaupp, Albert, and Engel
(2018); Schaupp and Engel (2019); Nafie (2020); Schaupp and Engel (2020);
Schaupp  _et al._ (2021); Schaupp and Engel (2022a, b). In recent years the
method of exact factorization has been developed and promoted (Abedi, Maitra,
and Gross (2010, 2012); Cederbaum (2013, 2014, 2015); Parashar, Sajeev, and
Ghosh (2015); Jecko, Sutcliffe, and Woolley (2015); Requist and Gross (2016);
Requist, Tandetzky, and Gross (2016); Curchod and Agostini (2017); Requist,
Proetto, and Gross (2017); Agostini and Curchod (2018); Lorin (2021);
Martinazzo and Burghardt (2023)). A notable application of this method is the
calculation of higher order or nonadiabatic corrections (Scherrer,
Vuilleumier, and Sabastiani (2013); Scherrer  _et al._ (2015); Schild,
Agostini, and Gross (2016)). Other approaches to higher order corrections
include Pachucki and Komasa (2008); Mátyus (2018); Mátyus and Teufel (2019);
Requist (2023). Of these, Mátyus and Teufel (2019); Requist (2023) have used a
version of what we are calling Moyal perturbation theory.

The advantages of Moyal perturbation theory, as we see them, for accessing
nonadiabatic corrections and for related topics are the following. The method
is clean and elegant, allowing one to visualize the structure of the theory as
a systematic procedure that can be carried to any order in the Born-
Oppenheimer parameter $\kappa=(m/M)^{1/4}$. The method automatically
incorporates terms in the effective Hamiltonian, such as $K_{22}$ that is
discussed below, that should be there but are usually neglected. The method
automatically accommodates large-amplitude motions, such as occur in
scattering, isomerization and photoexcitation problems, and is not restricted
to small-amplitude expansions about an equilibrium. It works in simple cases
such as motion on a single surface in the electrostatic model, as well as in
generalizations including multiple surfaces, fine-structure and spin. Many of
the first order results presented in this article, for example, the first
order corrections to the wave function and formulas for electronic momenta and
currents, are known in the case of the electrostatic model on a single
surface, but generalizations to multiple surfaces and spin are incomplete.
These are of current interest, for example, in hybrid classical-quantum models
used in surface hopping (Mannouch and Richardson (2023); Runeson and
Manolopoulos (2023); Bian  _et al._ (2021, 2023)). The method interacts well
with the exact symmetries of the system, including overall translations and
rotations. The case of rotational invariance was treated by Littlejohn,
Rawlinson, and Subotnik (2023), and translational invariance is covered in
this article. The method allows one to make simple statements about the
conservation of linear and angular momentum that are valid to all orders in
$\kappa$. The method leads to a perspective in which Born-Oppenheimer theory
consists of a series of dressing transformations, that is, unitary
transformations that change the physical meanings of the variables employed,
for example, “dressing” nuclear dynamical variables to include the effects of
the electrons. The notion of the Born-Oppenheimer transformation as a dressing
transformation was introduced by Cederbaum (2004), and we have developed the
idea considerably. We feel that it makes an important contribution to
perspectives on Born-Oppenheimer theory.

For example, in a recent study (Littlejohn, Rawlinson, and Subotnik (2023)) we
showed that in the usual Born-Oppenheimer representation the operator that
looks like the orbital angular momentum of the nuclei actually includes the
orbital angular momentum of the electrons. This is in the electrostatic model;
in fine-structure models the electron spin is included as well. These
statements are exact in the sense of being valid to all orders of the Born-
Oppenheimer expansion parameter $\kappa$. For another example, in this article
we display (in subsection V.10) a unitary transformation that has the effect
of separating the translational degrees of freedom, after which the the
operator that looks like the nuclear center of mass is actually, from a
physical standpoint, the molecular center of mass, including the electrons.

Finally, there are questions regarding three steps that may be distinguished
in the Born-Oppenheimer theory of isolated molecules, namely, the separation
of the translational degrees of freedom, the separation of the rotational
degrees of freedom, and the expansion in the Born-Oppenheimer parameter
$\kappa$. Note that the third step amounts to the separation of the nuclear
and electronic degrees of freedom. For example, there is the question of the
advantages and disadvantages of carrying out these steps in various orders. We
feel that treating the third step by Moyal perturbation theory clarifies the
analysis of this question.

The commutativity of the first and third steps (the separation of translations
and the $\kappa$-expansion) is explicitly addressed in this article. It is
sometimes stated that the Born-Oppenheimer approximation violates the
separation of translational degrees of freedom (Cederbaum (2013)), but this
does not mean that translational invariance or the associated conservation
laws are lost, simply that they take on different forms in the different
representations. Some of these are very interesting.

Moyal perturbation theory was first developed by Blount (1962a), who applied
it to various problems in condensed matter physics as well as to the Foldy-
Wouthuysen transformation (Blount (1962b); Foldy and Wouthuysen (1950);
Bjorken and Drell (1964)), the latter of which is the block-diagonalization of
the Dirac Hamiltonian. The method was rediscovered by Littlejohn and Flynn
(1991b, a), who applied it to the WKB theory of multi-component wave equations
and to several other problems in adiabatic or semiclassical theory, including
spin-orbit coupling and multi-dimensional Landau-Zener transitions (Littlejohn
and Flynn (1992, 1993)). Later Littlejohn and Weigert (1993) applied the
method to the adiabatic motion of a neutral spinning particle in an
inhomogeneous magnetic field, otherwise a Stern-Gerlach apparatus. That
problem is formally almost identical to the Born-Oppenheimer approximation in
molecules, in which the spin of the neutral particle corresponds to the
electronic degrees of freedom, while the spatial variables correspond to the
nuclear degrees of freedom. A notable outcome of that work was the production
of a second-order contribution to the Hamiltonian, what we call $K_{22}$ in
this article. This term had been missing in the analysis by other authors
(Aharonov and Stern (1992)) of the Stern-Gerlach problem.

Later Weigert and Littlejohn (1993) applied Moyal theory to the Born-
Oppenheimer approximation in molecules, producing among other things the
contribution $K_{22}$ to the Hamiltonian. This term is notable for the fact
that for large-amplitude motions it is of the same order of magnitude as the
well known diagonal Born-Oppenheimer correction, but it has received much less
attention than the latter.

Later Panati, Spohn, and Teufel (2002) and Teufel (2003) developed the
functional analysis of Moyal perturbation theory, making several applications
and innovations. One of the latter is the breaking of the problem into two
parts, the computation of the formally invariant subspace and the
representation of the dynamics on that subspace. Another is the use of
operator-valued Weyl symbols, which leads to a compact and largely coordinate-
and basis-independent formalism. This formalism was applied by Mátyus and
Teufel (2019) to Born-Oppenheimer theory, producing an expansion of the Born-
Oppenheimer Hamiltonian on one or a few potential energy surfaces out to order
$\kappa^{6}$ (described as “third order”).

The term $K_{22}$ is discussed more fully in Sec. III.10 below. Its history
and its role in Born-Oppenheimer theory have been discussed by Mátyus (2018),
whose citations place its discovery in the 1960’s, although several of the
earlier references are specific to small vibrations about an equilibrium in
diatomic or triatomic molecules. To put this term into general language,
applying to a molecule with any number of atoms, it is necessary to work in
lab or Cartesian coordinates. In addition, the importance of this term grows
with energy, and it is more important for large-amplitude motions than for
small vibrations. Therefore to treat this term in generality it is necessary
not to expand about an equilibrium position. As far as we can see the first
authors to derive this term in such generality were Moody, Shapere, and
Wilczek (1989), who used a path-integral formalism for treating the Born-
Oppenheimer approximation and its higher order corrections and who noted that
this term modifies the kinetic-energy metric, making it non-Euclidean at
second order. In addition, there is the independent discovery by Goldhaber
(2005), who derived $K_{22}$ in a discussion that links Berry’s phase,
geometry, and averaging techniques in quantum field theory. In this article we
will show that the contribution $K_{22}$ to the Hamiltonian is necessary in
the Born-Oppenheimer representation to account for the proper (rectilinear)
evolution of the center of mass of the molecule.

The original treatments of Born-Oppenheimer theory (Born and Oppenheimer
(1927); Born and Huang (1954)) and large parts of the vast literature on this
subject envision an expansion in powers of $\kappa=(m/M)^{1/4}$, in which the
nuclear displacements from an equilibrium position are small, of order
$\kappa$, while nuclear momenta are of order $\kappa^{-1}$. We call this the
“small-amplitude” ordering. It was pointed out by Littlejohn and Flynn (1992),
however, that large-amplitude motions such as occur in scattering theory
required a different ordering, in which nuclear momenta are of order
$\kappa^{-2}$. Later Mead (2006) discussed the large-amplitude ordering in
detail and its contrast to the small-amplitude ordering. Another analysis has
been given recently by Mátyus and Teufel (2019), and we have returned to the
subject in Sec. II.4 below. To summarize the situation, Moyal perturbation
theory is an expansion in powers of $\kappa^{2}$, it accommodates large-
amplitude motions, and it is equivalent to a semiclassical theory in the
nuclear degrees of freedom, that is, it is an expansion in the “nuclear
$\hbar$” (but not the electronic).

The outline of this paper is as follows. Section II contains the setup of the
problem, including notation, the relevant Hamiltonians, the electronic bases
that will be used (adiabatic, diabatic, etc.), the coordinate transformations
that bring about the exact elimination of the translational degrees of
freedom, and the several representations we will use. In this article the word
“representation” is used in a particular sense, implying a space of wave
functions, a mapping between physical states (for example, the ground state of
the molecule or an eigenstate of energy and angular momentum) and wave
functions in that space, and a mapping between physical observables (energy,
momentum, etc.) and specific linear operators that act on those wave
functions. Different representations are connected by unitary operators.

In Sec. III we describe Moyal perturbation theory in a manner that is
hopefully as painless as possible, based on the approach of Littlejohn and
Flynn (1991b); Littlejohn and Weigert (1993); Weigert and Littlejohn (1993).
Much the same territory has been covered by Mátyus and Teufel (2019), based on
the formalism of Panati, Spohn, and Teufel (2002); Teufel (2003). In
comparison, our formalism is more coordinate-based and less abstract. Our
discussion of the relation between the generator of the dressing
transformation and the usual Landau-Zener transmission probability seems to be
new, as is our presentation in subsection III.13 of the second order,
nonadiabatic corrections to the wave function. For simplicity in this section
we treat only single-surface problems.

In Sec. IV we generalize parts of the single-surface calculation in Sec. III
to multiple surfaces, and we review how the Moyal theory is applied to the
calculation of electronic momenta and currents. We offer an approach in which
the higher order corrections are obtained, not by finding corrections to the
wave function, but by transforming (“dressing”) the operators.

In Sec. V we begin by treating a number of issues regarding the translational
degrees of freedom, in a manner similar to our treatment of rotations in
Littlejohn, Rawlinson, and Subotnik (2023). In particular, we establish phase
conventions for the electronic basis states, both adiabatic and diabatic, in
effect showing how they transform under translations. We also give the
transformation laws for the electronic Hamiltonian and the derivative
couplings under translations. These results are placed in the geometrical
context of the translational fiber bundle, which we explain as a certain
decomposition of the nuclear configuration space. A notable consequence of
this work is that in the Born-Oppenheimer representations, both the original
and all the dressed versions of it, the operator that looks like the total
momentum of the nuclei actually includes, from a physical standpoint, the
momentum of the electrons. A similar statement regarding angular momentum was
a principal point of Littlejohn, Rawlinson, and Subotnik (2023).

We then address several topics concerning the translational degrees of freedom
and Moyal perturbation theory. First we study the molecular center of mass
(including the electrons) and show how its uniform, rectilinear motion is
expressed in the dressed Born-Oppenheimer representation. It turns out that
the term $K_{22}$ in the Hamiltonian is required for this to work out right.
Next we study the translational components of the electronic current and the
derivative couplings, in which the resolvent operator can be eliminated and
results obtained in closed form. The existence of the translational components
of the derivative couplings leads to some apparent paradoxes, which we resolve
by engineering a unitary transformation that removes them. This transformation
is revealed by the use of Moyal perturbation theory, but ultimately is
expressed in closed form (not as a power series). This produces an operator
that has the effect of exactly separating the translational degrees of
freedom. As far as we can tell this is a new perspective on the separation of
the translational degrees of freedom, in that the procedure can be seen as a
dressing transformation.

Finally, in Sec. VI we present some conclusions and plans for future work.

A note on notation and terminology. Starting in subsection III.4 we use hats
on nuclear operators, for example ${\hat{\mathbf{P}}}_{\alpha}$, to
distinguish them from classical quantities or $c$-numbers. See the
introduction to Sec. III. The Hamiltonian $H_{\rm BO}$ is the usual (scalar)
Born-Oppenheimer Hamiltonian on a single surface; see (53). An overbar is used
to indicate fully dressed operators, for example the Hamiltonian ${\bar{H}}$
(see (58)) or ${\bar{\mathbf{p}}}_{i}$. The scalar Hamiltonian $K$ is the
diagonal element of ${\bar{H}}$, the effective Hamiltonian that replaces
$H_{\rm BO}$ (see above (93)). In this article, “diabatic” means what is
sometimes called “quasi-diabatic,” that is, there is no implication that the
diabatic basis eliminates the derivative couplings (which is usually
impossible anyway).

##  II Preliminaries

In this section we describe the setup of the problem and establish notation.
We adopt the electrostatic model for the electronic Hamiltonian (2). For
simplicity, we treat the nuclei as distinguishable, spinless particles.

###  II.1 Hamiltonians and Notation

We consider a molecule with $N$ nuclei. The molecular Hamiltonian in the lab
frame is

$H=\sum_{\alpha=1}^{N}\frac{{\mathbf{P}}_{\alpha}^{2}}{2M_{\alpha}}+H_{e}(x;{\mathbf{r}},{\mathbf{p}}),$
(1)

where the electronic Hamiltonian is

$H_{e}(x;{\mathbf{r}},{\mathbf{p}})=\sum_{i=1}^{N_{e}}\frac{{\mathbf{p}}_{i}^{2}}{2m_{e}}+V_{\rm
Coul}({\mathbf{X}},{\mathbf{r}}).$ (2)

Here $N_{e}$ is the number of electrons; ${\mathbf{X}}_{\alpha}$,
${\mathbf{P}}_{\alpha}=-i\hbar\nabla_{\alpha}$ and $M_{\alpha}$ are
respectively the nuclear positions, momenta and masses; and ${\mathbf{r}}_{i}$
and ${\mathbf{p}}_{i}=-i\hbar\nabla_{i}$ are respectively the electron
positions and momenta. Gradient operators are
$\nabla_{\alpha}=\partial/\partial{\mathbf{X}}_{\alpha}$ and
$\nabla_{i}=\partial/\partial{\mathbf{r}}_{i}$. All positions and momenta are
taken with respect to a lab frame. Indices $\alpha$, $\beta$ label nuclei and
range over $1,\ldots,N$, while indices $i$, $j$ label electrons and range over
$1,\ldots,N_{e}$. The electron mass is $m_{e}$. The bare symbols
${\mathbf{X}}$ and ${\mathbf{r}}$ stand for the collections
$({\mathbf{X}}_{1},\ldots,{\mathbf{X}}_{N})$ and
$({\mathbf{r}}_{1},\ldots,{\mathbf{r}}_{N_{e}})$, respectively, and similarly
for the momenta ${\mathbf{P}}$ and ${\mathbf{p}}$. The potential $V_{\rm
Coul}$ contains all the Coulomb interactions among all the particles, nuclei
and electrons. The nuclear configuration space is the space upon which the
${\mathbf{X}}_{\alpha}$ are coordinates; it is $\mathbb{R}^{3N}$. We denote a
point of nuclear configuration space by the symbol $x$, whose coordinates are
${\mathbf{X}}=({\mathbf{X}}_{1},\ldots,{\mathbf{X}}_{N})$. In our usage, $x$
and ${\mathbf{X}}$ mean almost the same thing ($x$ is the point,
${\mathbf{X}}$, the coordinates of that point).

###  II.2 The Adiabatic Basis and the Working Basis

We define the “adiabatic basis” as the energy eigenbasis of the electronic
Hamiltonian,

$H_{e}(x)\,|ax;k\rangle=\epsilon_{k}(x)\,|ax;k\rangle.$ (3)

The adiabatic basis vectors in ket language are $|ax;k\rangle$, where $a$
stands for “adiabatic,” $x$ indicates the parametric dependence on the nuclear
configuration, and $k$ is the quantum number (a sequencing number for energy
eigenvalues). In (3) we write simply $H_{e}(x)$ for the electronic
Hamiltonian, indicating the $x$-dependence of the operator, which elsewhere is
written $H_{e}(x;{\mathbf{r}},{\mathbf{p}})$. The energy eigenvalues
themselves are $\epsilon_{k}(x)$. The energy eigenfunctions in wave function
language will be denoted by

$\phi_{ak}(x;{\mathbf{r}})=\langle{\mathbf{r}}|ax;k\rangle.$ (4)

Although we shall treat the index $k$ notationally as if it were discrete, in
fact above a threshold the spectrum of $H_{e}(x)$ becomes continuous and
formal sums over $k$ that employ the adiabatic basis must be interpreted as
discrete sums over the discrete spectrum plus integrals over the continuous
spectrum. In addition, one should note that the continuum threshold is a
function of $x$.

In the following we shall make use of a basis of electronic states that we
call the “working basis,” which we denote simply by $|x;k\rangle$ (without the
$a$). It is defined to be either the adiabatic basis or else a substitute,
which will be defined as the need arises. We define the working basis wave
functions by

$\phi_{k}(x;{\mathbf{r}})=\langle{\mathbf{r}}|x;k\rangle$ (5)

(again, without the $a$).

We choose the phase conventions of the adiabatic basis states so that these
states are invariant under time reversal. We also assume that the coefficients
of the linear transformation taking us from the adiabatic basis to the working
basis are real, so that the working basis states are also invariant under time
reversal. In the electrostatic model, this just means that the basis wave
functions $\phi_{ak}(x;{\mathbf{r}})$ or $\phi_{k}(x;{\mathbf{r}})$ are real.

Here and below when we use bra-ket notation for scalar products and matrix
elements, it is understood that we are integrating only over the electronic
coordinates ${\mathbf{r}}$, not the nuclear coordinates, except as noted. For
example, we have the orthonormality relations,

$\langle x;k|x;l\rangle=\delta_{kl}.$ (6)

We introduce the derivative couplings,

${\mathbf{F}}_{\alpha;kl}(x)=\langle
x;k|\nabla_{\alpha}|x;l\rangle=-(\nabla_{\alpha}\langle x;k|)|x;l\rangle,$ (7)

where the alternative form follows from the orthogonality relations (6). The
derivative couplings are defined to the extent that the basis states
$|x;k\rangle$ are differentiable with respect to $x$. The orthonormality
relations (6) imply that ${\mathbf{F}}_{\alpha;kl}(x)$, regarded as a matrix
in $kl$, is anti-Hermitian,

${\mathbf{F}}_{\alpha;kl}(x)=-{\mathbf{F}}^{*}_{\alpha;lk}(x).$ (8)

In addition, the invariance of the basis vectors under time reversal implies
that ${\mathbf{F}}_{\alpha;kl}(x)$ is real, so that as a matrix in $(kl)$, it
is real and antisymmetric. Thus, the derivative couplings vanish on the
diagonal,

${\mathbf{F}}_{\alpha;kk}(x)=0.$ (9)

###  II.3 Two Representations

We introduce two representations for the quantum state of the molecule, the
“molecular” and the “Born-Oppenheimer.” These are described in more detail in
Sec. V of Littlejohn, Rawlinson, and Subotnik (2023). In the molecular
representation the state of the molecule is represented by a wave function
$\Psi({\mathbf{X}},{\mathbf{r}})$ that depends on the coordinates of both
nuclei and electrons. In the Born-Oppenheimer representation the state of the
molecule is represented by an infinite-dimensional vector of purely nuclear
wave functions, $\psi_{k}({\mathbf{X}})$, indexed by $k$ and defined by

$\Psi({\mathbf{X}},{\mathbf{r}})=\sum_{k}\psi_{k}({\mathbf{X}})\,\phi_{k}({\mathbf{X}};{\mathbf{r}}).$
(10)

The Born-Oppenheimer representation depends on the choice made for the working
basis, and, in particular, on the phase conventions of that basis. We denote
the association between the two representations of the state of the molecule
by

$\Psi({\mathbf{X}},{\mathbf{r}})\longleftrightarrow\psi_{k}({\mathbf{X}}),$
(11)

where the notation $\longleftrightarrow$ is a reminder that the relationship
is one-to-one and that no information is lost by transforming to the Born-
Oppenheimer representation. For time-dependent problems both $\Psi$ and
$\psi_{k}$ depend on $t$, but the basis kets $|x;k\rangle$ or wave functions
$\phi_{k}({\mathbf{X}};{\mathbf{r}})$ do not.

Similarly, let $A$ be an operator acting in the molecular representation, say,
$\Psi^{\,\prime}({\mathbf{X}},{\mathbf{r}})=(A\Psi)({\mathbf{X}},{\mathbf{r}})$,
and let $\Psi\longleftrightarrow\psi_{k}$ and
$\Psi^{\,\prime}\longleftrightarrow\psi^{\,\prime}_{k}$. Then in the Born-
Oppenheimer representation $A$ is represented by an infinite-dimensional
matrix of operators $A_{kl}$ that act on nuclear wave functions, such that

$\psi^{\,\prime}_{k}({\mathbf{X}})=\sum_{l}(A_{kl}\,\psi_{l})({\mathbf{X}}).$
(12)

We denote the association between these two representations of the operator by
$A\longleftrightarrow A_{kl}$. The components $A_{kl}$ of the operator in the
Born-Oppenheimer representation are purely nuclear operators, that can be
thought of as functions of ${\mathbf{X}}$ and ${\mathbf{P}}$.

Several examples of this association are discussed in Littlejohn, Rawlinson,
and Subotnik (2023), including the following. A purely multiplicative operator
in the molecular representation depending only on the nuclear coordinates
${\mathbf{X}}_{\alpha}$, call it $f(x)$, corresponds, in the Born-Oppenheimer
representation, to a multiple of the identity matrix, that is,
$f(x)\longleftrightarrow f(x)\,\delta_{kl}$. In particular, this applies to
the nuclear coordinates, the components of ${\mathbf{X}}_{\alpha}$. A purely
electronic operator in the molecular representation, that is, a function of
${\mathbf{r}}$ and ${\mathbf{p}}$, corresponds to a purely multiplicative
operator in the Born-Oppenheimer representation, that is, a matrix whose
matrix elements are functions of ${\mathbf{X}}$ only. These are nothing but
the matrix elements of the original operator in the working basis. For
example, in the case of the momentum ${\mathbf{p}}_{i}$ of one of the
electrons, we have ${\mathbf{p}}_{i}\longleftrightarrow\langle
x;k|{\mathbf{p}}_{i}|x;l\rangle$. This also applies to electronic operators
that are parameterized by $x$, such as the electronic Hamiltonian. We define

$\langle x;k|H_{e}(x)|x;l\rangle=W_{kl}(x),$ (13)

so that $H_{e}(x)\longleftrightarrow W_{kl}(x)$. Note that if the working
basis is the adiabatic basis (at least for some range of indices $k$, $l$),
then

$W_{kl}(x)=\epsilon_{k}(x)\,\delta_{kl}.$ (14)

The momentum ${\mathbf{P}}_{\alpha}=-i\hbar\,\nabla_{\alpha}$, which in the
molecular representation represents the kinetic momentum of nucleus $\alpha$,
is transformed to the Born-Oppenheimer representation according to

${\mathbf{P}}_{\alpha}\longleftrightarrow{\mathbf{P}}_{\alpha}\,\delta_{kl}-i\hbar\,{\mathbf{F}}_{\alpha;kl}(x).$
(15)

Equation (15) illustrates an important theme of this article. The operator
${\mathbf{P}}_{\alpha}$ is the differential operator $-i\hbar\nabla_{\alpha}$
in all representations, but it only represents physically the kinetic momentum
of nucleus $\alpha$ in the molecular representation. In the Born-Oppenheimer
representation, there is a correction term, as shown (which, it turns out, is
related to the electronic momentum).

Now suppose that we have operators $A$, $B$ and $C$, and that
$A\longleftrightarrow A_{kl}$, $B\longleftrightarrow B_{kl}$ and
$C\longleftrightarrow C_{kl}$. Then we have some theorems. First, if
$B=A^{\dagger}$, then $B_{kl}=(A_{lk})^{\dagger}$. Next, if $C=AB$, then

$C_{kl}=\sum_{p}A_{kp}\,B_{pl}.$ (16)

Next, if $\Psi_{1}\longleftrightarrow\psi_{1k}$ and
$\Psi_{2}\longleftrightarrow\psi_{2k}$, then

$\langle\Psi_{1}|\Psi_{2}\rangle_{Xr}=\sum_{k}\langle\psi_{1k}|\psi_{2k}\rangle_{X}$
(17)

where the subscripts $X$, $Xr$ indicate the variables that are integrated over
in taking the scalar product. (If no subscript is given, integration over
${\mathbf{r}}$ alone is implied, as in (6).) Finally, matrix elements of
operators transform between the two representations according to

$\langle\Psi_{1}|A|\Psi_{2}\rangle_{Xr}=\sum_{kl}\langle\psi_{1k}|A_{kl}|\psi_{2l}\rangle_{X}.$
(18)

As for the molecular Hamiltonian (1), when transformed to the Born-Oppenheimer
representation, $H\longleftrightarrow H_{kl}$, it becomes

$H_{kl}=\sum_{\alpha=1}^{N}\frac{1}{2M_{\alpha}}\sum_{p}\bigl{[}{\mathbf{P}}_{\alpha}\,\delta_{kp}-i\hbar\,{\mathbf{F}}_{\alpha;kp}(x)\bigr{]}\cdot\bigl{[}{\mathbf{P}}_{\alpha}\,\delta_{pl}-i\hbar\,{\mathbf{F}}_{\alpha;pl}(x)\bigr{]}+W_{kl}(x),$
(19)

where the sum on $p$ indicates a multiplication of matrices of operators, as
in (16). We expand this according to

$H_{kl}=H_{0,kl}+\hbar\,H_{1,kl}+\hbar^{2}\,H_{2,kl},$ (20)

where

$\displaystyle H_{0,kl}$ $\displaystyle=$
$\displaystyle\left(\sum_{\alpha=1}^{N}\frac{{\mathbf{P}}_{\alpha}^{2}}{2M_{\alpha}}\right)\delta_{kl}+W_{kl}({\mathbf{X}}),$
(21a) $\displaystyle H_{1,kl}$ $\displaystyle=$
$\displaystyle-i\sum_{\alpha=1}^{N}\frac{1}{2M_{\alpha}}[{\mathbf{P}}_{\alpha}\cdot{\mathbf{F}}_{\alpha;kl}({\mathbf{X}})+{\mathbf{F}}_{\alpha;kl}({\mathbf{X}})\cdot{\mathbf{P}}_{\alpha}],$
(21b) $\displaystyle H_{2,kl}$ $\displaystyle=$
$\displaystyle-\sum_{\alpha=1}^{N}\sum_{p}\frac{1}{2M_{\alpha}}{\mathbf{F}}_{\alpha;kp}({\mathbf{X}})\cdot{\mathbf{F}}_{\alpha;pl}({\mathbf{X}}).$
(21c)

Note that, since ${\mathbf{F}}_{\alpha;kl}(x)$ is real, the two terms in the
square brackets in (21b) are Hermitian conjugates of each other, so the sum is
Hermitian. Overall, we have $(H_{1,kl})^{\dagger}=H_{1,lk}$, so that $H_{1}$,
as a matrix of nuclear operators, is Hermitian.

The molecular Schrödinger equation, which in the molecular representation is

$H\,\Psi({\mathbf{X}},{\mathbf{r}})=E\,\Psi({\mathbf{X}},{\mathbf{r}}),$ (22)

or, in its time-dependent version,

$H\,\Psi({\mathbf{X}},{\mathbf{r}},t)=i\hbar\frac{\partial\Psi({\mathbf{X}},{\mathbf{r}},t)}{\partial
t},$ (23)

becomes, in the Born-Oppenheimer representation,

$\sum_{l}H_{kl}\,\psi_{l}({\mathbf{X}})=E\,\psi_{k}({\mathbf{X}}),$ (24)

or, in its time-dependent version,

$\sum_{l}H_{kl}\,\psi_{l}({\mathbf{X}},t)=i\hbar\frac{\partial\psi_{k}({\mathbf{X}},t)}{\partial
t}.$ (25)

.

###  II.4 Born-Oppenheimer Ordering

The standard Born-Oppenheimer ordering parameter is $\kappa=(m_{e}/M)^{1/4}$,
where $M$ is an average or typical nuclear mass (Born and Oppenheimer (1927);
Born and Huang (1954)). By this definition, $\kappa$ has the numerical value
of roughly $10^{-1}$. A different approach is to scale the nuclear masses by
$\kappa^{-4}$, that is, to replace the masses $M_{\alpha}$ in the Hamiltonian
by $M_{\alpha}/\kappa^{4}$, where $\kappa$ is a scaling parameter that ranges
from $\kappa=0$ (for complete adiabatic separation of the length and time
scales) to $\kappa=1$ (which gives the physical values). Actually, the limit
$\kappa\to 0$ is singular since the nuclear de Broglie wavelength goes to zero
and the nuclear wave function is no longer defined. But many physical
quantities are well defined in the limit $\kappa\to 0$. By making $\kappa$ a
variable the Hamiltonian becomes parameterized by $\kappa$, and thereby
represents a family of dynamical systems. We can think of a curve through the
space of Hamiltonians, parameterized by $\kappa$. Then we have a solvable or
at least simplified system as $\kappa\to 0$.

In the following we shall adopt the second approach, in which $\kappa$ is a
variable. Then we can speak of the order of magnitude of various quantities in
terms of the power of $\kappa$ by which they scale as $\kappa\to 0$. For
example, since the electronic Hamiltonian does not depend on the nuclear
masses, it is independent of $\kappa$ and everything associated with it (the
energy eigenvalues $\epsilon_{k}(x)$, the derivative couplings
${\mathbf{F}}_{\alpha;kl}(x)$, etc.), are independent of $\kappa$ and
therefore of order $\kappa^{0}=1$.

To say a quantity is of order $\kappa^{0}=1$ in this sense does not
necessarily imply that its numerical value is close to unity, since there may
be other reasons (besides the values of the nuclear masses or momenta) why it
is large or small. For example, in a scattering problem
${\mathbf{X}}_{\alpha}$ may take on large values, independent of $\kappa$, and
the derivative couplings in the adiabatic basis diverge as a degeneracy or
conical intersection is approached. In terms of the $\kappa$ scaling, however,
these quantities are of order unity. On the other hand, if there is no other
reason why a given quantity should be large or small, and if it is of order
$\kappa^{n}$ in the scaling sense, then it should have a numerical value in
atomic units of roughly $10^{-n}$.

In addition to the scaling of the nuclear masses, we must also indicate how
the dynamical variables scale with $\kappa$ as we take $\kappa\to 0$. For
small vibrations about an equilibrium, the displacement from equilibrium
${\mathbf{X}}_{\alpha}-{\mathbf{X}}_{0,\alpha}$ is of order
$\kappa^{1}=\kappa$, as discussed by Born and Huang (1954), while the nuclear
momentum ${\mathbf{P}}_{\alpha}$ is order $\kappa^{-1}$ and the nuclear
velocity ${\mathbf{V}}_{\alpha}={\mathbf{P}}_{\alpha}/M_{\alpha}$ is of order
$\kappa^{3}$. Thus, the time require to execute an oscillation of amplitude of
order $\kappa$ is $\kappa^{-2}$, and the vibrational frequency is of order
$\kappa^{2}$. In this analysis we are imagining that as we move through the
space of dynamical systems by taking $\kappa\to 0$, we are also changing the
space of wave functions or quantum states that we are considering. For the
problem of small vibrations, we take this limit in such a way that the
vibrational quantum number is constant, which causes the spatial extent of the
wave function to shrink in proportion to $\kappa$. A nice discussion of this
(small-amplitude) ordering is given by Mead and Moscowitz (1967).

Modern applications such as isomerization, scattering and photoexcitation
problems involve motions in which the nuclei are not restricted to a small
range around an equilibrium, but rather move over distances of the order of an
atomic unit or larger. For these applications a different assumption about the
scaling of the dynamical variables is required. The modified scaling
assumptions required for such applications have been discussed in detail by
Mead (2006). See also Mátyus and Teufel (2019). In the following we add to
this discussion and show how it applies to our approach to Moyal perturbation
theory.

We shall call the modified scaling law the “large-amplitude” scaling. In this
scaling we imagine that the kinetic or potential energy of the nuclei remains
constant as $\kappa\to 0$. For example, if a nucleus climbs half way up a
potential well when $\kappa=1$, we require that it continue to do so as we
take $\kappa\to 0$. We do this because we want the qualitative features of the
simplified problem as $\kappa\to 0$ to be the same as those of the physical
problem at $\kappa=1$. This means, for example, that bound state quantum
numbers, to the extent those are meaningful for large-amplitude motions, go to
infinity as $\kappa\to 0$.

If the kinetic energy ${\mathbf{P}}_{\alpha}^{2}/2M_{\alpha}$ is independent
of $\kappa$ as $\kappa\to 0$, as required by the large-amplitude ordering,
then, given the scaling of $M_{\alpha}$, we conclude that
${\mathbf{P}}_{\alpha}$ is order $\kappa^{-2}$ (therefore larger by a factor
of $1/\kappa$ than in the small-amplitude scaling). Nuclear velocities
${\mathbf{V}}_{\alpha}={\mathbf{P}}_{\alpha}/M_{\alpha}$ are therefore of
order $\kappa^{2}$ and the time required for a nucleus to travel a distance of
the order of one atomic unit is of order $\kappa^{-2}$. Thus nuclear motions
are executed ever more slowly as $\kappa\to 0$, as expected in a system that
can be considered “adiabatic.”

The spatial scale of an electronic energy eigenfunction
$\phi_{ak}({\mathbf{X}};{\mathbf{r}})$ for fixed $k$, that is, the electronic
de Broglie wave length, is of order $\kappa^{0}=1$, since the electronic
Hamiltonian is independent of the nuclear masses, as explained. The same
applies to the operator $\nabla_{i}$. The same also applies to the
${\mathbf{X}}$-dependence of the electronic energy eigenfunctions and the
operator $\nabla_{\alpha}$ when acting on them. We shall assume that the
transformation from the adiabatic basis to the working basis is independent of
$\kappa$, so the same scale lengths (in both ${\mathbf{r}}$ and
${\mathbf{X}}$) apply to the working basis functions
$\phi_{k}({\mathbf{X}};{\mathbf{r}})$. The nuclear wave functions
$\psi_{k}({\mathbf{X}})$ in the Born-Oppenheimer representation, however, have
a de Broglie wave length of order $\kappa^{2}$, since the operator
${\mathbf{P}}_{\alpha}=-i\hbar\,\nabla_{\alpha}$ is of order $\kappa^{-2}$.
Thus, in the expansion (10) the scale lengths of the two factors on the right
($\psi_{k}$ and $\phi_{k}$) with respect to ${\mathbf{X}}$ are quite
different.

In the large-amplitude scaling, the nuclear de Broglie wave length is of order
$\kappa^{2}$ times smaller than the scale length of the environment (the
potential, derivative couplings, etc.) in which the nuclei move. These are
otherwise the conditions for WKB theory, which would be more useful for
molecular dynamics were it not for the usual difficulties of WKB theory in
multidimensional problems (chaotic orbits and the lack of global, invariant
Lagrangian manifolds in phase space). Nevertheless, the basic method of this
paper, the use of Moyal perturbation theory, is a part of semiclassical
theory.

In this paper we use the large-amplitude scaling exclusively. For simplicity
we will not insert explicit factors of $\kappa^{n}$ to indicate the order of
various terms. Instead, the ordering of any term can be obtained by applying
the two rules, $M_{\alpha}=O(\kappa^{-4})$ and
${\mathbf{P}}_{\alpha}=O(\kappa^{-2})$.

In the large-amplitude ordering, the three terms in (21), $H_{0,kl}$,
$H_{1,kl}$ and $H_{2,kl}$, are of order $\kappa^{0}=1$, $\kappa^{2}$ and
$\kappa^{4}$, respectively. This ordering will be the basis of the
perturbation theory we will carry out in this paper. It is not, however,
ordinary perturbation theory, which is based solely on the magnitudes of the
unperturbed Hamiltonian and the perturbation, because the nuclear wave
functions also have a dependence on $\kappa$, that is, the nuclear de Broglie
wave length is of order $\kappa^{2}$.

###  II.5 Translational Degrees of Freedom

The Hamiltonian (1) commutes with translations and thus conserves total linear
momentum. The manner in which this invariance and conservation law are
manifested is straightforward in the molecular representation but more obscure
in the Born-Oppenheimer representation. Therefore one is motivated to carry
out the exact elimination of the translational degrees of freedom before
switching to the Born-Oppenheimer representation. This is a standard procedure
(Sutcliffe (2000)) but it requires some algebra and attention to notation.

We begin with the nuclear configuration space, the $3N$-dimensional space upon
which the coordinates are ${\mathbf{X}}_{\alpha}$, $\alpha=1,\ldots,N$. This
is the space upon which purely nuclear wave functions $\psi({\mathbf{X}})$ are
defined. A coordinate transformation on this space that is useful for
separating the translational degrees of freedom is
${\mathbf{X}}_{\alpha}\to({\mathbf{X}}_{\rm CM},{\mathbf{Y}}_{\kappa})$, where
${\mathbf{X}}_{\rm CM}$ is the center of mass of the nuclei,

${\mathbf{X}}_{\rm
CM}=\frac{1}{M_{t}}\sum_{\alpha=1}^{N}M_{\alpha}{\mathbf{X}}_{\alpha},$ (26)

and where ${\mathbf{Y}}_{\kappa}$, $\kappa=1,\ldots,N-1$, are a set of Jacobi
vectors (see Appendix B). Here $M_{t}$ is the total nuclear mass,

$M_{t}=\sum_{\alpha=1}^{N}M_{\alpha}.$ (27)

We let the momenta conjugate to $({\mathbf{X}}_{\rm
CM},{\mathbf{Y}}_{\kappa})$ be $({\mathbf{P}}_{t},{\mathbf{Q}}_{\kappa})$.
Then it turns out that ${\mathbf{P}}_{t}$, the momentum conjugate to
${\mathbf{X}}_{\rm CM}$, is the total momentum of the nuclei,

${\mathbf{P}}_{t}=\sum_{\alpha=1}^{N}{\mathbf{P}}_{\alpha}.$ (28)

The other details of this coordinate transformation are summarized in Appendix
B.

When this transformation is applied to the nuclear kinetic energy, we find

$\sum_{\alpha=1}^{N}\frac{{\mathbf{P}}_{\alpha}^{2}}{2M_{\alpha}}=\frac{{\mathbf{P}}_{t}^{2}}{2M_{t}}+\sum_{\kappa=1}^{N-1}\frac{{\mathbf{Q}}_{\kappa}^{2}}{2\mu_{\kappa}},$
(29)

where $\mu_{\kappa}$ are reduced masses associated with Jacobi coordinates
${\mathbf{Y}}_{\kappa}$.

We require a different coordinate transformation on the molecular
configuration space, the $3(N+N_{e})$-dimensional space with coordinates
$({\mathbf{X}}_{\alpha},{\mathbf{r}}_{i})$ upon which the molecular wave
function $\Psi({\mathbf{X}},{\mathbf{r}})$ is defined. This can be seen as the
composition of two transformations,
$({\mathbf{X}}_{\alpha},{\mathbf{r}}_{i})\to({\mathbf{X}}_{\rm
CM},{\mathbf{Y}}_{\kappa},{\mathbf{r}}_{i})\to({\mathbf{R}}_{\rm
CM},{\mathbf{Y}}_{\kappa},{\mathbf{s}}_{i})$, where in the first step we
transform the nuclear coordinates as described above, leaving the electronic
coordinates ${\mathbf{r}}_{i}$ alone, while in the second step we leave the
nuclear Jacobi vectors ${\mathbf{Y}}_{\kappa}$ alone while introducing two new
sets of coordinates defined by

${\mathbf{R}}_{\rm CM}=\frac{M_{t}{\mathbf{X}}_{\rm
CM}+N_{e}m_{e}{\mathbf{r}}_{\rm CM}}{M_{\rm mol}}$ (30)

and

${\mathbf{s}}_{i}={\mathbf{r}}_{i}-{\mathbf{X}}_{\rm CM}.$ (31)

Here $M_{\rm mol}$ is the total molecular mass,

$M_{\rm mol}=M_{t}+N_{e}m_{e},$ (32)

and ${\mathbf{r}}_{\rm CM}$ is the center of mass of the electrons,

${\mathbf{r}}_{\rm CM}=\frac{1}{N_{e}}\sum_{i=1}^{N_{e}}{\mathbf{r}}_{i}.$
(33)

Evidently ${\mathbf{R}}_{\rm CM}$ is the molecular center of mass (not to be
confused with ${\mathbf{X}}_{\rm CM}$, the nuclear center of mass) and
${\mathbf{s}}_{i}$ is the position of electron $i$ relative to the nuclear
center of mass.

We denote the momenta conjugate to coordinates $({\mathbf{R}}_{\rm
CM},{\mathbf{Y}}_{\kappa},{\mathbf{s}}_{i})$ by $({\mathbf{P}}_{\rm
CM},{\mathbf{Q}}_{\kappa},{\mathbf{q}}_{i})$, where ${\mathbf{Q}}_{\kappa}$,
the momentum conjugate to Jacobi vector ${\mathbf{Y}}_{\kappa}$, is the same
as produced by the first (purely nuclear) transformation (see (250)). As for
${\mathbf{P}}_{\rm CM}$, the momentum conjugate to ${\mathbf{R}}_{\rm CM}$, it
turns out to be the total linear momentum of the molecule,

${\mathbf{P}}_{\rm CM}={\mathbf{P}}_{t}+{\mathbf{p}}_{t},$ (34)

where ${\mathbf{P}}_{t}$, the total nuclear momentum, is given by (28), and
${\mathbf{p}}_{t}$, the total electronic momentum, is given by

${\mathbf{p}}_{t}=\sum_{i=1}^{N_{e}}{\mathbf{p}}_{i}.$ (35)

The other details of this coordinate transformation are given in Appendix C.

When we transform the molecular Hamiltonian (1) to the new coordinates, we
obtain

$H=\frac{{\mathbf{P}}_{\rm CM}^{2}}{2M_{\rm mol}}+H_{\rm
TR}({\mathbf{Y}},{\mathbf{s}},{\mathbf{Q}},{\mathbf{q}}),$ (36)

where the “translation-reduced” molecular Hamiltonian is

$H_{\rm
TR}({\mathbf{Y}},{\mathbf{s}},{\mathbf{Q}},{\mathbf{q}})=\sum_{\kappa=1}^{N-1}\frac{{\mathbf{Q}}_{\kappa}^{2}}{2\mu_{\kappa}}+H_{{\rm
TR}e}(y;{\mathbf{s}},{\mathbf{q}}),$ (37)

and the translation-reduced electronic Hamiltonian is

$H_{{\rm
TR}e}(y;{\mathbf{s}},{\mathbf{q}})=\sum_{i=1}^{N_{e}}\frac{{\mathbf{q}}_{i}^{2}}{2m_{e}}+\sum_{i,j=1}^{N_{e}}\frac{{\mathbf{q}}_{i}\cdot{\mathbf{q}}_{j}}{2M_{t}}+V_{\rm
Coul}({\mathbf{Y}},{\mathbf{s}}).$ (38)

Here $y$ means a point of the translation-reduced configuration space whose
coordinates are ${\mathbf{Y}}_{\kappa}$, and $V_{\rm Coul}$ is independent of
${\mathbf{R}}_{\rm CM}$ because it is translationally invariant. The
translation-reduced configuration space will be discussed more fully in
subsection V.2 (see also Fig. 1). The middle term on the right in (38) is the
mass-polarization term.

###  II.6 The Translation-Reduced Molecular and Born-Oppenheimer
Representations

The Hamiltonian $H_{\rm TR}$ in (37) gives us a “translation-reduced”
molecular representation, in which the wave function $\Psi_{\rm
TR}({\mathbf{Y}},{\mathbf{s}})$ is defined by

$\Psi({\mathbf{X}},{\mathbf{r}})=\Psi({\mathbf{R}}_{\rm
CM},{\mathbf{Y}},{\mathbf{s}})=\exp(i{\mathbf{K}}\cdot{\mathbf{R}}_{\rm
CM})\,\Psi_{\rm TR}({\mathbf{Y}},{\mathbf{s}}),$ (39)

which is an eigenfunction of ${\mathbf{P}}_{\rm CM}$ with eigenvalue
$\hbar{\mathbf{K}}$. This is in the time-independent version; in the time-
dependent version we define $\Psi_{\rm TR}$ by

$\Psi({\mathbf{X}},{\mathbf{r}},t)=\Psi({\mathbf{R}}_{\rm
CM},{\mathbf{Y}},{\mathbf{s}},t)=\exp[i({\mathbf{K}}\cdot{\mathbf{R}}_{\rm
CM}-\hbar{\mathbf{K}}^{2}t/2M_{\rm mol})]\,\Psi_{\rm
TR}({\mathbf{Y}},{\mathbf{s}},t).$ (40)

The wave function $\Psi_{\rm TR}$ satisfies a translation-reduced version of
the Schrödinger equation,

$H_{\rm TR}\,\Psi_{\rm TR}({\mathbf{Y}},{\mathbf{s}})=E_{\rm TR}\,\Psi_{\rm
TR}({\mathbf{Y}},{\mathbf{s}}),$ (41)

or its time-dependent version,

$H_{\rm TR}\,\Psi_{\rm
TR}({\mathbf{Y}},{\mathbf{s}},t)=i\hbar\,\frac{\partial\Psi_{\rm
TR}({\mathbf{Y}},{\mathbf{s}},t)}{\partial t}.$ (42)

Once we have achieved the translation-reduced molecular representation of the
system we can expand in an electronic basis and create a translation-reduced
Born-Oppenheimer representation, in which the state of the system is
represented by an infinite-dimensional vector of nuclear wave functions
$\psi_{{\rm TR},k}({\mathbf{Y}})$, defined on the translation-reduced
configuration space. The potential energy surfaces defined by the eigenvalues
$\epsilon_{{\rm TR},k}(y)$ are not the same as before, because the
translation-reduced electronic Hamiltonian (37) contains mass-polarization
terms, in contrast to the original electronic Hamiltonian (2), whose
eigenvalues are $\epsilon_{k}(x)$.

|  |  | Dressed
---|---|---|---
| Molecular | Born-Opp. | Born-Opp.
| $\Psi({\mathbf{X}},{\mathbf{r}})$ | $\psi_{k}({\mathbf{X}})$ | ${\bar{\psi}}_{k}({\mathbf{X}})$
Lab | $H$ | $H_{kl}$ | ${\bar{H}}_{kl}$
| (1) | (10), (11) | (58)
| $\Psi_{\rm TR}({\mathbf{Y}},{\mathbf{s}})$ | $\psi_{{\rm TR},k}({\mathbf{Y}})$ | ${\bar{\psi}}_{{\rm TR},k}({\mathbf{Y}})$
Translation- | $H_{\rm TR}$ | $H_{{\rm TR},kl}$ | ${\bar{H}}_{{\rm TR},kl}$
Reduced | (37), (39), (40) | — | —
Table 1: Six representations of the quantum mechanics of a molecule. In each
case, the relevant wave function and Hamiltonian are shown, as well as
references to relevant equations. In the molecular representation (first
column) the wave function $\Psi$ or $\Psi_{\rm TR}$ is a function defined on
the combined nuclear-electronic configuration space. In the Born-Oppenheimer
representation (second column), the wave function $\psi_{k}$ or $\psi_{{\rm
TR},k}$ is an infinite-dimensional vector of wave functions defined on the
nuclear configuration space. The same is true of the dressed representations
(third column); in this representation, the Hamiltonian is block-diagonal.
Sometimes we will call the Born-Oppenheimer representation (second column) the
“original” Born-Oppenheimer representation, to distinguish it from the dressed
representations in the third column.

Thus we now have several representations of the quantum mechanics of the
molecule, including the original or “lab” molecular and Born-Oppenheimer
representations, and now “translation-reduced” versions of these. The
situation is summarized in Table 1. Once we have achieved the Born-Oppenheimer
representation, either the lab or translation-reduced, the Hamiltonian is an
infinite-dimensional matrix of nuclear operators, which we can block-
diagonalize by means of Moyal perturbation theory. This carries us to the last
column of Table 1, the “dressed” or block-diagonalized Born-Oppenheimer
representation. We turn now to a detailed description of this last step.

##  III Moyal Perturbation Theory on a Single Surface

The purpose of this section is to describe in detail the Moyal method for the
block diagonalization of the Born-Oppenheimer Hamiltonian. For a somewhat
different formalism, see Panati, Spohn, and Teufel (2002); Mátyus and Teufel
(2019).

For this purpose we can work either with the lab representations or the
translation-reduced. In some ways it is preferable to carry out the exact
separation of the translational degrees of freedom first, and then to block-
diagonalize the Hamiltonian. But much of the current literature, for example,
that devoted to surface hopping, uses the lab representations (so also the
treatment of Mátyus and Teufel (2019)). Therefore in this section we shall
illustrate Moyal diagonalization of the Born-Oppenheimer Hamiltonian in the
lab representations, effectively working across the first row of Table 1 from
left to right. For simplicity we will do this for a single surface problem,
deferring multiple surfaces to subsection IV.1. We select one surface $k=n$,
so that $n$ is fixed for the duration of this section (Sec. III).

In subsections III.1–III.3 we use unaccented symbols for operators, for
example, $H$ for the Hamiltonian. Starting in subsection III.4, however, and
continuing to the end of the paper, it is necessary to distinguish between
operators, which are represented with a hat, and their Weyl symbols, which are
$c$-number functions on the classical phase space and which are denoted
without the hat. For example, starting in subsection III.4, the Hamiltonian
operator is ${\hat{H}}$ while its Weyl symbol is $H$. Weyl symbols are
explained in Appendix A.

###  III.1 Basis Vectors and the Privileged Subspace

We let $|x;n\rangle=|ax;n\rangle$, so that in the one case $k=n$ the working
basis vector and the adiabatic basis vector are the same. We define the
“privileged subspace” of the electronic Hilbert space as

${\mathcal{S}}(x)=\mathop{\rm span}\\{|ax;n\rangle\\}$ (43)

and we define ${\mathcal{S}}^{\perp}(x)$ as the complementary, orthogonal
subspace. The working basis vectors $|x;k\rangle$ for $k\neq n$ are either the
adiabatic basis vectors $|ax;k\rangle$ for $k\neq n$ or another choice to be
specified. Both the working basis vectors $|x;k\rangle$ and the adiabatic
basis vectors $|ax;k\rangle$ for $k\neq n$ span ${\mathcal{S}}^{\perp}(x)$, so
they are linear combinations of one another. The subspace ${\mathcal{S}}(x)$
is one-dimensional, while ${\mathcal{S}}^{\perp}(x)$ is infinite-dimensional.
The projection operator onto the privileged subspace is

$P(x)=|ax;n\rangle\langle ax;n|,$ (44)

while the complementary projector is $Q(x)=1-P(x)$. Equation (44) gives $P(x)$
in the molecular representation; on transforming to the Born-Oppenheimer
representation it becomes

$P(x)\longleftrightarrow P_{kl}=\delta_{nk}\,\delta_{nl}.$ (45)

###  III.2 The Born-Oppenheimer Approximation

Since $|x;n\rangle=|ax;n\rangle$, the diagonal matrix element $W_{nn}$ of the
electronic Hamiltonian satisfies

$W_{nn}(x)=\epsilon_{n}(x).$ (46)

In addition, since $|x;n\rangle$ is orthogonal to ${\mathcal{S}}^{\perp}(x)$,
we have $W_{nk}(x)=W_{kn}(x)=0$ for $k\neq n$, and the matrix $W_{kl}(x)$ is
block-diagonal. The blocks in question are the $1\times 1$ block and the
infinite-dimensional block corresponding to ${\mathcal{S}}(x)$ and
${\mathcal{S}}^{\perp}(x)$, respectively. An operator that is block-diagonal
does not couple ${\mathcal{S}}(x)$ with ${\mathcal{S}}^{\perp}(x)$. The terms
in the Hamiltonian (21) that do couple ${\mathcal{S}}(x)$ and
${\mathcal{S}}^{\perp}(x)$ are the off-block-diagonal elements of $H_{1}$ and
$H_{2}$, that is, elements with indices $(nk)$ or $(kn)$ for $k\neq n$. These
matrix elements involve the derivative couplings and are small in the $\kappa$
ordering, either of order $\kappa^{2}$ or $\kappa^{4}$. In particular, they
are small in comparison to the diagonal element at zeroth order,

$H_{0,nn}=\sum_{\alpha=1}^{N}\frac{{\mathbf{P}}_{\alpha}^{2}}{2M_{\alpha}}+\epsilon_{n}(x),$
(47)

which is of order $\kappa^{0}$.

In the literature the “Born-Oppenheimer approximation” usually means
approximating the molecular wave function as

$\Psi({\mathbf{X}},{\mathbf{r}})=\psi({\mathbf{X}})\,\phi_{an}({\mathbf{X}};{\mathbf{r}}),$
(48)

that is, a version of (10) with a single term $k=n$ and employing the
adiabatic, energy eigenfunction $\phi_{an}({\mathbf{X}};{\mathbf{r}})$. But if
we throw away the small, off-block-diagonal terms of the molecular Hamiltonian
(21), then the molecular Schrödinger equation (22) has solutions of the form
(48), that is, (24) has solutions in which $\psi_{k}$ is nonzero only in the
one slot $k=n$. The same is true of the time-dependent versions of the
molecular Schrödinger equation, (23) and (25), which have solutions of the
form (48) in which both $\Psi({\mathbf{X}},{\mathbf{r}},t)$ and
$\psi({\mathbf{X}},t)$ depend on time but the adiabatic basis state
$\phi_{an}({\mathbf{X}};{\mathbf{r}})$ does not.

Therefore we shall regard the “Born-Oppenheimer approximation,” in the context
of a single-surface problem, as one of throwing away the off-block-diagonal
matrix elements of $H_{kl}$. In this approximation, the molecular Schrödinger
equation (24) or (25) on the privileged subspace ${\mathcal{S}}(x)$ becomes

$H_{\rm
BO}({\mathbf{X}},{\mathbf{P}})\,\psi({\mathbf{X}})=E\,\psi({\mathbf{X}}),$
(49)

or its time-dependent version,

$H_{\rm
BO}({\mathbf{X}},{\mathbf{P}})\,\psi({\mathbf{X}},t)=i\hbar\frac{\partial\psi({\mathbf{X}},t)}{\partial
t},$ (50)

where $\psi=\psi_{n}$ is the one, nonvanishing component of $\psi_{k}$, and
where $H_{\rm BO}=H_{nn}=H_{0,nn}+\hbar\,H_{1,nn}+\hbar^{2}\,H_{2,nn}$ is the
diagonal element of the Hamiltonian matrix, the effective Hamiltonian operator
on this subspace.

To write out $H_{\rm BO}$ explicitly we note first that $H_{1,nn}=0$ because
of (9). As for $H_{2,nn}$, we take the divergence of
${\mathbf{F}}_{\alpha;nn}(x)=\langle x;n|\nabla_{\alpha}|x;n\rangle=0$ to
obtain the identity,

$\nabla_{\alpha}\cdot{\mathbf{F}}_{\alpha;nn}(x)=(\nabla_{\alpha}\langle
x;n|)\cdot(\nabla_{\alpha}|x;n\rangle)+\langle
x;n|\nabla_{\alpha}^{2}|x;n\rangle=0.$ (51)

Then using (7) we have

$\displaystyle\sum_{p}{\mathbf{F}}_{\alpha;np}(x)\cdot{\mathbf{F}}_{\alpha;pn}(x)$
$\displaystyle=$ $\displaystyle-\sum_{p}(\nabla_{\alpha}\langle
x;n|)|x;p\rangle\cdot\langle x;p|(\nabla_{\alpha}|x;n\rangle)$ (52)
$\displaystyle=$ $\displaystyle-(\nabla_{\alpha}\langle
x;n|)\cdot(\nabla_{\alpha}|x;n\rangle)=\langle
x;n|\nabla_{\alpha}^{2}|x;n\rangle,$

where in the first step we use (7) and in the last step, (51). Altogether, we
obtain

$H_{\rm
BO}=H_{nn}=\sum_{\alpha=1}^{N}\frac{{\mathbf{P}}_{\alpha}^{2}}{2M_{\alpha}}+\epsilon_{n}(x)-\sum_{\alpha=1}^{N}\frac{1}{2M_{\alpha}}\langle
x;n|\nabla_{\alpha}^{2}|x;n\rangle.$ (53)

This is a standard Hamiltonian in Born-Oppenheimer theory. It is the same as
the projection of the original molecular Hamiltonian (1) onto the subspace
${\mathcal{S}}(x)$, that is,

$P(x)HP(x)\longleftrightarrow\delta_{kn}\,\delta_{nl}\,H_{\rm BO}.$ (54)

The final term in (53) is the “diagonal Born-Oppenheimer correction” (Valeev
and Sherrill (2003); Gauss  _et al._ (2006); Meek and Levine (2016); Gherib
_et al._ (2016); Culpitt  _et al._ (2022)). It is otherwise $H_{2,nn}$. It
appears to be a contribution to the potential energy, but physically it is a
part of the kinetic energy. Away from conical intersections, it is small,
being of order $\kappa^{4}$, and many authors take liberties with it (throwing
it away, or throwing part of it away, etc.). In the adiabatic basis it
diverges at conical intersections, for the same reason as the derivative
couplings (Matsunaga and Yarkony (1998); Yarkony (1997b)), but in our
treatment of single-surface problems we shall work in regions far enough from
conical intersections that this term is small. See Sec. III.6 for more
details. Near conical intersections we must use a multi-surface approach (see
subsection IV.1).

The terms in the molecular Hamiltonian that we have thrown away to obtain the
Born-Oppenheimer approximation are of order $\kappa^{2}$ or higher and thus
are small compared to the diagonal element $H_{0,nn}$ (see (47)). On the other
hand, they are large in comparison to $H_{2,nn}$, the diagonal Born-
Oppenheimer correction. Therefore there is a question of whether $H_{2,nn}$ is
even the correct, second order contribution to the Hamiltonian when the wave
function has a single component $\psi_{n}$. In fact, a perturbative treatment
of the off-block-diagonal terms, one that transforms them away rather than
throwing them away, generates a contribution at second order, what we call
$K_{22}$, that is of the same order in $\kappa$ as $H_{2,nn}$. See Sec.
III.10.

###  III.3 Dressing Transformations

As an alternative to the Born-Oppenheimer approximation, that is, just
throwing away terms in the molecular Hamiltonian (21), we will transform them
away by means of a series of unitary transformations. Since the terms to be
eliminated are small, each unitary transformation is close to the identity and
can be conveniently represented as the exponential of another operator, the
generator.

In the first step we transform the Hamiltonian by the unitary operator
$U_{1}=\exp(\hbar G_{1})$, where for convenience we split off a factor of
$\hbar$ from the anti-Hermitian generator $G_{1}$. The transformed Hamiltonian
$H^{\prime}$ can be written as an exponential series of iterated commutators,

$H^{\prime}=U_{1}HU_{1}^{\dagger}=e^{\hbar G_{1}}\,H\,e^{-\hbar
G_{1}}=H+\hbar[G_{1},H]+\frac{\hbar^{2}}{2!}[G_{1},[G_{1},H]]+\ldots$ (55)

Although we are writing $H$, $G_{1}$ and $U_{1}$ without matrix indices, we
will work in the Born-Oppenheimer representation and all such symbols should
be understood to be matrices of nuclear operators, for example, $H_{kl}$,
$G_{1,kl}$, etc. Likewise, commutators of operators are interpreted as matrix
commutators, for example,

$[G_{1},H]_{kl}=\sum_{p}\bigl{(}G_{1,kp}\,H_{pl}-H_{kp}\,G_{1,pl}\bigr{)}.$
(56)

The products that appear in (56) are operator products, as in (16), and the
ordering must be respected.

The strategy will be to choose $G_{1}$ so as to eliminate the off-block-
diagonal terms in $H$ to order $\kappa^{2}$. As it turns out, this means that
$G_{1}$ itself is of order $\kappa^{2}$. This leaves behind off-block-diagonal
terms in $H^{\prime}$ that are of order $\kappa^{4}$. To eliminate these, we
perform a second unitary transformation,
$H^{\prime\prime}=U_{2}\,H^{\prime}\,U_{2}^{\dagger}$, where
$U_{2}=\exp(\hbar^{2}G_{2})$. Operator $G_{2}$ is a second generator that
turns out to be of order $\kappa^{4}$. Proceeding this way we envision an
infinite sequence of unitary operators,

$U=\ldots U_{3}U_{2}U_{1},$ (57)

and a transformed Hamiltonian ${\bar{H}}$ that is formally free of off-block-
diagonal terms to all orders,

${\bar{H}}=UHU^{\dagger}=\ldots
U_{3}U_{2}U_{1}\,H\,U_{1}^{\dagger}U_{2}^{\dagger}U_{3}^{\dagger}\ldots$ (58)

These unitary transformations create a sequence of new representations for the
state of the system, what we will call “dressed” representations. Thus we can
distinguish between the “original” (that is, undressed) Born-Oppenheimer
representation, and successive dressed versions of it. See Table 1. If the
wave function is $\psi$ in the original Born-Oppenheimer representation (where
$\psi$ is short for the infinite-dimensional vector $\psi_{k}$), then the wave
function in the final, dressed representation is ${\bar{\psi}}=U\psi$ (that
is, it is the vector ${\bar{\psi}}_{k}=\sum_{l}U_{kl}\psi_{l}$). Then the
Schrödinger equation in the final, dressed representation is
${\bar{H}}{\bar{\psi}}=E{\bar{\psi}}$ or its time-dependent version, where
${\bar{H}}$ has no off-block-diagonal terms. Therefore there are solutions
${\bar{\psi}}$ that are nonzero only in the one slot $k=n$, and the effective
Hamiltonian for this component is the diagonal element ${\bar{H}}_{nn}$. This
is the (scalar) Hamiltonian driving this one component, which takes the place
of $H_{\rm BO}$ within dressing theory. We now explain the dressing
transformations in detail.

###  III.4 First Order Transformation and Weyl Symbols

In this subsection we begin by assuming that the working basis is the
adiabatic basis, that is, $|x;k\rangle=|ax;k\rangle$ for all $k$ (not just
$k=n$). This assumption causes (14) to hold for all $k$, $l$. It simplifies
the derivation somewhat but leads to difficulties that we will explain when we
encounter them.

In this subsection and henceforth we must change notation in order to
distinguish between nuclear operators and their corresponding Weyl symbols.
See Appendix A for the definition of the Weyl symbol of an operator and for
the invertible mapping between operators and their corresponding Weyl symbols.
As in Appendix A, we will place hats on nuclear operators such as
${\hat{\mathbf{X}}}_{\alpha}$ or ${\hat{\mathbf{P}}}_{\alpha}$, while omitting
the hats on the corresponding classical quantities or Weyl symbols
${\mathbf{X}}_{\alpha}$ or ${\mathbf{P}}_{\alpha}$. The notation ${\hat{x}}$
or ${\hat{p}}$ of Appendix A is to be identified with ${\hat{\mathbf{X}}}$ or
${\hat{\mathbf{P}}}$ here. We use the Weyl symbol correspondence only for
nuclear operators, not electronic ones, so we omit the hats on purely
electronic operators such as the electronic Hamiltonian $H_{e}$.

Thus, we should go back and put hats on operators in equations such as (21) to
conform with the new notation. We will not do this, but we note that when the
(hatted, operator) versions of (21) are converted into Weyl symbols, we obtain

$\displaystyle H_{0,kl}$ $\displaystyle=$
$\displaystyle\left(\sum_{\alpha=1}^{N}\frac{{\mathbf{P}}_{\alpha}^{2}}{2M_{\alpha}}\right)\delta_{kl}+W_{kl}({\mathbf{X}}),$
(59a) $\displaystyle H_{1,kl}$ $\displaystyle=$
$\displaystyle-i\sum_{\alpha=1}^{N}\frac{1}{M_{\alpha}}{\mathbf{P}}_{\alpha}\cdot{\mathbf{F}}_{\alpha;kl}({\mathbf{X}}),$
(59b) $\displaystyle H_{2,kl}$ $\displaystyle=$
$\displaystyle-\sum_{\alpha=1}^{N}\sum_{p}\frac{1}{2M_{\alpha}}{\mathbf{F}}_{\alpha;kp}({\mathbf{X}})\cdot{\mathbf{F}}_{\alpha;pl}({\mathbf{X}}),$
(59c)

where now $W_{kl}(x)=\delta_{kl}\,\epsilon_{k}(x)$ on account of our
assumptions about the working basis. These look almost the same as (21) but
now they are symbol equations and note that the expression for $H_{1,kl}$ has
simplified. In fact, this is an example of (233), illustrating the convenient
manner in which Weyl symbols handle operator ordering. Weyl symbols are
$c$-numbers, and can be multiplied in any order. Equations (59) contain
exactly the same information as (21); no information is lost on transforming
to Weyl symbols.

Likewise, in the case of (55), to conform with the new notation we replace
$H$, $H^{\prime}$ and $G_{1}$ by ${\hat{H}}$, ${\hat{H}}^{\prime}$ and
${\hat{G}}_{1}$, respectively, indicating operators (actually, matrices of
operators). Then converting these operator equations into symbols, we have the
symbol equation,

$H^{\prime}=H+\hbar[G_{1},H]_{*}+\frac{\hbar^{2}}{2}[G_{1},[G_{1},H]_{*}]_{*}+\ldots,$
(60)

where the notation $[G_{1},H]_{*}$ is explained in Appendix A. Note in
particular the definition (234) of the $*$-product and the discussion in
subsection A.3 of commutators of matrices of operators, especially (241).

Let us now write out the expansion of the symbol of the molecular Hamiltonian,
$H=H_{0}+\hbar\,H_{1}+\hbar^{2}\,H_{2}$, where the individual terms are given
by (59), and substitute it into (60). We obtain

$H^{\prime}=H_{0}+\hbar\bigl{(}H_{1}+[G_{1},H_{0}]_{*}\bigr{)}+\hbar^{2}\left(H_{2}+[G_{1},H_{1}]_{*}+\frac{1}{2}[G_{1},[G_{1},H_{0}]_{*}]_{*}\right)+\ldots$
(61)

The commutators shown contain an implicit $\hbar$-ordering. To expand this out
it is convenient to define

$T_{1}=[G_{1},H_{0}]_{*}=T_{10}+\hbar\,T_{11}+\ldots,$ (62)

where the $1$-subscript on $T_{1}$ is a reminder that this term is associated
with generator $G_{1}$, and where the expansion follows (241). In particular,
we have

$\displaystyle T_{10,kl}$ $\displaystyle=$
$\displaystyle[G_{1},H_{0}]_{kl}=\sum_{p}(G_{1,kp}\,H_{0,pl}-H_{0,kp}\,G_{1,pl}),$
(63a) $\displaystyle T_{11,kl}$ $\displaystyle=$
$\displaystyle\frac{i}{2}\sum_{p}\bigl{(}\\{G_{1,kp},H_{0,pl}\\}-\\{H_{0,kp},G_{1,pl}\\}\bigr{)},$
(63b)

where $T_{10}=[G_{1},H_{0}]$ is the ordinary commutator of matrices of Weyl
symbols (without the $*$) and where the brackets in (63b) are Poisson
brackets. Then (61) becomes

$H^{\prime}=H_{0}+\hbar\left(H_{1}+T_{10}\right)+\hbar^{2}\left(H_{2}+[G_{1},H_{1}]+T_{11}+\frac{1}{2}[G_{1},T_{10}]\right)+\ldots$
(64)

where we have dropped the $*$-subscript in commutators at order $\hbar^{2}$
since the ordinary, matrix commutator is the leading term in $\hbar$ (see
subsection A.3) and the higher order corrections are beyond the order to which
we are working.

The Moyal expansion of the $*$-product (236) obviously generates a power
series in $\hbar$. As pointed out at the end of subsection II.4, however, each
power of $\hbar$ in the $\hbar$-ordering of the molecular Hamiltonian (21) is
coupled to one power of $\kappa^{2}$. That is, terms $H_{0}$, $H_{1}$ and
$H_{2}$ are of order $\kappa^{0}$, $\kappa^{2}$ and $\kappa^{4}$,
respectively. Likewise, in the Moyal expansion (238), which is nominally an
expansion in $\hbar$, each successive bracket involves one higher derivative
with respect to momentum (see (237)). In this work the operators of interest
are polynomials in the nuclear momentum ${\hat{\mathbf{P}}}_{\alpha}$, which
is of order $\kappa^{-2}$, and their Weyl symbols are corresponding
polynomials in ${\mathbf{P}}_{\alpha}$. Therefore each momentum derivative
kills one factor of ${\mathbf{P}}_{\alpha}$ and raises the order of $\kappa$
by 2\. In other words, the Moyal expansion couples every power of $\hbar$ with
a power of $\kappa^{2}$. The result is that every term of order $\hbar^{n}$
that we encounter in the Moyal expansion behaves in the large-amplitude
scaling as of order $\kappa^{2n}$. In the following when we say, “first
order,” we shall mean, order $\hbar$ or $\kappa^{2}$, while “second order”
means order $\hbar^{2}$ or $\kappa^{4}$, etc.

The fact that the $\hbar$ ordering is closely related to the $\kappa$ ordering
should not be surprising, in view of the close connection between Born-
Oppenheimer theory, especially in the large-amplitude ordering, with
semiclassical theory (for example, the fact that the nuclear de Broglie
wavelength is of order $\kappa^{2}$). But it is a semiclassical theory only in
the nuclear variables, not the electronic.

###  III.5 The Generator $G_{1}$

Returning now to (64) we wish to choose $G_{1}$ to kill the off-block-diagonal
elements of $H_{1}$, so that we will have
$H^{\prime}_{1,nk}=H^{\prime}_{1,kn}=0$ for $k\neq n$. Thus, we require

$[G_{1},H_{0}]_{nk}=-H_{1,nk},\qquad(k\neq n).$ (65)

Of the two terms in $H_{0}$ (see (59a)) the kinetic energy is a multiple of
the identity matrix which commutes with the matrix $G_{1}$, while for the
potential energy we use (14) since we are working (for the time being) in the
adiabatic basis. Then the commutator in (65) is easy to evaluate and we
obtain,

$G_{1,nk}(\epsilon_{k}-\epsilon_{n})=-H_{1,nk},\qquad(k\neq n),$ (66)

or,

$G_{1,nk}=-i\sum_{\alpha=1}^{N}\frac{1}{M_{\alpha}}\frac{{\mathbf{P}}_{\alpha}\cdot{\mathbf{F}}_{\alpha;nk}({\mathbf{X}})}{\epsilon_{n}({\mathbf{X}})-\epsilon_{k}({\mathbf{X}})},\qquad(k\neq
n).$ (67)

If we carry out the same procedure for the $(kn)$-block, we obtain

$G_{1,kn}=-i\sum_{\alpha=1}^{N}\frac{1}{M_{\alpha}}\frac{{\mathbf{P}}_{\alpha}\cdot{\mathbf{F}}_{\alpha;kn}({\mathbf{X}})}{\epsilon_{k}({\mathbf{X}})-\epsilon_{n}({\mathbf{X}})},\qquad(k\neq
n),$ (68)

so that $G_{1,kn}=-G_{1,nk}^{*}$, as required if ${\hat{G}}_{1}$ is to be
anti-Hermitian. These are symbol equations; if we want the operator
${\hat{G}}_{1}$ then we can convert back to operators. For example, we find

${\hat{G}}_{1,nk}=-i\sum_{\alpha=1}^{N}\frac{1}{2M_{\alpha}}\left[{\hat{\mathbf{P}}}_{\alpha}\cdot\frac{{\mathbf{F}}_{\alpha;nk}({\hat{\mathbf{X}}})}{\epsilon_{n}({\hat{\mathbf{X}}})-\epsilon_{k}({\hat{\mathbf{X}}})}+\frac{{\mathbf{F}}_{\alpha;nk}({\hat{\mathbf{X}}})}{\epsilon_{n}({\hat{\mathbf{X}}})-\epsilon_{k}({\hat{\mathbf{X}}})}\cdot{\hat{\mathbf{P}}}_{\alpha}\right],\qquad(k\neq
n).$ (69)

in which the ordering of products must be respected. It will be appreciated
that Weyl symbols are easier to work with than operators.

Equation (65) is all we require of $G_{1}$; since it does not determine the
on-block components of $G_{1}$, we set these to zero, that is, $G_{1,nn}=0$
and $G_{1,kl}=0$ for $k\neq n$ and $l\neq n$. The result is that $G_{1}$ is
purely off-block-diagonal.

###  III.6 Energy Denominators, the Excluded Region, and Landau-Zener

The solution (67) diverges where the energy denominator
$\epsilon_{n}(x)-\epsilon_{k}(x)$ goes to zero. Since $k\neq n$ this can only
happen if the energy denominator goes to zero for the nearest neighbor,
$k=n\pm 1$. We define the “degeneracy manifold” as the subset of nuclear
configuration space where $\epsilon_{n}(x)=\epsilon_{k}(x)$ for $k=n\pm 1$,
that is, where level $n$ is degenerate with one of its immediate neighbors.
(If $n$ is the ground state, then it has only one neighbor, the first excited
state.) Degeneracy manifolds are usually called “seams.” As $x$ approaches the
degeneracy manifold, $G_{1}$ diverges and the perturbation expansion (55)
breaks down. Actually, not only does $G_{1}$ diverge, but so do the derivative
couplings and therefore $H_{1}$ itself. This is shown by the Feynman-Hellman
formula,

${\mathbf{F}}_{\alpha;nk}(x)=-\frac{\langle
ax;n|\nabla_{\alpha}H_{e}(x)|ax;k\rangle}{\epsilon_{n}(x)-\epsilon_{k}(x)},\qquad(k\neq
n),$ (70)

which applies in the adiabatic basis, as shown. Substituting this into (67),
we see that $G_{1}$ diverges as $(\epsilon_{n}-\epsilon_{k})^{-2}$ as we
approach the degeneracy manifold.

If the nearest neighbor energy denominator $\epsilon_{n}(x)-\epsilon_{k}(x)$
in (67) or in (70) is of order $\kappa^{0}=1$, then the size of the generators
$G_{1,nk}$ is determined by the factor ${\mathbf{P}}_{\alpha}/M_{\alpha}$, the
nuclear velocity, which is of order $\kappa^{2}$. In this case the generators
$G_{1,nk}$ (all of them, for all $k$) are of order $\kappa^{2}$, and therefore
small. But the nearest neighbor $\epsilon_{n}(x)-\epsilon_{k}(x)$ is of order
unity if the distance from the degeneracy manifold, call it $\Delta x$, is
also of order unity. (This is true in the sense of the $\kappa$-scaling.) At
such distances, the generators $G_{1,nk}$ are small and the first few terms of
the series (55) are rapidly decreasing. This is all we require of the series
for practical purposes; we do not require convergence. Thus, Moyal
perturbation theory, applied to a single-surface problem for the adiabatic
separation of the privileged surface $n$ from the others, is valid when
$\Delta x$ is of order $\kappa^{0}=1$.

It is also valid for smaller $\Delta x$, but not too much. As we approach the
degeneracy manifold, the nearest-neighbor $G_{1,nk}$ grows and when it becomes
of order unity the series (55) breaks down, that is, the terms are no longer
rapidly decreasing. Since $G_{1,nk}$ goes as
$(\epsilon_{n}-\epsilon_{k})^{-2}$, the breakdown occurs when the energy
difference is of order $\kappa$, which implies that $\Delta x$ is also of
order $\kappa$. This defines, at least in terms of its $\kappa$-scaling, the
region that must be excluded for the validity of single-surface methods.
Inside the excluded region it will be necessary to use multiple-surface
methods, which are discussed in subsection IV.1.

Since the nuclear de Broglie wave length is of order $\kappa^{2}$, the
excluded region, which is within a distance of order $\kappa$ from the
degeneracy manifold, contains on the order of $\kappa^{-1}$ de Broglie wave
lengths. The actual number in a practical situation depends on the magnitudes
of various matrix elements, the actual value of the nuclear momentum and other
quantities. But making a first stab based on the $\kappa$-scaling alone, we
estimate that the excluded region extends roughly 0.1 atomic unit away from
the degeneracy manifold, and that it contains roughly 10 nuclear de Broglie
wavelengths in this distance. This is all in the large-amplitude scaling.

In the case of small molecules, it is possible to visualize the excluded
region geometrically, but it is necessary to distinguish between the nuclear
configuration space and “shape space.” Nuclear configuration space is the
space $\mathbb{R}^{3N}$ upon which the nuclear positions in the lab frame,
${\mathbf{X}}_{\alpha}$, are coordinates. Nuclear configuration space contains
information about the center-of-mass position and the orientation of the
molecule. Shape space, on the other hand, is the $(3N-6)$-dimensional space in
which the center-of-mass position and orientation of the molecule are ignored.
Coordinates on shape space are translationally and rotationally invariant
quantities such as bond lengths and angles. In the case of triatomics ($N=3$),
shape space is topologically one half of $\mathbb{R}^{3}$ (Littlejohn and
Reinsch (1995)), while the codimension 2 degeneracy manifold is a curve in
this space (Yarkony (1996, 2001, 2004b, 2004a)). The excluded region is
therefore a tube whose width is of order $\kappa$ surrounding this curve. This
curve may bifurcate (Gordon, Glezaku, and Yarkony (1998); Yarkony (1997a)),
and, where it does, the excluded region surrounds it like the sleeves of a
shirt. The geometry of shape space and its relation to nuclear configuration
space is explained in detail by Littlejohn and Reinsch (1997).

The size of the excluded region depends on the nuclear momentum, and grows in
proportion to the square root of the latter. If we change the scaling of the
nuclear momentum, we find other scalings for the size of the excluded region.
In the small-amplitude scaling ${\mathbf{P}}_{\alpha}$ is of order
$\kappa^{-1}$, as explained above. This means that the size of the excluded
region is of order $\kappa^{3/2}$, or roughly 0.03 atomic units. Conversely,
if we make ${\mathbf{P}}_{\alpha}$ larger than the large-amplitude value of
$\kappa^{-2}$, we can make the excluded region larger. If we make
${\mathbf{P}}_{\alpha}$ large enough, we can make the size of the excluded
region as large as order $\kappa^{0}=1$, which means that the excluded region
is everywhere and the perturbation expansion fails everywhere. This occurs
when ${\mathbf{P}}_{\alpha}$ is of order $\kappa^{-4}$, so that the nuclear
velocity ${\mathbf{P}}_{\alpha}/M_{\alpha}$ is of order $\kappa^{0}=1$. This
means that the nuclear velocity is of the same order as the electron velocity,
at which point naturally there is no adiabatic separation of time scales. Such
momenta are too large to be of much interest for the usual chemical processes,
but it is nice to see that the theory breaks down completely where it should.

The generator $G_{1}$ has the same structure as the one-dimensional, Landau-
Zener transition probability (Landau (1932); Zener (1932); Stückelberg (1932);
Nikitin and Umanskii (1984)), which can be written as

$\exp\left(-\frac{2\pi V_{12}^{2}}{\hbar
v_{0}|V^{\prime}_{11}-V^{\prime}_{22}|}\right).$ (71)

Apart from the notation and numerical constants, the dimensionless quantity in
the exponential has the same structure as the reciprocal of $\hbar\,G_{1,nk}$,
which is also dimensionless. To see this we identify $v_{0}$ in (71), which
refers to the velocity with which the one-dimensional transition point is
passed, with ${\mathbf{P}}_{\alpha}/M_{\alpha}$; we identify $V_{12}$ in (71),
which is otherwise one half the minimum energy difference between two levels
at the transition point, with $(\epsilon_{n}-\epsilon_{k})/2$; and we identify
$V^{\prime}_{11}-V^{\prime}_{22}$ in (71), which is the $x$-derivative of the
difference in the diagonal elements of the Hamiltonian matrix, with $\langle
x;n|\nabla_{\alpha}H_{e}|x;k\rangle$. Thus, the quantity
$\exp[-\pi/(2\hbar|G_{1,nk}|)]$, for nearest neighbor $k$, behaves something
like a local Landau-Zener transition probability, and the excluded region is
where this quantity is of order unity.

The excluded region defines what we may think of as the “size of a conical
intersection,” that is, the distance over which it exerts its influence in
preventing the adiabatic separation of potential energy surfaces.

###  III.7 Diabatic Basis on ${\mathcal{S}}^{\perp}(x)$

When nondegenerate, the adiabatic basis vector $|ax;k\rangle$ is only
determined by (3) to within a phase, and when degenerate the set of degenerate
vectors is only determined to within the choice of an orthonormal frame inside
the degenerate eigenspace. Phase and frame conventions can be chosen at each
and every point $x$ of nuclear configuration space, and some such choice is
implied in the notation $|ax;k\rangle$, including in the continuum. But Born-
Oppenheimer theory generates derivatives of the basis vectors, for example,
inside the derivative couplings, so we must ask whether the phase and frame
conventions for the basis $\\{|ax;k\rangle\\}$ can be made as smooth functions
of $x$ over some region of interest of the nuclear configuration space.

The answer is no, for topological reasons it is impossible, in general, to
make a smooth assignment of phase and frame conventions for the vectors of the
adiabatic basis, and singularities (places where the basis vectors are not
differentiable with respect to $x$) are inevitable. Such singularities occur
on surfaces of two different types. First, when level $k$ is degenerate with
level $k\pm 1$, the derivatives $\nabla_{\alpha}|ax;k\rangle$ are not defined.
The manifolds where this occurs are like the degeneracy manifolds introduced
earlier, except they apply to any $k$, not just $k=n$. In the electrostatic
model these manifolds have codimension 2. As we approach such a manifold along
a curve, the frame has a limit, but the limit depends on the direction of
approach so derivatives are not defined on the manifold itself. (The
singularity is like that of the orthonormal frame in 3-dimensional space in
spherical coordinates,
$({\hat{\mathbf{r}}},{\hat{\bm{\theta}}},{\hat{\bm{\phi}}})$, as we approach
the $z$-axis.) The second type of surface upon which the adiabatic basis is
singular consists of codimension-1 manifolds that emanate from the degeneracy
manifolds. These manifolds can be moved about by gauge transformations like a
branch cut in complex variable theory. The adiabatic basis vectors suffer a
change in sign when crossing such a manifold and so are not differentiable
there.

These manifolds (of both types) proliferate as $k$ increases, and present
serious conceptual complications when derivatives of the adiabatic basis
arise. In addition there is the problem of the continuum, where the concept of
a nearest neighbor breaks down. For these reasons we advocate replacing the
adiabatic basis on ${\mathcal{S}}^{\perp}(x)$ with a diabatic basis, which is
a choice of an orthonormal frame $\\{|dx;k\rangle,k\neq n\\}$, with $d$ for
“diabatic,” that is a smooth function of $x$. Diabatic bases are usually used
only in the subspace ${\mathcal{S}}(x)$, but here we are using them in
${\mathcal{S}}^{\perp}(x)$. The fact that a Hilbert space always possesses a
discrete basis means that the diabatic basis in ${\mathcal{S}}^{\perp}(x)$ can
be chosen to be discrete; there is no continuum. The fact that the diabatic
basis on ${\mathcal{S}}^{\perp}(x)$ can be chosen to be a smooth function of
$x$, at least in some region of nuclear configuration space, is discussed by
Littlejohn, Rawlinson, and Subotnik (2022). (For example, one can use parallel
transport to create such a basis.) In addition, one can choose a diabatic
basis that is invariant under time reversal, so that the derivative couplings
${\mathbf{F}}_{\alpha;kl}(x)$ are real.

Therefore we now change the definition of the working basis, defining

$|x;k\rangle=\begin{cases}|ax;n\rangle&k=n,\\\
|dx;k\rangle&\text{otherwise.}\end{cases}$ (72)

In the working basis, the electronic Hamiltonian has matrix elements,

$\langle x;k|H_{e}(x)|x;l\rangle=\begin{cases}\epsilon_{n}(x)&k=l=n,\\\
W_{kl}(x)&\text{$k\neq n$ and $l\neq n$},\\\ 0&\text{otherwise}.\end{cases}$
(73)

The derivative couplings are still defined by (7) but are now taken with
respect to the (new) working basis, and the expressions (21) for the
Hamiltonian are still valid.

The adiabatic basis is widely used in the literature on Born-Oppenheimer
theory, in which formal sums over an infinite number of basis states,
including the continuum, are of frequent occurrance. Our scruples over the use
of this basis have nothing to do with the practical impossibility of the
numerical evaluation of such sums, but rather with the singularities of the
adiabatic basis as $x$ is varied. These present such a confusing picture that
it is hard to know what one is really talking about when derivatives with
respect to $x$ are taken. Our purpose in switching to a discrete, diabatic
basis is to obtain a formalism in which we can be confident that the answers
we obtain are correct.

On the other hand, for a fixed value of $x$ the adiabatic basis states do form
a complete set, and there is no harm in expanding vectors or operators in
terms of them. Therefore for $k\neq n$ we will consider it permissible to use
the adiabatic basis vectors $|ax;k\rangle$ in cases where no derivatives of
the basis vectors appear and the results are independent of the phase and
frame conventions of those vectors, for example, in (76) below. In cases where
derivatives of the basis vectors do appear we shall require the diabatic basis
(or the working basis, as we are now defining it).

The fact that the adiabatic basis is complete follows from the self-
adjointness of the electronic Hamiltonian $H_{e}(x)$. It has been claimed that
the electronic Hamiltonian is not self-adjoint (Sutcliffe (2012)), but we
believe the arguments for this are not well founded and that the matter is
settled by Jecko (2014).

However, since we are using the one adiabatic basis vector $|ax;n\rangle$ on
the privileged subspace, we should address its singularites. We will assume
that we are working in the region where the Moyal perturbation theory is
valid, so that $\Delta x$ is larger than order $\kappa$ and so that we avoid
the degeneracy manifold and its singularities. We can also avoid the $-1$
jumps in the phase of the adiabatic basis if our region is not so large as to
encircle a conical intersection; this can be achieved if our region is simply
connected (Juanes-Marcos, Althorpe, and Wrede (2005); Althorpe (2006, 2012);
Littlejohn, Rawlinson, and Subotnik (2023)). In the following we will assume
that $|ax;n\rangle$ is free of singularities in the region under
consideration, so that its derivatives are well defined.

###  III.8 The Resolvent $R(\epsilon,x)$

Now that we have changed the working basis, we must return to the calculation
of the generator $G_{1}$, starting with (65), and revise it. The
$(nk)$-component of the commutator for $k\neq n$ is now

$\displaystyle[G_{1},H_{0}]_{nk}$ $\displaystyle=$
$\displaystyle\sum_{p}G_{1,np}\,H_{0,pk}-\sum_{p}H_{0,np}\,G_{1,pk}=\sum_{p\neq
n}G_{1,np}\,W_{pk}(x)-\epsilon_{n}(x)G_{1,nk}$ (74) $\displaystyle=$
$\displaystyle\sum_{p\neq
n}G_{1,np}\bigl{[}W_{pk}(x)-\epsilon_{n}(x)\,\delta_{pk}\bigr{]}=-H_{1,nk},$

where again the kinetic energy part of $H_{0}$ does not contribute and where
we use (73) for the potential energy. In addition, in the first sum only terms
$p\neq n$ contribute since $W_{pk}$ is block-diagonal, and in the second sum
only the term $p=n$ contributes, for the same reason.

We now introduce the matrix of the resolvent operator in the working basis,
$R_{kl}(\epsilon,x)$ for $k,l\neq n$, defined by

$\sum_{p\neq
n}\bigl{[}\epsilon\,\delta_{kp}-W_{kp}(x)\bigr{]}R_{pl}(\epsilon,x)=\delta_{kl},\qquad(k,l\neq
n),$ (75)

so that $R(\epsilon,x)$ is basically the inverse of $\epsilon-H_{e}(x)$ on the
space ${\mathcal{S}}^{\perp}(x)$. We also define $R_{kl}(\epsilon,x)=0$ when
$k=n$ or $l=n$, so that $R(\epsilon,x)=0$ on ${\mathcal{S}}(x)$ and the matrix
$R_{kl}(\epsilon,x)$ is block-diagonal. The operator $R(\epsilon,x)$ is given
explicitly in the adiabatic basis by

$R(\epsilon,x)=\sum_{k\neq n}\frac{|ax;k\rangle\langle
ax;k|}{\epsilon-\epsilon_{k}(x)};$ (76)

by $R_{kl}(\epsilon,x)$ we mean the matrix elements of this operator in the
working basis. This operator is defined as long as $\epsilon$ is not equal to
any of the eigenvalues of $H_{e}(x)$ on ${\mathcal{S}}^{\perp}(x)$. In
particular, $R(\epsilon_{n},x)$ is defined over the region in which we are
working, which avoids the degeneracy manifold.

Now we may solve (74) for $G_{1,nk}$,

$\displaystyle G_{1,nk}$ $\displaystyle=$ $\displaystyle\sum_{l\neq
n}H_{1,nl}\,R_{lk}(\epsilon_{n},x)=-i\sum_{l\neq
n}\sum_{\alpha=1}^{N}\frac{1}{M_{\alpha}}{\mathbf{P}}_{\alpha}\cdot\langle
x;n|\nabla_{\alpha}|x;l\rangle\langle x;l|R(\epsilon_{n},x)|x;k\rangle$ (77)
$\displaystyle=$
$\displaystyle+i\sum_{\alpha=1}^{N}\frac{1}{M_{\alpha}}{\mathbf{P}}_{\alpha}\cdot\bigl{(}\nabla_{\alpha}\langle
x;n|\bigr{)}R(\epsilon_{n},x)|x;k\rangle,\qquad(k\neq n).$

In the last step we have used (7) to transfer the derivative to $\langle
x;n|$, we have extended the $l$-sum to include the term $l=n$ which vanishes
anyway, and then removed the resolution of the identity.

The Feynman-Hellman formula (70) applies only in the adiabatic basis but we
can obtain an equivalent result in the diabatic basis. First we note that the
resolvent satisfies

$\displaystyle P(x)R(\epsilon,x)$ $\displaystyle=$ $\displaystyle
R(\epsilon,x)P(x)=0,$ (78a) $\displaystyle R(\epsilon,x)$ $\displaystyle=$
$\displaystyle Q(x)R(\epsilon,x)=R(\epsilon,x)Q(x)=Q(x)R(\epsilon,x)Q(x),$
(78b)

as well as

$R(\epsilon,x)\bigl{[}\epsilon-H_{e}(x)\bigr{]}=\bigl{[}\epsilon-
H_{e}(x)\bigr{]}R(\epsilon,x)=Q(x).$ (79)

Next we write $[H_{e}(x)-\epsilon_{n}(x)]|x;n\rangle=0$ and apply
$\nabla_{\alpha}$ to obtain

$\bigl{[}\nabla_{\alpha}H_{e}(x)-\nabla_{\alpha}\epsilon_{n}(x)\bigr{]}|x;n\rangle+\bigl{[}H_{e}(x)-\epsilon_{n}(x)\bigr{]}\nabla_{\alpha}|x;n\rangle=0.$
(80)

We multiply this by $R(\epsilon_{n},x)$ to obtain

$R(\epsilon_{n},x)\,\nabla_{\alpha}H_{e}(x)|x;n\rangle=Q(x)\nabla_{\alpha}|x;n\rangle,$
(81)

where we use (78b) and where $R$ moves past the $c$-number
$\nabla_{\alpha}\epsilon_{n}(x)$ and annihilates $|x;n\rangle$. This can be
regarded as a result of first-order perturbation theory, in which $\delta
H_{e}=\sum_{\alpha}\bm{\xi}_{\alpha}\cdot\nabla_{\alpha}H_{e}(x)$, where
$\bm{\xi}_{\alpha}$ is a small displacement in nuclear configuration space. It
is often needed in Hermitian conjugated form,

$\bigl{(}\nabla_{\alpha}\langle x;n|\bigr{)}Q(x)=\langle
x;n|\nabla_{\alpha}H_{e}(x)\,R(\epsilon_{n},x).$ (82)

For example, the off-block-diagonal derivative couplings can now be written

$\displaystyle{\mathbf{F}}_{\alpha;nk}(x)$ $\displaystyle=$
$\displaystyle-\bigl{(}\nabla_{\alpha}\langle
x;n|\bigr{)}|x;k\rangle=-\bigl{(}\nabla_{\alpha}\langle
x;n|\bigr{)}Q(x)|x;k\rangle$ (83) $\displaystyle=$ $\displaystyle-\langle
x;n|\nabla_{\alpha}H_{e}(x)\,R(\epsilon_{n},x)|x;k\rangle,\qquad(k\neq n).$

This is the analog of the Feynman-Hellman formula (70) in the diabatic basis.
The other off-diagonal block is also useful,

${\mathbf{F}}_{\alpha;kn}(x)=-{\mathbf{F}}_{\alpha;nk}(x)^{*}=\langle
x;k|R(\epsilon_{n},x)\,\nabla_{\alpha}H_{e}(x)|x;n\rangle,\qquad(k\neq n).$
(84)

Similarly, we can replace $R$ by $QR$ in (77) and then substitute (82) to
obtain

$G_{1,nk}=i\sum_{\alpha=1}^{N}\frac{1}{M_{\alpha}}{\mathbf{P}}_{\alpha}\cdot\langle
x;n|\nabla_{\alpha}H_{e}(x)\,R(\epsilon_{n},x)^{2}|x;k\rangle,\qquad(k\neq
n).$ (85)

Equations (77) and (85) give us alternative versions of the $(nk)$-components,
for $k\neq n$, of $G_{1}$ in the diabatic basis. The derivation has nowhere
relied on dubious derivatives of the adiabatic basis vectors $|ax;k\rangle$
for $k\neq n$. The components $G_{1,kn}=-G^{*}_{1,nk}$ of the other off-
diagonal block are obtained by the anti-Hermiticity of $G_{1}$. The block-
diagonal components of $G_{1}$ vanish, as in the calculation using the
adiabatic basis, so that $G_{1}$ is purely off-block-diagonal.

The operator $R(\epsilon_{n},x)$, when acting on an arbitrary wave function,
has an effect whose order of magnitude is its largest eigenvalue, which by
(76) is the nearest neighbor value of $1/(\epsilon_{n}-\epsilon_{k})$. Thus
(85) shows that $G_{1,nk}$ in the diabatic basis behaves like the nearest
neighbor value of $(\epsilon_{n}-\epsilon_{k})^{-2}$, as in the adiabatic
basis.

###  III.9 The Second Order Hamiltonian

We have found the generator ${\hat{G}}_{1}$, whose Weyl symbol is given by
(77) or (85), that kills the off-block-diagonal components of ${\hat{H}}_{1}$.
To describe this situation it is convenient to say that an operator or its
symbol is “even” if it is block-diagonal, and “odd” if it is purely off-block-
diagonal. This terminology is borrowed from the Foldy-Wouthuysen
transformation (Foldy and Wouthuysen (1950); Bjorken and Drell (1964)), which
closely resembles the procedure we are following here in block-diagonalizing
the Born-Oppenheimer Hamiltonian. An arbitrary symbol or operator matrix is
neither even nor odd but can be uniquely decomposed into its even and odd
parts, which for example in the case of a symbol matrix we write as
$A=A^{e}+A^{o}$. When we multiply matrices of operators or their Weyl symbols,
the product of an even times an even or an odd times an odd matrix is even,
while the product of an odd times an even matrix or vice versa is odd. That
is, $(AB)^{e}=A^{e}B^{e}+A^{o}B^{o}$ and $(AB)^{o}=A^{e}B^{o}+A^{o}B^{e}$. The
same rule applies to commutators of matrices.

Concerning the operators or symbol matrices discussed so far, $H_{0}$ is
purely even, $H_{0}=H_{0}^{e}$, and $G_{1}$ is purely odd, $G_{1}=G_{1}^{o}$.
As for $H_{1}$ (see (59b)), it vanishes on the ${\mathcal{S}}$-block due to
the time-reversal invariance of the basis states, that is, $H_{1,nn}=0$, but
it has nonvanishing elements in the ${\mathcal{S}}^{\perp}$-block. Altogether,
it has both even and odd parts, $H_{1}=H_{1}^{e}+H_{1}^{o}$. Likewise we have
$H_{2}=H_{2}^{e}+H_{2}^{o}$. The symbol matrix
$T_{1}=[G_{1},H_{0}]_{*}=[G_{1}^{o},H_{0}^{e}]_{*}$ (see (62)) is purely odd,
$T_{1}=T_{1}^{o}$, as are the terms $T_{10}$, $T_{11}$, etc., of its
$\hbar$-expansion. The condition that ${\hat{G}}_{1}$ kill the off-block-
diagonal terms in the Hamiltonian through order $\kappa^{2}$ can be written as

$T_{10}=T_{10}^{o}=[G_{1},H_{0}]=-H_{1}^{o}.$ (86)

This makes $H^{\prime}_{1}$ (see (64)) purely even,

$H^{\prime}_{1}=H_{1}+T_{10}=H_{1}+[G_{1},H_{0}]=H_{1}^{e},$ (87)

which was the point of the first transformation.

As for $H^{\prime}$ at second order, we can use (64) and (86) to write

$H^{\prime}_{2}=H_{2}+[G_{1},H_{1}^{e}]+\frac{1}{2}[G_{1},H_{1}^{o}]+T_{11},$
(88)

which decomposes into even and odd parts,

$\displaystyle H^{\prime\,e}_{2}$ $\displaystyle=$ $\displaystyle
H_{2}^{e}+\frac{1}{2}[G_{1},H_{1}^{o}],$ (89a) $\displaystyle
H^{\prime\,o}_{2}$ $\displaystyle=$ $\displaystyle
H_{2}^{o}+[G_{1},H_{1}^{e}]+T_{11}.$ (89b)

The unitary transformation ${\hat{U}}_{1}=\exp(\hbar{\hat{G}}_{1})$ has killed
the off-block-diagonal elements of the Hamiltonian at first order, but not at
second order, as we see from (89b). To kill these we apply a second unitary
transformation, ${\hat{U}}_{2}=\exp(\hbar^{2}{\hat{G}}_{2})$, defining a
second transformed Hamiltonian,

${\hat{H}}^{\prime\prime}=e^{\hbar^{2}{\hat{G}}_{2}}\,{\hat{H}}^{\prime}\,e^{-\hbar^{2}{\hat{G}}_{2}}={\hat{H}}^{\prime}+\hbar^{2}[{\hat{G}}_{2},{\hat{H}}^{\prime}]+\ldots$
(90)

Converting this operator equation to symbols and expanding in $\hbar$, we find
$H^{\prime\prime}_{0}=H^{\prime}_{0}$, $H^{\prime\prime}_{1}=H^{\prime}_{1}$
and

$H^{\prime\prime}_{2}=H^{\prime}_{2}+[G_{2},H_{0}].$ (91)

Thus, to kill the off-block-diagonal terms in the Hamiltonian to second order
we require

$[G_{2},H_{0}]=-H^{\prime\,o}_{2},$ (92)

a condition on $G_{2}$ that may be compared to condition (86) on $G_{1}$.
Combined with (91), (92) shows that $H^{\prime\prime}_{2}=H^{\prime\,e}_{2}$,
which is given by (89a).

Thus we do not need $G_{2}$ in order to get $H^{\prime\prime}_{2}$, the
transformed Hamiltonian at second order, which is the more interesting object.
Note also that the higher order transformations (${\hat{U}}_{3}$, etc.) will
not change the Hamiltonian through second order, so that ${\bar{H}}$ and
$H^{\prime\prime}$ agree through second order (see (58)). Therefore we will
express the answers in terms of ${\bar{H}}$, not $H^{\prime\prime}$. Finally,
we are only interested in the $(nn)$-component of ${\bar{H}}$, which is the
effective (scalar) Hamiltonian driving the wave function on surface $n$.

We will call the Weyl symbol of this scalar Hamiltonian $K$, so that
$K={\bar{H}}_{nn}$. (The operator itself is ${\hat{K}}$.) Then we have
$K_{0}=H_{0,nn}$, which is given by (47), and $K_{1}=H_{1,nn}=0$. As for
$K_{2}={\bar{H}}_{2,nn}$, it is the $(nn)$-component of (89a), which has two
terms. We write $K_{2}=K_{21}+K_{22}$ for these, of which the first term,
$K_{21}=H_{2,nn}$, is the diagonal Born-Oppenheimer correction seen in (53).
In fact, we have $K=H_{BO}+K_{22}$ to second order, so that $K_{22}$ is
another correction, in addition to the diagonal Born-Oppenheimer correction.
The new correction $K_{22}$ is

$K_{22}=\frac{1}{2}[G_{1},H_{1}^{o}]_{nn}=\frac{1}{2}\sum_{k\neq
n}G_{1,nk}\,H_{1,kn}-\frac{1}{2}\sum_{k\neq n}H_{1,nk}\,G_{1,kn}.$ (93)

The sums are only over $k\neq n$ because $G_{1,nn}=H_{1,nn}=0$.

We work on the first sum. The quantity $G_{1,nk}$ is given by (85), while
$H_{1,kn}$ is given by (59b) with $l\to n$, into which we substitute (84).
This gives

$H_{1,kn}=-i\sum_{\beta}\frac{1}{M_{\beta}}{\mathbf{P}}_{\beta}\cdot\langle
x;k|R(\epsilon_{n},x)\nabla_{\beta}H_{e}(x)|x;n\rangle.$ (94)

With this expression for $H_{1,kn}$ the first sum in (93) can be extended to
all $k$ since the term $k=n$ vanishes anyway, after which we have a resolution
of the identity. In this way we evaluate the first sum in (93). When we
evaluate the second sum, including the minus sign, we find that it is equal to
the first sum.

Overall, we obtain

$\displaystyle K_{22}$ $\displaystyle=$
$\displaystyle\sum_{\alpha,\beta}\frac{1}{M_{\alpha}M_{\beta}}{\mathbf{P}}_{\alpha}\cdot\langle
x;n|\nabla_{\alpha}H_{e}(x)\,R(\epsilon_{n},x)^{3}\,\nabla_{\beta}H_{e}(x)|x;n\rangle\cdot{\mathbf{P}}_{\beta}$
(95) $\displaystyle=$
$\displaystyle\sum_{\alpha,\beta}\frac{1}{M_{\alpha}M_{\beta}}{\mathbf{P}}_{\alpha}\cdot\bigl{(}\nabla_{\alpha}\langle
x;n|\bigr{)}R(\epsilon_{n},x)\bigl{(}\nabla_{\beta}|x;n\rangle\bigr{)}\cdot{\mathbf{P}}_{\beta},$

where in the second form we have used (78b) and (81). By inserting resolutions
of the identity in the working basis this can be written

$K_{22}=-\sum_{\alpha,\beta}\frac{1}{M_{\alpha}M_{\beta}}\sum_{kl}[{\mathbf{P}}_{\alpha}\cdot{\mathbf{F}}_{\alpha;nk}(x)]R_{kl}(\epsilon_{n},x)[{\mathbf{F}}_{\beta;ln}(x)\cdot{\mathbf{P}}_{\beta}],$
(96)

where the derivative couplings are evaluated in the working basis. Finally,
transforming this to the adiabatic basis we obtain

$K_{22}=-\sum_{\alpha,\beta}\frac{1}{M_{\alpha}M_{\beta}}\sum_{k\neq
n}\frac{[{\mathbf{P}}_{\alpha}\cdot{\mathbf{F}}_{\alpha;nk}(x)][{\mathbf{F}}_{\beta;kn}(x)\cdot{\mathbf{P}}_{\beta}]}{\epsilon_{n}(x)-\epsilon_{k}(x)},$
(97)

where now the derivative couplings are evaluated in the adiabatic basis. This
is the result we would have obtained if we had worked strictly in the
adiabatic basis, without worrying about its singularities.

###  III.10 Discussion of $K_{22}$

Equations (95)–(97) give the symbol of the operator ${\hat{K}}_{22}$. The
operator itself can be obtained by replacing ${\mathbf{X}}_{\alpha}$ and
${\mathbf{P}}_{\alpha}$ by ${\hat{\mathbf{X}}}_{\alpha}$ and
${\hat{\mathbf{P}}}_{\alpha}$, in the order shown, plus commutator terms which
are of higher order in $\hbar$ and can be neglected. The resulting operator is
Hermitian.

The term $K_{22}$ should be included in Born-Oppenheimer descriptions of
nuclear motion, but it has usually been neglected, in spite of the fact that
it is of the same order of magnitude as the diagonal Born-Oppenheimer
correction, which has received extensive attention. (At least, this is true in
the large-amplitude ordering.) Admittedly both terms are small, and in many
applications it may be reasonable to neglect one or both. Nevertheless,
$K_{22}$ grows with nuclear momentum while the diagonal Born-Oppenheimer
correction does not, and at sufficiently large nuclear momenta it will become
important.

This term has been discussed in depth by Mátyus (2018), who calls it the
“nonadiabatic mass correction.” The idea is that this term modifies the
expression for the kinetic energy of the nuclei by adding a contribution due
to the electrons that are carried along by them during nuclear motions. We
refer to Mátyus (2018) for more details, mentioning also Kutzelnigg (2007).
Here we add a few further comments about $K_{22}$.

One is that this term is more difficult to evaluate numerically than the other
terms in the Born-Oppenheimer Hamiltonian. Mátyus (2018) has given ways of
evaluating it, basically by inverting the operator $\epsilon_{n}-H_{e}$
numerically. This requires careful attention to the choice of the subspace on
which the inversion is carried out, and to numerical issues due to the fact
that $\epsilon_{n}-H_{e}$ has a nontrivial kernel.

Actually the resolvent operator is ubiquitous in higher order Born-Oppenheimer
calculations, appearing not only in $K_{22}$ but also in matrix elements of
the electronic momentum and angular momentum (Nafie (1983); Freedman and Nafie
(1986)). These are needed in studies of radiative processes such as
vibrational circular dichroism (Stephens and Lowe (1985); Buckingham, Fowler,
and Galwas (1987); Nafie (1992)), as well as in expressions for electronic
currents (Barth  _et al._ (2009); Okuyama and Takatsuka (2009)). A method for
evaluating such matrix elements has been implemented by Stephens (1985) and
similar ideas have been proposed in a more general context by Patchkovskii
(2012). Basically, one calls on the fact that matrix elements involving the
resolvent occur in perturbation theory. By concocting a perturbation of the
right form, one can make the energy shifts in perturbation theory or their
derivatives agree with the desired matrix elements.

Another remark is that including $K_{22}$, the entire $K={\bar{H}}_{nn}$, that
is, the Weyl symbol of the operator ${\hat{K}}$, is the eigenvalue of the
symbol matrix $H_{kl}$, computed through second order. Since $H_{kl}$ is
already diagonal to lowest order in $\hbar$ its eigenvalues to higher order
can be computed by means of perturbation theory, taking into account the off-
diagonal contributions from $H_{1}$ and $H_{2}$. At second order this gives
both $K_{21}$ and $K_{22}$. This is interesting in terms of the work of
Littlejohn and Flynn (1991b), who show that the Weyl symbol of the operator
driving a given mode or polarization of a system of coupled wave equations is
the eigenvalue of the symbol matrix of the wave operator, plus two
corrections. The first is a correction due to Berry’s phase, and the second,
sometimes called the “no-name” term (Emmrich and Weinstein (1996)), involves
the Poisson bracket of the polarization vectors. In the Born-Oppenheimer case,
the Berry’s phase term vanishes on a single surface, due to time-reversal
invariance, while the Poisson bracket or no-name term is of third order in
$\hbar$ and so does not appear. This gives a complete accounting through
second order of the Hamiltonian $K$ relevant for a single surface in Born-
Oppenheimer theory, including $K_{22}$.

###  III.11 Tensor Notation and the Generator $G_{2}$

We now introduce tensor notation for the kinetic energy metric, as explained
in Appendix D, which is convenient especially for second-order calculations.
Note that henceforth there will be an implied summation on repeated indices
$\mu$, $\nu$, etc (one lower, one upper), where $\mu,\nu=1,\ldots,3N$. Note
also that $P_{\mu}$ (with a lower or covariant index) is the usual momentum,
which is of order $\kappa^{-2}$, while $P^{\mu}$ (with an upper or
contravariant index) is the same as the velocity, which is of order
$\kappa^{2}$. That is, $P^{\mu}={\dot{X}}^{\mu}$. In the new notation Eqs.
(59) become

$\displaystyle H_{0,kl}$ $\displaystyle=$
$\displaystyle\left(\frac{1}{2}P^{\mu}P_{\mu}\right)\delta_{kl}+W_{kl}(x),$
(98a) $\displaystyle H_{1,kl}$ $\displaystyle=$
$\displaystyle-i\,P^{\mu}\,F_{\mu;kl}(x)=-i\,P^{\mu}\langle
x;k|\partial_{\mu}|x;l\rangle,$ (98b) $\displaystyle H_{2,kl}$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{p}F^{\mu}{}_{;kp}(x)\,F_{\mu;pl}(x)=\frac{1}{2}(\partial^{\mu}\langle
x;k|)(\partial_{\mu}|x;l\rangle),$ (98c)

where we write the derivative couplings as

$F_{\mu;kl}(x)=\langle x;k|\partial_{\mu}|x;l\rangle=-(\partial_{\mu}\langle
x;k|)|x;l\rangle,$ (99)

with

$F^{\mu}{}_{;kl}(x)=G^{\mu\nu}\,F_{\nu;kl}(x).$ (100)

Notice that $F_{\mu;kl}=O(1)$ while $F^{\mu}{}_{;kl}=O(\kappa^{4})$. The
metric $G^{\mu\nu}$ or $G_{\mu\nu}$ (see (256 and (257)) is not to be confused
with the generators $G_{1}$, $G_{2}$, etc.

In the new notation the Feynman-Hellman formula (81) becomes

$R(\epsilon_{n},x)\,[\partial_{\mu}H_{e}(x)]\,|x;n\rangle=Q(x)\partial_{\mu}|x;n\rangle,$
(101)

and the generator $G_{1}$ becomes

$G_{1,nk}=i\,P^{\mu}(\partial_{\mu}\langle
x;n|)\,R(\epsilon_{n},x)|x;k\rangle,\qquad(k\neq n),$ (102)

which is equivalent to (77), while the diagonal element
$K={\bar{H}}_{nn}=K_{0}+\hbar K_{1}+\hbar^{2}K_{2}+\ldots$ of the diagonalized
Hamiltonian is given by

$\displaystyle K_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{2}P^{\mu}P_{\mu}+\epsilon_{n}(x),$ (103a)
$\displaystyle K_{1}$ $\displaystyle=$ $\displaystyle 0,$ (103b)
$\displaystyle K_{2}$ $\displaystyle=$ $\displaystyle K_{21}+K_{22}.$ (103c)

Here $K_{21}=H_{2,nn}$ is the diagonal Born-Oppenheimer correction,

$K_{21}=\frac{1}{2}(\partial^{\mu}\langle
x;n|)(\partial_{\mu}|x;n\rangle)=-\frac{1}{2}\langle
x;n|\partial^{\mu}\partial_{\mu}|x;n\rangle,$ (104)

which is equivalent to the last term in (53), and $K_{22}$ is the extra
correction,

$K_{22}=P^{\mu}P^{\nu}(\partial_{\mu}\langle
x;n|)R(\epsilon_{n},x)(\partial_{\nu}|x;n\rangle),$ (105)

which is a version of (95).

The second generator $G_{2}$ is worked out in Appendix E. It is purely odd
(off-block-diagonal) and anti-Hermitian. Its $(nk)$-block for $k\neq n$ is
given by

$\displaystyle G_{2,nk}=$ $\displaystyle-$ $\displaystyle
P^{\mu}P^{\nu}\,\partial_{\nu}[(\partial_{\mu}\langle
x;n|)R(\epsilon_{n},x)]R(\epsilon_{n},x)|x;k\rangle$ (106) $\displaystyle+$
$\displaystyle[\partial^{\mu}\epsilon_{n}(x)]\,(\partial_{\mu}\langle
x;n|)\,R(\epsilon_{n},x)^{2}\,|x;k\rangle$ $\displaystyle+$
$\displaystyle\frac{1}{2}(\partial_{\mu}\langle
x;n|)\partial^{\mu}[R(\epsilon_{n},x)\,|x;k\rangle],\qquad(k\neq n).$

Notice that the final term contains derivatives of the basis vectors
$|x;k\rangle$ on ${\mathcal{S}}^{\perp}$. The other off-diagonal block is
given by $G_{2,kn}=-G_{2,nk}^{*}$.

###  III.12 Changes of Diabatic Basis in ${\mathcal{S}}^{\perp}$

On ${\mathcal{S}}^{\perp}$, the working basis vectors are the same as the
diabatic basis vectors, that is, $|x;k\rangle=|dx;k\rangle$ for $k\neq n$ (see
(72)). But the diabatic basis on ${\mathcal{S}}^{\perp}$ is highly arbitrary,
that is, we can change it according to

$|dx;k\rangle=\sum_{l\neq
n}|d^{\prime}x;l\rangle\,V^{\perp}_{lk}(x),\qquad(k\neq n),$ (107)

where $V^{\perp}_{lk}$ is a smooth function of $x$. Here $|dx;k\rangle$ is the
old diabatic basis on ${\mathcal{S}}^{\perp}$ and $|d^{\prime}x;k\rangle$ is
the new one. We restrict $V^{\perp}_{lk}$ to be unitary so that the new basis
is orthonormal, and restrict it further to be real and orthogonal so that the
new basis is invariant under time reversal (we are assuming that the old basis
is invariant under time reversal).

According to (10), under the change of basis (107) the wave function
transforms according to

$\psi^{\prime}_{k}(x)=\sum_{l\neq
n}V^{\perp}_{kl}(x)\,\psi_{l}(x),\qquad(k\neq n).$ (108)

The wave function $\psi_{k}$ is not invariant under the change of diabatic
basis on ${\mathcal{S}}^{\perp}$, but, as we shall say, it transforms as a
vector.

Physical results cannot depend on the choice of diabatic basis vectors on
${\mathcal{S}}^{\perp}$. Note that the second order Hamiltonians $K_{21}$ and
$K_{22}$, given by (104) and (105), make no reference to the diabatic basis on
${\mathcal{S}}^{\perp}$ and thus are invariant under a change of that basis.
The wave function $\psi_{k}$ is not invariant but transforms in a simple
manner (as a vector) under a change of basis on ${\mathcal{S}}^{\perp}$. This
makes it easy to construct invariants, for example, $\sum_{k\neq
n}|\psi_{k}(x)|^{2}$. The components $G_{1,nk}$ of the generator $G_{1}$ are
not invariant but do transform as a vector in the index $k$, as seen by (102).
This means that physical results will never involve individual components
$G_{1,nk}$, but rather sums over $k\neq n$ of $G_{1,nk}$ times something that
transforms as a vector in $k$, in order to make an invariant.

On the other hand, it is clear from (107) that derivatives of the diabatic
basis vectors in ${\mathcal{S}}^{\perp}$, things that look like
$\partial_{\mu}|x;k\rangle$ for $k\neq n$, have a more complicated
transformation law, one that involves derivatives of the transformation matrix
$V^{\perp}_{kl}(x)$. Derivatives of the basis vectors on
${\mathcal{S}}^{\perp}$ do not occur in the components $G_{1,nk}$ of the
generator $G_{1}$, but they do occur in $G_{2}$, as shown in Appendix E. Such
terms complicate the construction of quantities that are invariant under the
change of basis (107).

This complication arises because we are using ordinary derivatives
$\partial_{\mu}$ instead of a properly defined covariant derivative, which is
the ordinary derivative plus a correction term. The correction term is,
however, of higher order in $\kappa^{2}$ than the ordinary derivative, so the
use of covariant derivatives mixes up the expansion in powers of $\kappa^{2}$.
We will stick with ordinary derivatives, to keep the formalism as simple as
possible.

###  III.13 Transformation of the Wave Function

Following the discussion of Sec. III.3, we write
$\hat{\mbox{${\bar{H}}$}}{\bar{\psi}}=E{\bar{\psi}}$ or
$\hat{\mbox{${\bar{H}}$}}{\bar{\psi}}=i\hbar\,\partial{\bar{\psi}}/\partial t$
for the Schrödinger equation in the completely dressed representation (with
hats on operators), where ${\bar{\psi}}$ is an infinite-dimensional vector and
$\hat{\mbox{${\bar{H}}$}}$ an infinite-dimensional matrix. But
$\hat{\mbox{${\bar{H}}$}}$ is block-diagonal and the wave function
${\bar{\psi}}$ is nonzero only in the slot $k=n$, so we write
${\bar{\psi}}_{k}=\delta_{kn}\,f$ where $f$ is the single wave function on
surface $n$, satisfying ${\hat{K}}f=Ef$ or ${\hat{K}}f=i\hbar\,\partial
f/\partial t$, where ${\hat{K}}=\hat{\mbox{${\bar{H}}$}}_{nn}$. The symbol $K$
of operator ${\hat{K}}$ is given by (103). If we solve this equation, we
obtain $f$ and hence ${\bar{\psi}}$ in the fully dressed Born-Oppenheimer
representation.

The wave function in the original Born-Oppenheimer representation is given by
$\psi={\hat{U}}^{\dagger}{\bar{\psi}}$, that is,
$\psi_{k}=\sum_{l}({\hat{U}}^{\dagger})_{kl}{\bar{\psi}}_{l}$. All the
components $\psi_{k}$ are nonzero. The component $k=n$ is the largest and
approximately equal to $f$. Keeping this one and neglecting all the others is
equivalent to the usual Born-Oppenheimer approximation. Thus, finding the
other components, $\psi_{k}$ for $k\neq n$, as well as the corrections to
$\psi_{n}$, amounts to finding corrections to the Born-Oppenheimer
approximation. It is feasible to calculate these corrections by hand to second
order, as we shall now do.

Since ${\bar{\psi}}_{k}=\delta_{nk}\,f$, the transformation equations are

$\displaystyle\psi_{n}$ $\displaystyle=$
$\displaystyle({\hat{U}}^{\dagger})_{nn}\,f=({\hat{U}}_{nn})^{\dagger}f=({\hat{U}}^{e}_{nn})^{\dagger}f,$
(109a) $\displaystyle\psi_{k}$ $\displaystyle=$
$\displaystyle({\hat{U}}^{\dagger})_{kn}\,f=({\hat{U}}_{nk})^{\dagger}f=({\hat{U}}^{o}_{nk})^{\dagger}f,\qquad(k\neq
n),$ (109b)

where, as we see, only the even part of ${\hat{U}}$ contributes to $\psi_{n}$
while only the odd part contributes to $\psi_{k}$ for $k\neq n$. On the other
hand, we have

${\hat{U}}=\ldots{\hat{U}}_{2}{\hat{U}}_{1}=\ldots\exp(\hbar^{2}{\hat{G}}_{2})\exp(\hbar{\hat{G}}_{1})=1+\hbar{\hat{G}}_{1}+\hbar^{2}\left({\hat{G}}_{2}+\frac{1}{2}{\hat{G}}_{1}^{2}\right)+\ldots,$
(110)

so, taking into account the fact that ${\hat{G}}_{1}$ and ${\hat{G}}_{2}$ are
purely odd, we can write (109) as

$\displaystyle\psi_{n}$ $\displaystyle=$ $\displaystyle
f+\frac{\hbar^{2}}{2}[({\hat{G}}_{1}^{2})_{nn}]^{\dagger}f,$ (111a)
$\displaystyle\psi_{k}$ $\displaystyle=$
$\displaystyle\hbar({\hat{G}}_{1,nk})^{\dagger}f+\hbar^{2}({\hat{G}}_{2,nk})^{\dagger}f,\qquad(k\neq
n),$ (111b)

which are valid through second order. Note that the first correction to
$\psi_{n}=f$ occurs at second order.

We work first on the term involving ${\hat{G}}_{1}^{2}$. Since
${\hat{G}}_{1}^{2}$ is Hermitian (as a matrix of operators) its diagonal
elements are Hermitian (as scalar operators), and the $\dagger$ is unnecessary
in (111a). Converting operators into symbols, we have

$\displaystyle({\hat{G}}_{1}^{2})_{nn}$ $\displaystyle=$
$\displaystyle\sum_{k}{\hat{G}}_{1,nk}{\hat{G}}_{1,kn}\longleftrightarrow\sum_{k}G_{1,nk}*G_{1,kn}=-\sum_{k}|G_{1,nk}|^{2}+\ldots$
(112) $\displaystyle=$ $\displaystyle-P^{\mu}P^{\nu}\,(\partial_{\mu}\langle
x;n|)\,R(\epsilon_{n},x)^{2}(\partial_{\nu}|x;n\rangle)+\ldots,$

where we use $\longleftrightarrow$ as in Appendix A, where we have used
$G_{1,kn}=-G_{1,nk}^{*}$ and where the ellipsis indicates higher order terms
in the expansion of the $*$-product, which we neglect since we are working at
second order. In the final step we have used (102) and evaluated the sum. When
we convert the final expression, a Weyl symbol, back into an operator, we are
free to put the factors ${\hat{P}}^{\mu}$, ${\hat{P}}^{\nu}$ in any order,
since the differences among the choices are commutators that are of higher
order in $\hbar$. With this understanding, we can write (111a) as

$\psi_{n}=f-\frac{\hbar^{2}}{2}(\partial_{\mu}\langle
x;n|)\,R(\epsilon_{n},x)^{2}(\partial_{\nu}|x;n\rangle){\hat{P}}^{\mu}{\hat{P}}^{\nu}f.$
(113)

As for (111b), we will need

$\displaystyle({\hat{G}}_{1,nk})^{\dagger}$ $\displaystyle\longleftrightarrow$
$\displaystyle G_{1,nk}^{*}=-i\,P^{\mu}\langle
x;k|\,R(\epsilon_{n},x)\,(\partial_{\mu}|x;n\rangle)$ (114)
$\displaystyle\longleftrightarrow$ $\displaystyle-i\,\langle
x;k|\,R(\epsilon_{n},x)\,(\partial_{\mu}|x;n\rangle)\,{\hat{P}}^{\mu}-\frac{\hbar}{2}\partial^{\mu}[\langle
x;k|\,R(\epsilon_{n},x)\,(\partial_{\mu}|x;n\rangle)],$

where we use (102) and a version of (231b). Since $G_{1}$ occurs at first
order, we must keep the commutator term (the final term) in (114) to obtain
results valid to second order.

A similar treatment of $G_{2}$ (see (106)) gives rise to three terms, one of
which is

$\frac{1}{2}\\{\partial^{\mu}[\langle
x;k|\,R(\epsilon_{n},x)]\\}(\partial_{\mu}|x;n\rangle),$ (115)

the term that depends on the derivatives of $\langle x;k|$. But so does the
commutator term in (114), and when the two are combined the derivatives of
$\langle x;k|$ drop out. Thus (111b) becomes

$\displaystyle\psi_{k}$ $\displaystyle=$ $\displaystyle-i\hbar\,\langle
x;k|\,R(\epsilon_{n},x)\,(\partial_{\mu}|x;n\rangle)\,{\hat{P}}^{\mu}f+\hbar^{2}\\{-\langle
x;k|\,R(\epsilon_{n},x)\,\partial_{\nu}[R(\epsilon_{n},x)\,(\partial_{\mu}|x;n\rangle)]\,{\hat{P}}^{\mu}{\hat{P}}^{\nu}$
(116) $\displaystyle+$ $\displaystyle[\partial^{\mu}\epsilon_{n}(x)]\,\langle
x;k|\,R(\epsilon_{n},x)^{2}(\partial_{\mu}|x;n\rangle)-\frac{1}{2}\langle
x;k|\,R(\epsilon_{n},x)\,(\partial^{\mu}\partial_{\mu}|x;n\rangle)\\}f,\quad(k\neq
n).$

We see that through second order, $\psi_{k}$ transforms as a vector, as it
must.

We can also write out the expansion of the wave function
$\Psi({\mathbf{X}},{\mathbf{r}})$, the solution of the original molecular
Schrödinger equation, (22) or (23). We write

$\Psi({\mathbf{X}},{\mathbf{r}})=\Psi_{0}({\mathbf{X}},{\mathbf{r}})+\Psi_{1}({\mathbf{X}},{\mathbf{r}})+\ldots,$
(117)

where the subscripts 0, 1, etc., indicate the power of $\kappa^{2}$. Then
$\Psi_{0}({\mathbf{X}},{\mathbf{r}})=\phi_{n}({\mathbf{X}};{\mathbf{r}})\,f({\mathbf{X}})=\langle{\mathbf{r}}|x;n\rangle
f({\mathbf{X}})$, which is the usual Born-Oppenheimer approximation, and

$\displaystyle\Psi_{1}({\mathbf{X}},{\mathbf{r}})$ $\displaystyle=$
$\displaystyle-i\hbar\sum_{k\neq
n}\phi_{k}({\mathbf{X}};{\mathbf{r}})\,\langle
x;k|\,R(\epsilon_{n},x)(\partial_{\mu}|x;n\rangle)\,{\hat{P}}^{\mu}f$ (118)
$\displaystyle=$
$\displaystyle-i\hbar\langle{\mathbf{r}}|\,R(\epsilon_{n},x)(\partial_{\mu}|x;n\rangle)\,{\hat{P}}^{\mu}f$
$\displaystyle=$
$\displaystyle-i\hbar\sum_{\alpha=1}^{N}\frac{1}{M_{\alpha}}\langle{\mathbf{r}}|\,R(\epsilon_{n},x)(\nabla_{\alpha}|x;n\rangle)\cdot{\hat{\mathbf{P}}}_{\alpha}f$

is the first order correction, written in several different ways. For brevity
we omit the second order correction, but it is easily written out with the
help of (113) and (116).

The first order correction to the Born-Oppenheimer wave function was
apparently first derived by Nafie (1983), in the context of small vibrations,
and has been rederived many times since. Our results apply in the more general
context of large-amplitude motions, and our second order results (113) and
(116) seem to be new.

##  IV Generalizations and Applications

In this section we outline the generalization of Moyal perturbation theory to
multisurface problems and we discuss some applications including the
calculation of electronic currents.

###  IV.1 Multiple Surfaces

For multi-surface problems we define an index set $I$ of $N_{l}$ adjacent
electronic energy levels,

$I=\\{k_{0},k_{0}+1,\ldots,k_{0}+N_{l}-1\\},$ (119)

which are regarded as “strongly coupled,” at least somewhere in the region of
nuclear configuration space under consideration. These levels define a
privileged subspace ${\mathcal{S}}(x)$,

${\mathcal{S}}(x)=\mathop{\rm span}\\{|ax;k\rangle,k\in I\\},$ (120)

which replaces (43). As before, the complementary, orthogonal subspace is
${\mathcal{S}}^{\perp}(x)$. Most of the formulas of Sec. III, which was
devoted to single surface problems in which the privileged subspace was the
single level $k=n$, are generalized to the multi-surface case by replacing
$k=n$ and $k\neq n$ by $k\in I$ and $k\notin I$, respectively. For example,
the resolvent is still given by (76), except now summed over $k\notin I$, and
the projector onto ${\mathcal{S}}(x)$ becomes

$P(x)=\sum_{k\in I}|ax;k\rangle\langle ax;k|.$ (121)

Now we can distinguish between “internal degeneracies,” those that occur among
levels $k\in I$, and “external degeneracies,” those that occur between a level
$k\in I$ and a level $k\notin I$. Now the region of nuclear configuration
space must be restricted to exclude external degeneracies and a margin around
them, while internal degeneracies are allowed.

Ever since (72) we have been using a working basis which is the diabatic basis
on ${\mathcal{S}}^{\perp}$. Now we must introduce a diabatic basis also on
${\mathcal{S}}$, that is, a basis $|dx;k\rangle$ for $k\in I$, so that
derivatives of the basis states are defined. Now we define the working basis
by

$|x;k\rangle=\begin{cases}|dx;k\rangle&k\in I,\\\ |dx;k\rangle&k\notin
I.\end{cases}$ (122)

To be clear about the notation, there are two diabatic bases, one on
${\mathcal{S}}$ and one on ${\mathcal{S}}^{\perp}$. The first is a linear
combination of the adiabatic basis vectors $|ax;k\rangle$ for $k\in I$, and
the second is a linear combination of those for $k\notin I$. For the
transformation between the diabatic and adiabatic bases on ${\mathcal{S}}$ we
write

$|dx;k\rangle=\sum_{l\in I}|ax;l\rangle\,V_{lk},\qquad k\in I,$ (123)

where $V_{lk}$ is an $N_{l}\times N_{l}$ orthogonal matrix,

$V_{lk}=\langle ax;l|dx;k\rangle.$ (124)

Now “even” and “odd” are defined relative to the sets of indices, $k\in I$ and
$k\notin I$.

The derivation proceeds much as in the case of single surface. The generator
$G_{1}$ (that is, its Weyl symbol) is purely odd and satisfies
$[G_{1},H_{0}]=-H_{1}^{o}$. This can be solved for one off-diagonal block,

$G_{1,kp}=\sum_{l,m\in I}\sum_{q\notin
I}(V^{\dagger})_{kl}\,V_{lm}\,H_{1,mq}\,R_{qp}(\epsilon_{l},x),\qquad\hbox{$k\in
I$, $p\notin I$},$ (125)

while the other is determined by the anti-Hermiticity of $G_{1}$. The presence
of the matrix $V$ in this formula means that effectively we have carried out
the linear algebra of the solution for $G_{1}$ in the adiabatic basis, then
transformed the result to the working (that is, diabatic) basis. The
calculation would be formally simpler if we had defined the working basis to
be the adiabatic basis on ${\mathcal{S}}$, but we wish to avoid any
derivatives of these adiabatic basis vectors since they are not defined at
internal degeneracies (that is, on the seams of the single surface problem).
Our presentation of the multisurface problem is designed so that such
derivatives never occur (but derivatives of the diabatic basis states on
${\mathcal{S}}$ are allowed).

Equation (125) raises the question of whether the solution $G_{1,kp}$ is a
smooth function of $x$, since the matrix $V$ is not smooth. (Matrix $V$
connects the diabatic basis, which is smooth, with the adiabatic basis, which
is not.) To answer this question we use the logic explained more fully by
Littlejohn, Rawlinson, and Subotnik (2022), that the singularities of the
adiabatic basis are caused by the topological impossibility of making smooth
phase conventions for energy eigenstates in the neighborhood of a degeneracy,
and that correspondingly expressions that are independent of such phase
conventions are smooth. It is easily seen that $G_{1,kp}$ given by (125) is
independent of the phase conventions for the adiabatic basis, and so is a
smooth function of $x$. A more sophisticated argument can be based on a
contour integral (see Appendix A of Littlejohn, Rawlinson, and Subotnik
(2022)). Now substituting (98b) into (125) and carrying out some manipulations
we find

$G_{1,kp}=-iP^{\mu}\sum_{l\in I}(V^{\dagger})_{kl}\,\langle
ax;l|[\partial_{\mu}R(\epsilon_{l},x)]|x;p\rangle,\qquad\hbox{$k\in I$,
$p\notin I$}.$ (126)

Results equivalent to (125) have previously been reported by Mátyus and Teufel
(2019), who work with an electronic operator that in our notation could be
written as

$\sum_{k\in I}\sum_{p\notin I}|x;k\rangle\,G_{1,kp}\,\langle x;p|+\hbox{\rm
h.c.}$ (127)

This is part of their formalism of operator-valued Weyl symbols. The use of
this formalism does not obviate questions of smoothness, for example, the
projector onto individual energy eigenstates, which appears in their formulas,
is not smooth at internal degeneracies.

The Hamiltonian ${\hat{K}}$ on the privileged subspace is an $N_{l}\times
N_{l}$ matrix, that is, with indices $k,l\in I$, whose symbol matrix can be
written $K=K_{0}+\hbar K_{1}+\hbar^{2}K_{2}$, to second order. As for $K_{0}$,
it is a version of (98a),

$K_{0,kl}=\left(\frac{1}{2}P^{\mu}P_{\mu}\right)\delta_{kl}+W_{kl},\qquad
k,l\in I,$ (128)

where $W_{kl}$ is a full matrix because we are using the diabatic basis for
$k\in I$. As for $K_{1}$, it is the same as $H_{1}$ on the diagonal block
corresponding to ${\mathcal{S}}$,

$K_{1,kl}=-i\,P^{\mu}\langle x;k|\partial_{\mu}|x;l\rangle,\qquad k,l\in I.$
(129)

Note that in multi-surface problems, $K_{1}$ has non-vanishing, off-diagonal
components with derivative couplings. Finally, $K_{2}=K_{21}+K_{22}$ consists
of two terms, of which $K_{21}$ is a version of (98c),

$K_{21,kl}=\frac{1}{2}(\partial^{\mu}\langle
x;k|)(\partial_{\mu}|x;l\rangle),\qquad k,l\in I,$ (130)

which is a multi-surface version of the diagonal Born-Oppenheimer correction.
The second correction $K_{22}$ can be written

$K_{22,kl}=\frac{1}{2}P^{\mu}P^{\nu}\sum_{m,n\in I}(V^{\dagger})_{km}\,\langle
ax;m|[\partial_{\mu}P(x)][R(\epsilon_{m},x)+R(\epsilon_{n},x)][\partial_{\nu}P(x)]|ax;n\rangle\,V_{nl},\qquad
k,l\in I.$ (131)

In this expression we have taken derivatives of the projection operator
$P(x)$, which is smooth, but not of the adiabatic basis states, which are not.
(The projection operator is not to be confused with the momenta $P^{\mu}$,
$P^{\nu}$.) A result equivalent to (131) has been previously reported by
Mátyus and Teufel (2019).

The second generator $G_{2}$ and the corrections to the wave function may also
be computed.

###  IV.2 Electronic Momentum

In single-surface problems the usual (lowest order) Born-Oppenheimer
approximation gives a vanishing result for the matrix elements of electron
momenta (linear and angular) and in the calculation of electronic currents.
For such problems it is necessary to go to higher order. The same applies to
the single-surface matrix elements of any Hermitian electronic operator that
is odd under time-reversal. In this subsection we explain how the issue may be
treated within Moyal perturbation theory.

We begin with the momentum ${\mathbf{p}}_{i}$ of one of the electrons, which
is odd under time reversal, $T{\mathbf{p}}_{i}T^{\dagger}=-{\mathbf{p}}_{i}$.
This is a useful practice case, although it has a special feature not shared
with other such operators, namely, the operator ${\mathbf{p}}_{i}$ can be
represented as an exact time derivative. That is, let
$\Psi({\mathbf{X}},{\mathbf{r}},t)$ be the time-dependent wave function of the
molecule. Then

$\langle\Psi|{\mathbf{p}}_{i}|\Psi\rangle_{Xr}=\frac{m_{e}}{i\hbar}\langle\Psi|[{\mathbf{r}}_{i},H]|\Psi\rangle_{Xr}=m_{e}\frac{d}{dt}\langle\Psi|{\mathbf{r}}_{i}|\Psi\rangle_{Xr}.$
(132)

This is exact in the electrostatic model, but if the Born-Oppenheimer
approximation (48) is used in it (where $\Psi$ and $\psi$ are allowed to
depend on $t$) then the left-hand side vanishes while the right-hand side does
not. To reconcile the two it is necessary to go to higher order in the Born-
Oppenheimer expansion in the evaluation of the left-hand side. Note that
${\mathbf{r}}_{i}$ is even under time reversal, so its matrix element does not
vanish at lowest order.

To evaluate $\langle\Psi|{\mathbf{p}}_{i}|\Psi\rangle$ we use a version of
(18),

$\langle\Psi|A|\Psi\rangle_{Xr}=\sum_{kl}\langle\psi_{k}|{\hat{A}}_{kl}|\psi_{l}\rangle_{X}=\sum_{kl}\langle{\bar{\psi}}_{k}|\hat{\mbox{${\bar{A}}$}}_{kl}|{\bar{\psi}}_{l}\rangle_{X},$
(133)

showing how the expectation value of an operator $A$ changes as the
representation moves from the molecular to the original Born-Oppenheimer to
the dressed Born-Oppenheimer (that is, across the first row of Table 1). Here
we define $\hat{\mbox{${\bar{A}}$}}={\hat{U}}{\hat{A}}{\hat{U}}^{\dagger}$, as
in (58), using hats to indicate nuclear operators or matrices thereof. If the
original $A$ is a purely electronic operator then ${\hat{A}}_{kl}=\langle
x;k|A|x;l\rangle$ (see subsection II.3), which is a purely multiplicative
nuclear operator (that is, a function of ${\mathbf{X}}$ only). In this case
there is hardly a difference between the operator ${\hat{A}}_{kl}$ and its
Weyl symbol $A_{kl}$.

We identify the original operator $A$ with ${\mathbf{p}}_{i}$, and write
$({\mathbf{p}}_{i})$ for the matrix representing it in the original Born-
Oppenheimer representation, so that

$({\mathbf{p}}_{i})_{kl}=\langle x;k|{\mathbf{p}}_{i}|x;l\rangle.$ (134)

If we wish to emphasize that this is a nuclear operator, we will write this
matrix as $\widehat{({\mathbf{p}}_{i})}$ or its components as
$\widehat{({\mathbf{p}}_{i})}_{kl}$, but these are purely multiplicative