Effective non-adiabatic Hamiltonians for the quantum nuclear motion over coupled electronic states

TL;DR

This paper develops a systematic method to derive effective non-adiabatic Hamiltonians for nuclear motion in molecules, including higher-order corrections and mass-correction terms, improving the accuracy of quantum nuclear dynamics models.

Contribution

It introduces a general framework for block-diagonalizing the electron-nuclear Hamiltonian up to arbitrary order, deriving explicit formulas for second- and third-order corrections, and incorporating previously neglected outlying electronic states.

Findings

Derived explicit second- and third-order effective Hamiltonians.

Revealed generalized mass-correction terms in multi-dimensional coupled electronic states.

Unified treatment of adiabatic and non-adiabatic nuclear motion.

Abstract

The quantum mechanical motion of the atomic nuclei is considered over a single- or a multi-dimensional subspace of electronic states which is separated by a gap from the rest of the electronic spectrum over the relevant range of nuclear configurations. The electron-nucleus Hamiltonian is block-diagonalized up to $\mathcal{O}(\varepsilon^{n+1})$ through a unitary transformation of the electronic subspace and the corresponding $n$th-order effective Hamiltonian is derived for the quantum nuclear motion. Explicit but general formulae are given for the second- and the third-order corrections. As a special case, the second-order Hamiltonian corresponding to an isolated electronic state is recovered which contains the coordinate-dependent mass-correction terms in the nuclear kinetic energy operator. For a multi-dimensional, explicitly coupled electronic band, the second-order Hamiltonian contains the usual BO terms and non-adiabatic corrections but generalized mass-correction terms appear as well. These, earlier neglected terms, perturbatively account for the outlying (discrete and continuous) electronic states not included in the explicitly coupled electronic subspace.

Equations

$$({\hat{\mathbf{H}}}^{(2)}{)}_{\alpha \beta }$$

$${\mathbf{M}}_{\alpha \beta ,ij}=\sum _{a,b=1}^d⟨{\psi }_{\alpha }|P_a({\partial }_{j}P)({\mathcal{R}}_a+{\mathcal{R}}_b)({\partial }_{i}P)P_b|{\psi }_{\beta }⟩,$$

$$\hat{H}$$

$$-\sum _{i=1}^N\frac{1}{2m_i}{\Delta }_{{\tilde{R}}_i}$$

$$=-\frac{{\epsilon }^2}{2}\Delta =-\frac{{\epsilon }^2}{2}(\epsilon \nabla )(\epsilon \nabla )=-\frac{1}{2}\sum _{j=1}^{3N}(\epsilon {\partial }_j{)}^2.$$

$$\hat{H}=-\frac{{\epsilon }^2}{2}\Delta +H_{\text{e}}(R)=:\hat{K}+H_{\text{e}}(R),$$

$$[\hat{p},f(R)]\psi =[\hat{p},f(R)]\psi =\hat{p}(f\psi )-f\hat{p}\psi =(\hat{p}f)\psi =-\mathrm{i}\epsilon (\nabla f)(R)\psi$$

$$\hat{H}=\left(\begin{array}{cc}P\hat{H}P& P\hat{H}{P}^{\perp }\\ P^{\perp }\hat{H}P& P^{\perp }\hat{H}{P}^{\perp }\end{array}\right),$$

$$P^{\perp }\hat{H}P$$

$$=P^{\perp }\hat{K}P([P,H_{\text{e}}]=0)$$

$$=P^{\perp }[\hat{K},P]P(PP=P\text{ and }{P}^{\perp }P=0)$$

$$=-\frac{{\epsilon }^2}{2}{P}^{\perp }[\Delta ,P]P$$

$$=-\frac{1}{2}{P}^{\perp }({\epsilon }^2\Delta P-P{\epsilon }^2\Delta )P$$

$$=-\frac{1}{2}{P}^{\perp }(\epsilon \nabla \cdot \epsilon (\nabla P)+\epsilon (\nabla P)\cdot \epsilon \nabla )P$$

$$=-\frac{\epsilon }{2}{P}^{\perp }\left(\epsilon \nabla \cdot P^{\prime }+P^{\prime }\cdot \epsilon \nabla \right)P$$

$$=:-\epsilon P^{\perp }{\hat{P}}_{\nabla }^{\prime }P,$$

$${\hat{A}}_{\nabla }^{\prime }:=\frac{1}{2}(\epsilon \nabla \cdot A^{\prime }+A^{\prime }\cdot \epsilon \nabla )=\frac{\mathrm{i}}{2}(\hat{p}\cdot A^{\prime }+A^{\prime }\cdot \hat{p}),$$

$$P\hat{H}{P}^{\perp }=\epsilon P{\hat{P}}_{\nabla }^{\prime }{P}^{\perp }.$$

$$\hat{Q}$$

$$P^{\perp }\hat{H}P=\epsilon P^{\perp }\hat{Q}P\text{and}P\hat{H}{P}^{\perp }=\epsilon P\hat{Q}{P}^{\perp },$$

$$\hat{H}=\left(\begin{array}{cc}{\hat{H}}_P^{(0)}& \epsilon P\hat{Q}{P}^{\perp }\\ \epsilon P^{\perp }\hat{Q}P& P^{\perp }\hat{H}{P}^{\perp }\end{array}\right)=\left(\begin{array}{cc}{\hat{H}}_P^{(0)}& 0\\ 0& P^{\perp }\hat{H}{P}^{\perp }\end{array}\right)+\mathcal{O}(\epsilon ).$$

$$\hat{H}=\left(\begin{array}{cc}\hat{\Pi }\hat{H}\hat{\Pi }& \hat{\Pi }\hat{H}{\hat{\Pi }}^{\perp }\\ {\hat{\Pi }}^{\perp }\hat{H}\hat{\Pi }& {\hat{\Pi }}^{\perp }\hat{H}{\hat{\Pi }}^{\perp }\end{array}\right)=\left(\begin{array}{cc}{\hat{H}}_P^{(n)}& 0\\ 0& {\hat{\Pi }}^{\perp }\hat{H}{\hat{\Pi }}^{\perp }\end{array}\right)+\mathcal{O}({\epsilon }^{n+1})$$

$$\hat{\Pi }={\mathrm{e}}^{\mathrm{i}\epsilon \hat{S}}P{\mathrm{e}}^{-\mathrm{i}\epsilon \hat{S}},$$

$${\hat{\mathcal{H}}}_P^{(n)}:={\mathrm{e}}^{-\mathrm{i}\epsilon \hat{S}}{\hat{H}}_P^{(n)}{\mathrm{e}}^{\mathrm{i}\epsilon \hat{S}}$$

$$={\mathrm{e}}^{-\mathrm{i}\epsilon \hat{S}}\hat{\Pi }\hat{H}\hat{\Pi }{\mathrm{e}}^{\mathrm{i}\epsilon \hat{S}}$$

$$=P{\mathrm{e}}^{-\mathrm{i}\epsilon \hat{S}}\hat{H}{\mathrm{e}}^{\mathrm{i}\epsilon \hat{S}}P,$$

$${\hat{A}}^{\mathrm{D}}:=P\hat{A}P+P^{\perp }\hat{A}{P}^{\perp }\text{mbox}and{\hat{A}}^{\mathrm{OD}}:=P\hat{A}{P}^{\perp }+P^{\perp }\hat{A}P,$$

$$[{\hat{A}}^{\mathrm{D}},P]=[{\hat{A}}^{\mathrm{D}},P^{\perp }]=0\text{mbox}and[{\hat{A}}^{\mathrm{OD}},P]=P^{\perp }\hat{A}P-P\hat{A}{P}^{\perp },$$

$$[{\hat{A}}^{\mathrm{D}},{\hat{B}}^{\mathrm{D}}{]}^{\mathrm{OD}}=0,[{\hat{A}}^{\mathrm{D}},{\hat{B}}^{\mathrm{OD}}{]}^{\mathrm{D}}=0,\text{mbox}and[{\hat{A}}^{\mathrm{OD}},{\hat{B}}^{\mathrm{OD}}{]}^{\mathrm{OD}}=0.$$

$${\hat{K}}^{\mathrm{OD}}=\epsilon \hat{Q}.$$

$$\hat{A}\mapsto {\mathcal{L}}_{\hat{B}}(\hat{A}):=-\mathrm{i}[\hat{B},\hat{A}],$$

$${\mathcal{L}}_{H_{\mathrm{e}}}(\hat{A})=-\mathrm{i}[H_{\mathrm{e}},\hat{A}].$$

$${\mathcal{L}}_{H_{\mathrm{e}}}({\hat{A}}^{\mathrm{D}})={\mathcal{L}}_{H_{\mathrm{e}}}(\hat{A}{)}^{\mathrm{D}}\text{mbox}and{\mathcal{L}}_{H_{\mathrm{e}}}({\hat{A}}^{\mathrm{OD}})={\mathcal{L}}_{H_{\mathrm{e}}}(\hat{A}{)}^{\mathrm{OD}}.$$

$${\mathcal{I}}_{H_{\mathrm{e}}}(B)=\mathrm{i}\sum _{a=1}^d\left({\mathcal{R}}_{a}B{P}_a-P_{a}B{\mathcal{R}}_a\right),$$

$${\mathcal{R}}_a(R):=(H_{\mathrm{e}}(R)-E_a(R){)}^{-1}{P}^{\perp }(R)$$

Full text

I Introduction

Molecules are central paradigms of chemistry. They acquire unique features as physical objects due to the three orders of magnitude difference in the mass of their constituent particles, the electrons and the atomic nuclei.

Significant improvements in the energy resolution of spectroscopy experiments Beyer et al. (2018); Cheng et al. (2018), developed or adapted for the molecular domain, provide us with new pieces of information which can be deciphered, if a similarly precise and accurate theoretical description becomes available. For this purpose, it is necessary to re-consider the usual approximations used in quantum chemistry, in particular, the Born–Oppenheimer (BO) and the non-relativistic approximations. These two approximations give rise to ‘effects’ which play a role at the presently available experimental energy resolution. Further ‘effects’ may also be visible, for example, the interaction between the molecule and the quantized photon field are more and more appreciated as significant corrections to the molecular energy at this resolution Piszczatowski et al. (2009). Molecules have a large number of sharp spectral transitions, which can be measured experimentally to high precision. The interplay of the many small (or often not so small) effects (may) show up differently for the transitions between the different dynamical domains, so we cannot rely on the cancellation of the small effects, but their explicit computation Piszczatowski et al. (2009); Puchalski et al. (2017); Wang and Yan (2018), and hence further development of molecular quantum theory is necessary.

As to the coupling of the quantum mechanical motion of the electrons and the atomic nuclei benchmark energies and wave functions can be obtained by the explicit, variational solution of the few-particle Schrödinger equation Suzuki and Varga (1998); Cafiero et al. (2003); Bubin et al. (2009); Mátyus and Reiher (2012); Mátyus (2013); Bubin et al. (2013); Mitroy et al. (2013); Pachucki and Komasa (2018); Muolo et al. (2018a, b); Mátyus (2019). We call this direction pre-Born–Oppenheimer (pre-BO) theory, because it completely avoids the BO separation, nor does it evoke the concept of a potential energy surface. Obviously, a pre-BO computation captures ‘all’ non-adiabatic ‘effects’. Although all bound and low-lying resonance states of the three-particle H${}_{2}^{+}=\{\text{p}^{+},\text{p}^{+},\text{e}^{-}\}$ molecular ion have been recently reported to an outstanding precision Korobov (2018), already for four- and five-particle systems Muolo et al. (2018b) the explicit many-particle solution is typically limited to a few selected states due to the increased computational cost and other methodological challenges.

In order to compute (reasonably) accurate energies and wave functions over a broad dynamical range, we look for effective non-adiabatic Hamiltonians.

There is a vast literature about dominant non-adiabatic features Mead and Truhlar (1982); Pacher et al. (1988); Yarkony (1996), such as conical intersections, the geometric phase effect, or Jahn–Teller systems. Practical diabatization procedures Pacher et al. (1988); Viel and Eisfeld (2004); Karman et al. (2016) have been developed which make it possible to couple close-coming electronic states (to a good approximation) without the explicit knowledge of non-adiabatic coupling vectors, which are tedious to compute and burdensome to interpolate for larger systems. These effects are sometimes called first-order non-adiabatic effects and represent qualitatively important features for the molecular dynamics.

In the case of an isolated electronic state the dynamics is well described using a single potential energy surface. In order to obtain more accurate results, one would need to couple an increasing number of electronic states. These additional, explicitly coupled states would give small but non-negligible contributions to the molecular energy when studied under high resolution. Tightly converging the rovibrational (rovibronic) energy by increasing the number of explicitly coupled electronic states is inpractical (or impossible, since one would need to include also continuum electronic states).

A correction which is often computed is the diagonal Born–Oppenheimer correction (DBOC), which gives a mass-dependent contribution to the potential energy surface (PES). It has been (empirically) observed that in rovibrational computations carried out on a single potential energy surface, it is ‘better’ (in comparison with experiments) to use the atomic mass, especially for heavier atoms, instead of the nuclear mass, which would have been rigorously dictated by the BO approximation. The difference between the atomic and the nuclear mass is small, it is the mass of the electrons. This empirical adjustment of the mass of the nuclei used in the rovibrational kinetic energy operator has been supplemented with the argument that attaching the electrons mass to the nuclear mass approximately accounts for small, ‘secondary’ non-adiabatic effects Kutzelnigg (2007). The empirical adjustment is motivated by the picture that the electrons ‘follow’ the atomic nuclei in their motion.

For an isolated electronic state, the effective rovibrational Hamiltonian including the rigorous mass-correction terms has been derived and re-derived in a number of independent and different (perturbative) procedures Bunker and Moss (1977, 1980); Aubert-Frecon et al. (1994); Schwenke (2001a, b); Pachucki and Komasa (2009); Scherrer et al. (2017) over the past decades and were numerically computed for a few systems Bunker et al. (1977); Aubert-Frecon et al. (1994); Bak et al. (2005); Holka et al. (2011); Pachucki and Komasa (2012); Sauer et al. (2005); Przybytek et al. (2017); Mátyus (2018a, b). These mass-correction terms perturbatively account for the effect of all other electronic states on the rovibrational motion.

We would like to have similar perturbative corrections, for a system which is governed by not only a single but by a few close-lying electronic states, which are explicitly, non-adiabatically coupled, but which are distant (separated by a finite gap) from the rest of the (discrete and continuous) electronic states over the relevant range of the nuclear coordinates.

To the best of our knowledge, the explicit formulae have never been derived for an electronic band which includes multiple electronic states, but all the necessary ideas and techniques have been available in the literature, in particular, in relation with the space-adiabatic theory of quantum mechanics Martinez and Sordoni (2002); Teufel (2003); Panati et al. (2007); Martinez and Sordoni (2009); Weigert and Littlejohn (1993), but also other rigorous approaches to compute higher-order corrections to the Born-Oppenheimer approximation have been developed, most notably in Ref. Hagedorn (1988). The techniques we use in this article are somewhat reminiscent of van Vleck’s perturbation theory and contact transformation often used in chemistry and physics. We use here a compact and powerful notation which will allow us to obtain not only second- but also third-order correction formulae for a single or multi-dimensional (non-adiabatically coupled) electronic subspace.

We believe that the explicit formulation of the effective non-adiabatic Hamiltonians for coupled electronic states, including the earlier missing kinetic (or mass) correction terms, will be useful for the chemical physics community. Their numerical application assume the computation of non-adiabatic coupling vectors. Earlier work in which the mass-correction terms were computed for a single electronic state, e.g., Refs. Mátyus (2018a, b), can be generalized for a multi-state band, so numerical applications will probably follow this theoretical work in the near future.

Summary of the main result.

At the end of this introduction, we summarize the main result of the paper to help orientation in the rather technical sections to come. In the following sections, we derive the general form of the effective non-adiabatic Hamiltonian $\widehat{\boldsymbol{H}}^{(2)}$ for a group of $d$ electronic levels $E_{1}(R),\ldots,E_{d}(R)$ that are separated by a gap from the rest of the spectrum. Our analysis implies for example that its eigenvalues approximate the eigenvalues of the full molecular Hamiltonian up to order $\varepsilon^{3}$, where $\varepsilon=\sqrt{\frac{m}{M}}$ is the square root of the mass ratio of electron and nuclear mass. It thus captures all second-order contributions. It is important to note that the perturbative expansion is carried out without assuming small nuclear momenta, so the nuclear kinetic energy $\|\varepsilon\nabla_{R}\psi\|^{2}$ is of order one and not of order $\varepsilon^{2}$.

After choosing $d$ electronic states $\psi_{1}(R),\ldots,\psi_{d}(R)$ that are smooth functions of $R$ and pointwise form an orthonormal basis of the selected electronic subspace, the projection onto which we denote by $P(R)$, (i.e., an adiabatic or diabatic basis set for the selected electronic subspace), the effective non-adiabatic Hamiltonian $\widehat{\boldsymbol{H}}^{(2)}$ takes the form of an operator acting on wave functions on the nuclear configuration space $\mathbb{R}^{3N}$ that take values in $\mathbb{C}^{d}$, and thus can be written as a $d\times d$-matrix of operators $(\widehat{\boldsymbol{H}}^{(2)})_{\alpha\beta}$ acting on functions on $\mathbb{R}^{3N}$:

[TABLE]

Here the boldface objects are $(d\times d)$ matrix-valued functions on the nuclear configuration space, with $(\mathbf{1})_{\alpha\beta}:=\delta_{\alpha\beta}$ denoting the identity matrix, and the others given as follows in terms of the electronic states $\psi_{1}(R),\ldots,\psi_{d}(R)$.

The coefficients of the non-abelian Berry connection, are as expected, $\mathbf{A}_{\alpha\beta,i}(R)=-\mathrm{i}\langle\psi_{\alpha}(R)|\partial_{i}\psi_{\beta}(R)\rangle$. The ‘diabatic’ electronic level matrix becomes $\mathbf{E}_{\alpha\beta}(R)=\langle\psi_{\alpha}(R)|{H_{\rm e}}(R)|\psi_{\beta}(R)\rangle$, where ${H_{\rm e}}(R)$ is the electron Hamiltonian for fixed nuclear configuration $R$. The second-order diagonal correction is $\mathbf{\Phi}_{\alpha\beta}(R)=\tfrac{1}{2}\sum_{i=1}^{3N}\langle\partial_{i}\psi_{\alpha}(R)|P^{\perp}(R)|\partial_{i}\psi_{\beta}(R)\rangle$, where $P^{\perp}(R)=1-P(R)$ projects on the orthogonal complement of the selected electronic subspace, i.e., on the orthogonal complement of the span of $\psi_{1}(R),\ldots,\psi_{d}(R)$.

While the matrix versions of terms discussed up to now could have been easily guessed from the single band ($d=1$) case, the determination of the second-order mass correction matrix requires the systematic perturbation approach developed in the following sections. The resulting expression is

[TABLE]

where for better readability we dropped the argument $R$ in all the functions. Here ${\mathcal{R}}_{a}(R):=({H_{\rm e}}(R)-E_{a}(R))^{-1}\,P^{\perp}(R)$ is the reduced resolvent of the level $E_{a}(R)$ acting as a bounded operator on the range of $P^{\perp}(R)$, and $P_{a}(R)$ is the projection onto the eigenspace of ${H_{\rm e}}(R)$ corresponding to the eigenvalue $E_{a}(R)$. In the special case that $\psi_{1}(R),\ldots,\psi_{d}(R)$ form an adiabatic basis set, i.e. ${H_{\rm e}}(R)\psi_{\alpha}(R)=E_{\alpha}(R)\psi_{\alpha}(R)$ for $\alpha=1,\ldots,d$, the expression for the mass correction term simplifies to $\mathbf{M}_{ab,ij}(R)=\langle\partial_{j}\psi_{a}(R)|{\mathcal{R}}_{a}(R)+{\mathcal{R}}_{b}(R)|\partial_{i}\psi_{b}(R)\rangle$.

II Hamiltonian, coupling parameter, and operator orders

Molecular Hamiltonian.

The molecular Hamiltonian is the sum of the nuclear kinetic energy acting on the nuclear coordinates $\tilde{R}$, the electron kinetic energy acting on the electronic coordinates $r$, and the Coulomb interaction terms, (in Hartree atomic units, $m_{\text{e}}=\hbar=1$)

[TABLE]

By absorbing the different masses of the atomic nuclei in the mass-scaled Cartesian coordinates $R_{i}=M_{i}^{-1/2}\tilde{R}_{i}$, where $M_{i}$ are the nuclear masses in atomic mass units Cohen et al. (2007), one can re-write the nuclear kinetic energy operator as

[TABLE]

and we label $\varepsilon^{2}=m_{\text{e}}/m_{\text{u}}\ll 1$ the conversion factor between the electronic and the atomic mass scale, which is characteristic to the three orders of magnitude mass difference of the electrons and the atomic nuclei.

Using the common notation $H_{\text{e}}(R):=-\sum_{k=1}^{n}\frac{1}{2}\Delta_{r_{k}}+V(R,r)$ for the electronic Hamiltonian we can re-write the electron-nucleus Hamiltonian into the compact form

[TABLE]

which highlights the $\varepsilon^{2}$ scale-separation (and coupling) between the electron-nucleus quantum mechanical motion. Also note that capital letters without hat label operators that are functions of $R$, i.e., they act fibrewise (pointwise in $R$) on the electronic Hilbert space, such as $H_{\text{e}}$. All other operators (which include differential operators of $R$) are labelled with a wide hat, and the nuclear kinetic energy ${\widehat{K}}$ is an example for this type of operators.

Counting operator orders.

During the course of this work, we will perform an asymptotic expansion of operators in powers of the small parameter $\varepsilon$. Since many of the operators appearing in the calculations are unbounded, we emphasize that we are interested in the action of operators on typical molecular wave functions with energies of order ${\mathcal{O}}(1)$. In particular, the nuclear kinetic energy $-\frac{\varepsilon^{2}}{2}\Delta=\frac{1}{2}\,\widehat{p}^{\,2}$ and thus also the nuclear momentum operator $\widehat{p}=-\mathrm{i}\varepsilon\nabla$ are of order ${\mathcal{O}}(1)$ (instead of ${\mathcal{O}}(\varepsilon)$, which one could naïvely think). This is because a typical molecular wave function shows oscillations with respect to the nuclear coordinates on a spatial scale of order $\varepsilon^{-1}$ and thus it has derivatives of order $\varepsilon^{-1}$. However, when $\widehat{p}$ acts on a smooth, perhaps operator-valued, function $f(R)$ of the nuclear coordinates, we have a quantity ${\mathcal{O}}(\varepsilon)$, indeed:

[TABLE]

These observations will be important to remember for the following calculations.

III The zeroth-order non-adiabatic Hamiltonian:

truncation error and a strategy for reducing this error

Let $\{E_{a}(R),a=1,\ldots,{d}\}$ be a finite set of eigenvalues of the electronic Hamiltonian ${H_{\rm e}}(R)$ that are isolated by a finite gap 111A precise formulation of the gap assumption is as follows. We say that there is a gap of size $g>0$ in the relevant region $\Omega\subset\mathbb{R}^{3N}$ of nuclear configuration space if there are two continuous functions $f_{\pm}:\Omega\to\mathbb{R}$ such that dist$(f_{\pm}(R),\sigma({H_{\rm e}}(R)))\geq g/2$ such that the interval $I(R):=[f_{-}(R),f_{+}(R)]$ satisfies $I(R)\cap\sigma({H_{\rm e}}(R))=\{E_{1}(R),\ldots,E_{d}(R)\}$. Then, all correction terms (Berry phase, diagonal correction, effective mass) are asymptotically small as $\varepsilon\to 0$ when choosing any smooth diabatic basis $\{\psi_{1}(R),\ldots,\psi_{d}(R)\}$. In a real problem, of course, $\varepsilon$ is fixed by the physical parameters and smallness of the correction terms is an indicator for a sufficiently large spectral gap. from the rest of the spectrum over the relevant range of nuclear configurations and denote by $P_{a}(R)$ the corresponding spectral projections. Then $P(R)=\sum_{a=1}^{d}P_{a}(R)$ projects onto the selected electronic subspace and $P^{\perp}:=1-P$ is the projection onto the orthogonal complement.

The full electron-nucleus Hamiltonian can be written in a block form as

[TABLE]

where ${{\widehat{H}}_{P}^{(0)}}:=P{\widehat{H}}P$ corresponds to the ‘usual’ non-adiabatic Hamiltonian 222Note that we follow here the physical chemistry terminology in which a description is called non-adiabatic if it concerns more than one electronic states which are coupled among each other. Compare this terminology with the mathematics’ naming conventions, in which BO or ‘adiabatic’ is commonly used for a description in which a (single or multi-dimensional) electronic subspace is not coupled with the rest of the electronic spectrum.. In numerical computations ${{\widehat{H}}_{P}^{(0)}}$ is usually represented over some adiabatic or diabatic basis, whereas the the off-diagonal $P{\widehat{H}}P^{\perp}$ and $P^{\perp}{\widehat{H}}P$ blocks have been neglected in earlier work.

The following calculation shows that the off-diagonal terms are indeed small (remember the gap condition at the beginning of this section), more precisely, of order $\varepsilon$:

[TABLE]

where for any fibred operator $A$ we use the abbreviations $A^{\prime}:=\nabla A$, and the symmetrized directional derivative of $A$ is

[TABLE]

i.e., $A^{\prime}$ is again a fibred operator, while $\widehat{A}^{\prime}_{\nabla}$ is a first-order differential operator in the nuclear coordinates of ${\mathcal{O}}(1)$. For the upper-right off-diagonal block we find similarly

[TABLE]

Hence, the off-diagonal part of the nuclear kinetic energy operator $\widehat{K}$ is of order $\varepsilon$ and we introduce the abbreviation

[TABLE]

which, as explained above, is an operator of ${\mathcal{O}}(1)$. Using ${\widehat{Q}}$ and the properties of $P$ and $P^{\perp}$, the off-diagonal blocks of the Hamiltonian are written in the compact form

[TABLE]

which are both of ${\mathcal{O}}(\varepsilon)$, and thus

[TABLE]

This expression confirms the well-known fact, the off-diagonal non-adiabatic couplings are small (for a group of bands separated by a gap from the rest of the electronic spectrum). As a consequence, the spectrum of ${{\widehat{H}}_{P}^{(0)}}$ provides an ${\mathcal{O}}(\varepsilon)$ approximation, at least locally in energy, to the spectrum of ${\widehat{H}}$. More precisely, within a neighbourhood of order $\varepsilon$ around any spectral value of ${{\widehat{H}}_{P}^{(0)}}$, there is also a spectral value of the full Hamiltonian ${\widehat{H}}$ 333One simple way to see this is as follows: Let $\psi$ be a normalized eigenfunction of ${{\widehat{H}}_{P}^{(0)}}$, i.e. $({{\widehat{H}}_{P}^{(0)}}-E)\psi=0$ and $\|\psi\|=\|P\psi\|=1$. Then $\chi:=({\widehat{H}}-E)\psi$ has norm $\|\chi\|=\varepsilon\|P^{\perp}{\widehat{Q}}\psi\|$ of order $\varepsilon$ and thus $\|({\widehat{H}}-E)^{-1}\frac{\chi}{\|\chi\|}\|=\frac{\|\psi\|}{\|\chi\|}\gtrsim\frac{1}{\varepsilon}$. Hence dist$(E,{\rm spec}({\widehat{H}}))=\frac{1}{\|({\widehat{H}}-E)^{-1}\|}\lesssim\varepsilon$. For spectral values in the continuous spectrum of ${{\widehat{H}}_{P}^{(0)}}$ one can use the same argument with a Weyl sequence instead of an eigenfunction..

To obtain a better approximation, we replace the projection $P$ by a slightly ‘tilted’ projection ${\widehat{\Pi}}$, such that the off-diagonal terms in the block-decomposition of ${\widehat{H}}$ with respect to ${\widehat{\Pi}}$ are lower order, namely

[TABLE]

for some $n\geq 1$. The projection ${\widehat{\Pi}}$ is obtained from $P$ through a near-identity unitary transformation,

[TABLE]

where the generator ${\widehat{S}}\approx{\widehat{A}}_{1}+\varepsilon{\widehat{A}}_{2}+\varepsilon^{2}{\widehat{A}}_{3}+\ldots$ will be determined exactly by the condition that the off-diagonal elements in Eq. (12) are of order $\varepsilon^{n+1}$. The $n$th-order effective Hamiltonian ${{\widehat{H}}_{P}^{(n)}}$ obtained in this way will give an ${\mathcal{O}}(\varepsilon^{n+1})$ approximation to the spectrum of ${\widehat{H}}$. The physical picture behind the projection ${\widehat{\Pi}}$ is the following. The range of the adiabatic projection $P$ is spanned by states of the form $\Psi(R,r)=\varphi(R)\psi_{a}(R,r)$, where $\psi_{a}(R,r)$ are eigenstates of the electronic Hamiltonian $H_{\text{e}}(R)$ for the clamped nuclear configuration $R$, $H_{\text{e}}(R)\psi_{a}(R,\cdot)=E_{a}(R)\psi_{a}(R,\cdot)$. However, since the nuclei are also moving, the molecular eigenstates are only approximately but not exactly of this local product form. Loosely speaking, the state of the electrons depends also on the momenta of the nuclei. This effect is taken care of by slightly tilting the projection $P$ into the projection ${\widehat{\Pi}}$.

For computing eigenvalues of the $n$th-order effective Hamiltonian ${{\widehat{H}}_{P}^{(n)}}$, it will be more appropriate to consider the unitarily equivalent (and thus isospectral) operator:

[TABLE]

which, by choosing a basis representation for $P$, will provide an $n$th-order effective Hamiltonian for the quantum nuclear motion corresponding to the selected electronic subspace. The eigenvectors of ${{\widehat{H}}_{P}^{(n)}}$ and ${{\widehat{\mathcal{H}}}_{P}^{(n)}}$ are related by the unitary transformation ${\mathrm{e}}^{\mathrm{i}\varepsilon{\widehat{S}}}$.

IV Calculation of the $n$th-order transformation matrix

We would like to achieve an ${\mathcal{O}}(\varepsilon^{n+1})$ block diagonalization of the electron-nucleus Hamiltonian by conjugating it with an appropriate unitary operator $\text{e}^{\text{i}\varepsilon{\widehat{S}}}$, a strategy that is somewhat reminiscent of what is known in chemical physics as van Vleck’s perturbation theory and contact transformation. The actual procedure we use is known as adiabatic perturbation theory in the mathematical literature Martinez and Sordoni (2002); Teufel (2003); Panati et al. (2007); Martinez and Sordoni (2009), but was applied before in physics also in the context of the Born–Oppenheimer approximation in Ref. Weigert and Littlejohn (1993). While the quoted references are, with the exception of Ref. Panati et al. (2007), use pseudodifferential calculus and are therefore quite technical and demanding, our approach is rather elementary and uses a general and compressed notation. As a consequence, we are able to derive with relatively little effort explicit expressions not only for ${{\widehat{\mathcal{H}}}_{P}^{(2)}}$ but also for ${{\widehat{\mathcal{H}}}_{P}^{(3)}}$ and not only for a single but for a finite number of coupled electronic states. Note that ${{\widehat{\mathcal{H}}}_{P}^{(2)}}$ for a single electronic state was explicitly calculated earlier, e.g., in Ref. Weigert and Littlejohn (1993); Teufel (2003); Panati et al. (2007); Bunker and Moss (1977); Pachucki and Komasa (2009); Scherrer et al. (2017), although in some cases without carefully paying attention to the subtleties of counting operator orders in relation with the nuclear momenta (see Sec. II.b).

The simplicity and generality of our derivation is based on the concept and simple algebraic properties of diagonal and off-diagonal operators and an explicit expression for the inverse of the quantum Liouvillian (vide infra) acting on operators which do not contain nuclear differential operators, as to be explained in the upcoming section.

IV.1 Technical preliminaries

Before constructing ${\widehat{S}}$ and ${{\widehat{\mathcal{H}}}_{P}^{(n)}}$, we introduce the concept of (off)diagonal operators and the (inverse) Liouvillian as well as some of their properties, which will become useful during the course of the calculations. Further useful relationships are collected in the Appendix.

IV.1.1 Diagonal and off-diagonal operators

We define the diagonal (${\rm D}$) and the off-diagonal (${\rm OD}$) parts of a linear operator $\widehat{A}$ with respect to the orthogonal projections $P$ and $P^{\perp}=1-P$ as

[TABLE]

respectively. From this definition, a number of simple relationships follow immediately. For example, it holds for all operators ${\widehat{A}}$ and $\widehat{B}$ that

[TABLE]

and

[TABLE]

Finally, note that ${\widehat{Q}}$ defined in Eq. (9) is related to the off-diagonal part of the nuclear kinetic energy as

[TABLE]

IV.1.2 Commutators, the quantum Liouvillian and its inverse

Given an operator $\widehat{B}$, let us define the linear mapping

[TABLE]

whenever the commutator is well defined (for the operators we consider in the following sections, no problems with operator domains occur). For $\widehat{B}={H_{\rm e}}$ we call this mapping the quantum Liouvillian of $\widehat{A}$:

[TABLE]

Since ${H_{\rm e}}$ commutes with $P$ and $P^{\perp}$, the Liouvillian $\mathscr{L}_{H_{\rm e}}$ preserves (off)diagonality, i.e.,

[TABLE]

Furthermore, $\mathscr{L}_{H_{\rm e}}$ is invertible on the space of fibred off-diagonal operators (labelled without any hat). Explicitly, for $B:=B^{\rm OD}$ its inverse is

[TABLE]

where

[TABLE]

is the reduced resolvent (with the nuclear coordinate dependence shown explicitly). Note that, due to the gap condition, $({H_{\rm e}}(R)-E_{a}(R))^{-1}$ is indeed a bounded operator when restricted to states in the orthogonal complement of the range of $P(R)$.

In the following lines, we check that Eq. (22) is indeed the inverse Liouvillian on off-diagonal fibred operators, but we carry out this calculation for a general operator, $\widehat{B}$, which may contain also the nuclear momentum operator (highlighted with a hat in the notation):

[TABLE]

where $E_{a}(R)$ is the electronic energy. For off-diagonal, fibred operators $\mathscr{I}_{H_{\rm e}}$ is indeed the (exact) inverse of $\mathscr{L}_{H_{\rm e}}$, since $[E_{a},B]=0$. Otherwise, it is an approximation to the inverse with an error depending on the value of the commutator $[E_{a},\widehat{B}]$ . For example, if $\widehat{B}={\widehat{K}}^{\rm OD}$, the commutator in Eq. (24) is $[E_{a}(R),\widehat{K}^{\rm OD}]=\varepsilon[E_{a}(R),{\widehat{Q}}]$, and thus $\mathscr{L}_{H_{\text{e}}}(\mathscr{I}_{H_{\text{e}}}({\widehat{K}}^{\rm OD}))={\widehat{K}}^{\rm OD}+{\mathcal{O}}(\varepsilon)$, so the inverse is obtained with an ${\mathcal{O}}(\varepsilon)$ error.

IV.2 Conditions for making the off-diagonal block of the Hamiltonian lower order

As explained in the previous section, in order to reduce the off-diagonal coupling between the selected electronic subspace $P$ and its orthogonal complement $P^{\perp}$, we will choose the self-adjoint operator ${\widehat{S}}={\widehat{A}}_{1}+\varepsilon{\widehat{A}}_{2}+\ldots$ such that the off-diagonal part (coupling) within the ${\widehat{\Pi}}$-block decomposition of the Hamiltonian ${\widehat{H}}$, Eq. (12), is small:

[TABLE]

Due to hermiticity, this condition also implies that ${\widehat{\Pi}}{\widehat{H}}{\widehat{\Pi}}^{\perp}$ is of ${\mathcal{O}}(\varepsilon^{n+1})$. To construct explicitly the operators ${\widehat{A}}_{1},{\widehat{A}}_{2},\ldots$, it is more practical to use the unitary transform

[TABLE]

of condition (25), where the transformed Hamiltonian is defined as

[TABLE]

In the language of diagonal and off-diagonal operators, the condition of Eq. (26) is fulfilled if the off-diagonal part of the transformed Hamiltonian, $\widehat{\mathcal{H}}^{\rm OD}=P^{\perp}\widehat{\mathcal{H}}P+P\widehat{\mathcal{H}}P^{\perp}$, is small:

[TABLE]

which will be our working equation to determine the operators ${\widehat{A}}_{1},\ldots,{\widehat{A}}_{n}$ up to the $(n+1)$st order in $\varepsilon$. Then, using these operators, an explicit expression will be derived for the relevant block of the transformed Hamiltonian, namely of ${{\widehat{\mathcal{H}}}_{P}^{(n)}}=P\,\widehat{\mathcal{H}}P$ as defined in Eq. (III) and explained in Section III.

IV.3 Reduction of the off-diagonal coupling: determination of ${\widehat{A}}_{1},\ldots,{\widehat{A}}_{n}$

The transformed electron-nucleus Hamiltonian is expanded in terms of increasing powers of $\varepsilon$ (see Appendix) as

[TABLE]

For the electronic Hamiltonian the expansion up to ${\mathcal{O}}(\varepsilon^{4})$ is

[TABLE]

Anticipating that commutators of the form $[{\widehat{A}}_{j},{\widehat{K}}]$ are of order $\varepsilon$, the expansion for the nuclear kinetic energy term up to the same order is

[TABLE]

Note that ${\widehat{K}}={\widehat{K}}^{\rm D}+\varepsilon\hat{Q}$, so its diagonal part is of leading order, while the off-diagonal part is ${\mathcal{O}}(\varepsilon)$.

Thus, the transformed Hamiltonian has an asymptotic expansion in powers of $\varepsilon$,

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

To reduce the off-diagonal coupling, we will now proceed by induction. Assuming that ${\widehat{A}}_{1},{\widehat{A}}_{2},\ldots,{\widehat{A}}_{n-1}$ have been chosen such that

[TABLE]

is ${\mathcal{O}}(\varepsilon^{n})$, we will fix ${\widehat{A}}_{n}$ such that $(\widehat{\mathcal{H}}^{(n)})^{\rm OD}=\sum_{i=0}^{n}\varepsilon^{i}\widehat{h}_{i}^{\rm OD}=:\varepsilon^{n+1}\widehat{B}_{n}$ is ${\mathcal{O}}(\varepsilon^{n+1})$.

We mention already at this point that fulfillment of this sequence of requirements will fix only the off-diagonal part ${\widehat{A}}_{i}^{\rm OD}$ of each ${\widehat{A}}_{i}$. The diagonal parts, ${\widehat{A}}_{i}^{\rm D}$, generate merely rotations within the subspaces $P$ and $P^{\perp}$, but do not affect the (de)coupling. Hence, we set ${\widehat{A}}_{i}^{\rm D}=0$ ($i=1,2,\ldots$). With this choice, ${\widehat{A}}_{1},{\widehat{A}}_{2},\ldots,{\widehat{A}}_{n}$ are completely determined by the requirement, Eq. (28).

Zeroth-order off-diagonal (OD) terms: The off-diagonal part of ${\widehat{h}}_{0}$, Eq. (33),

[TABLE]

is of order $\varepsilon$, hence decoupling is automatically fulfilled at this order.

First-order OD terms: In the next step we require

[TABLE]

where we used Eq. (21) and the fact that ${\widehat{A}}_{1}={\widehat{A}}_{1}^{\rm OD}$. This condition can be fulfilled by choosing

[TABLE]

Inserting Eq. (40) back into Eq. (39) and using Eq. (24), we find that

[TABLE]

This explicit expression for $\widehat{B}_{1}$ will be required for calculating the third-order terms in the effective Hamiltonian ${{\widehat{\mathcal{H}}}_{P}^{(3)}}$ (see Eqs. (57)–(67)).

Second-order OD terms: We require

[TABLE]

where in the third equality we used the algebraic relations of Eq. (17). Again, we solve this equation for ${\widehat{A}}_{2}$ using the approximate inverse Liouvillian, Eq. (22):

[TABLE]

The explicit expression for the remainder $\widehat{B}_{2}$ could be determined, if needed, through the calculation of the ${\mathcal{O}}(\varepsilon)$ error term, Eq. (24), from the approximate inversion:

[TABLE]

Third-order OD terms: By the same reasoning we require

[TABLE]

where we used that the last two terms in the first line are zero (off-diagonal part of diagonal operators). Again, we make the left hand side of Eq. (45) small by solving the equation for ${\widehat{A}}_{3}$ using the approximate inverse Liouvillian, Eq. (22):

[TABLE]

and the remainder term $\widehat{B}_{3}$ can be determined by the direct calculation of the ${\mathcal{O}}(\varepsilon)$ error of the inversion using Eq. (24).

It is obvious how to continue this induction to arbitrary orders. However, as we are only interested in explicit expressions for the effective Hamiltonians up to third order, we refrain from stating the general induction explicitly.

IV.4 Second- and third-order Hamiltonians

In this section we calculate the leading terms in the expansion of the $n$th-order effective Hamiltonian

[TABLE]

up to and including $P\,\widehat{h}_{3}P$, and thus obtain explicit expressions for the second and the third-order effective Hamiltonians, ${{\widehat{\mathcal{H}}}_{P}^{(2)}}$ and ${{\widehat{\mathcal{H}}}_{P}^{(3)}}$. To this end, we first calculate the diagonal parts $\widehat{h}_{j}^{\rm D}$ and then, in a second step, project onto the range of $P$.

Zeroth-order diagonal (D) terms:

[TABLE]

First-order D terms: Recalling that ${\widehat{A}}_{j}={\widehat{A}}_{j}^{\rm OD}$ for all $j\geq 1$ and the algebraic relations in Eq. (17), we find, in particular, that

[TABLE]

Second-order D terms: Similarly, since $\widehat{B}_{j}=\widehat{B}_{j}^{\rm OD}$ for all $j\geq 1$, we find

[TABLE]

Third-order D terms:

[TABLE]

By combining these expressions, we obtain the second- and the third-order effective Hamiltonians as

[TABLE]

and

[TABLE]

respectively. We note that in the third-order correction, the $-\tfrac{1}{2}\mathscr{L}_{{\widehat{A}}_{1}}(\widehat{B}_{1})$ remainder from second order, cancels the $\tfrac{1}{2}\mathscr{L}_{{\widehat{A}}_{1}}(\widehat{B}_{1})$ third-order term. In the following subsection, we continue with inserting the explicit formulae for ${\widehat{A}}_{1}$, ${\widehat{A}}_{2}$, and ${\widehat{Q}}$ into the compact expressions of ${{\widehat{\mathcal{H}}}_{P}^{(2)}}$ and ${{\widehat{\mathcal{H}}}_{P}^{(3)}}$ just obtained.

IV.5 More explicit expressions for the second- and

the third-order non-adiabatic Hamiltonian corrections

Second-order correction.

Using the explicit expression for the inverse Liouvillian, Eq. (22), Eq. (40) yields

[TABLE]

Thus, the ${\mathcal{O}}(\varepsilon^{2})$ correction term of the effective Hamiltonian in Eq. (52) is

[TABLE]

where we inserted Eq. (9) for ${\widehat{Q}}$ and used the fact that the reduced resolvent ${\mathcal{R}}_{a}$ acts only in the $P^{\perp}$ subspace. (Recall that $P^{\prime}_{\nabla}$ is defined by Eq. (7).) Since the commutator of $\varepsilon\nabla$ in $P^{\prime}_{\nabla}$ with smooth, fibred operators yields higher-order terms in $\varepsilon$, we can further simplify ${{\widehat{\mathcal{H}}}_{P}^{(2)}}$ to

[TABLE]

Note, however, that the ${\mathcal{O}}(\varepsilon^{3})$ term does contribute to the third-order effective Hamiltonian and can not be neglected when computing ${{\widehat{\mathcal{H}}}_{P}^{(3)}}$.

Third-order correction.

To obtain an explicit expression for the third-order correction, we need to derive an explicit expression for ${\widehat{A}}_{2}$, Eq. (43):

[TABLE]

which assumes the explicit knowledge of the first-order remainder term $\widehat{B}_{1}$, Eq. (41), too. The first term in Eq. (57) includes

[TABLE]

and its inverse Liouvillian, Eq. (22), is

[TABLE]

The inverse Liouvillian of $\widehat{B}_{1}$, Eq. (41), is

[TABLE]

Next, the explicit expression for ${\widehat{A}}_{2}$, obtained as the sum of Eqs. (59) and (60), is used to expand the third-order correction as

[TABLE]

By working out the second sum of Eq. (61), we obtain

[TABLE]

where the last expression exactly cancels the first term in Eq. (61), and thus the correction term at third order is

[TABLE]

The commutators can be evaluated as

[TABLE]

and

[TABLE]

with $P^{\prime\prime}_{\nabla^{2}}:=\varepsilon\sum_{i,j=1}^{3N}\partial_{j}(\partial_{j}\partial_{i}P)\varepsilon\partial_{i}$, and

[TABLE]

Inserting these identities into Eq. (63), we finally obtain for the third-order correction as

[TABLE]

Note that the first three lines are third order in the nuclear momentum, $\widehat{p}$, and the last line is linear in $\widehat{p}$.

When looking at the second- and the third-order corrections in Eq. (56) and (67), one might worry about singular expressions. Indeed, some of the summands become singular near points of the nuclear configuration space at which eigenvalues within the set $\{E_{a}\,|\,a=1,\ldots,d\}$ cross. At these points the single spectral projections $P_{a}$ might not be differentiable. However, as it can be seen from the original expressions, Eqs. (52) and (53), for ${{\widehat{\mathcal{H}}}_{P}^{(2)}}$ and ${{\widehat{\mathcal{H}}}_{P}^{(3)}}$, the complete expression (the full sum) remains bounded because the singularities in the different summands cancel each other. This property might require additional care in numerical computations.

V Basis representation and effective nuclear Hamiltonians

By choosing an electronic (e.g., adiabatic or diabatic) basis set $\{\psi_{\alpha},\alpha=1,\ldots,{d}\}$ for $P{\mathfrak{H}_{\rm e}}$, one can represent a molecular wave function $\Psi$ in the range of $P$, $\Psi\in P{\mathfrak{H}}$, as $\Psi(R,r)=\sum_{\alpha=1}^{{d}}\varphi_{\alpha}(R)\psi_{\alpha}(R,r)$ (where ${\mathfrak{H}_{\rm e}}$ and ${\mathfrak{H}}$ denote the electronic and the molecular Hilbert space, respectively). It is common practice to represent the zeroth-order effective Hamiltonian ${{\widehat{\mathcal{H}}}_{P}^{(0)}}={{\widehat{H}}_{P}^{(0)}}=P{\widehat{H}}P$ as a matrix operator ${\widehat{\boldsymbol{H}}}_{P}^{(0)}$ with respect to such a basis set, which then acts only on the nuclear functions $(\varphi_{1}(R),\ldots,\varphi_{{d}}(R))$. This yields, in particular, also the Berry phase and the diagonal BO correction terms (see below).

In what follows, we construct the matrix representation also for the second-order Hamiltonian ${{\widehat{\mathcal{H}}}_{P}^{(2)}}$. As a special case, the known mass-correction terms for a single, isolated electronic state will be recovered. The basis representation for the third-order correction, Eq. (67), can be worked out along the same lines.

V.1 Basis representation for the second-order, multi-state Hamiltonian

Let us choose an electronic (e.g., adiabatic or diabatic) basis set $\psi_{1}(R),\ldots,\psi_{d}(R)$ such that the $\psi_{\alpha}$ are smooth functions of $R$ and pointwise form an orthonormal basis of the range of $P(R)$, i.e.,

[TABLE]

Because of the gap condition, such a smooth diabatic basis set always exist 444Strictly speaking, such a smooth diabatic basis set exists over each contractible subset $\Omega\subseteq\mathbb{R}^{3N}$ of the nuclear configuration space $\mathbb{R}^{3N}$ on which the gap condition is satisfied. This is because the gap condition implies that the projections $P(R)$ define a smooth rank-$d$ vector bundle over $\Omega$, and a vector bundle over a contractible set always allows for a global trivialisation, that is for a diabatic basis set as above.. However, due to possible crossings within the set of eigenvalues $E_{1},\ldots,E_{d}$, it might not be possible to choose $\psi_{1}(R),\ldots,\psi_{d}(R)$ as smooth functions of $R$, and at the same time, as pointwise eigenfunctions of $H_{\text{e}}(R)$, i.e.,

[TABLE]

Then, the matrix representation of ${{\widehat{\mathcal{H}}}_{P}^{(2)}}$, Eq. (56), over $\psi_{\alpha}$, $\alpha=1,\ldots,{d}$, results in a matrix operator ${\widehat{\boldsymbol{H}}}_{P}^{(2)}$ for the quantum nuclear motion with matrix elements

[TABLE]

For the kinetic-energy part, we find

[TABLE]

By introducing the coefficient of the non-abelian Berry-connection

[TABLE]

we find

[TABLE]

With $\langle\psi_{\alpha}|\partial_{i}\psi_{\gamma}\rangle=-\langle\partial_{i}\psi_{\alpha}|\psi_{\gamma}\rangle$, the last term becomes

[TABLE]

and thus the kinetic-energy term can be written in the form

[TABLE]

with the matrix-valued Berry-connection coefficient $\mathbf{A}_{i}$ and the matrix-valued potential energy correction

[TABLE]

This latter quantity can be understood as an ${\mathcal{O}}(\varepsilon^{2})$ correction to the ‘diabatic’ electronic level matrix

[TABLE]

Note that $\mathbf{E}_{\alpha\beta}$ is a diagonal matrix if and only if all $\psi_{\alpha}$ are eigenvectors of ${H_{\rm e}}$. Defining the second-order mass-correction term as

[TABLE]

the matrix representation of ${{\widehat{\mathcal{H}}}_{P}^{(2)}}$ over $\psi_{\alpha}$, $\alpha=1,\ldots,{d}$, can be compactly written as

[TABLE]

We note that this is the complete second-order non-adiabatic Hamiltonian operator for the nuclear motion. It is important to remember the peculiarities of counting operator orders (Sec. II.c), which follow from not making the assumption that the nuclear momenta (when $-\text{i}\varepsilon\partial_{R}$ acts on the nuclear wave function) are small. Also note that we used Cartesian coordinates scaled with the nuclear mass, Eq. (2), so the derived expressions can be used also for heteronuclear systems by making this scaling factor explicit in the numerical computations.

In the special case where $\psi_{1},\ldots,\psi_{d}$ form an adiabatic basis set, i.e., satisfy Eq. (69), the expression for $\mathbf{M}_{\alpha\beta,ij}$ simplifies to

[TABLE]

To see this, one uses that the reduced resolvent, Eq. (23), contains a projection $P^{\perp}$, and

[TABLE]

and similarly for its adjoint, $P^{\perp}(\partial_{i}P)=\sum_{\gamma=1}^{d}P^{\perp}|\partial_{i}\psi_{\gamma}\rangle\langle\psi_{\gamma}|$.

In all our expressions the nuclear differential operators are written in terms of Cartesian coordinates. The operators can be transformed to curvilinear coordinates, necessary for efficient rovibrational computations, similarly to the transformation of the single-state non-adiabatic Hamiltonian as it was carried out in Ref. Mátyus (2018a) using the Jacobi and the metric tensors of the new coordinates.

For the special case of a single electronic state ($\psi_{1},E_{1}$), we are free to choose a real-valued, normalized electronic wave function $\psi_{1}$. Then, the effective operator of the atomic nuclei, Eq. (79), simplifies to

[TABLE]

So, we assume that $\psi_{1}$ is chosen such that $A_{1}=-\mathrm{i}\langle\psi_{1}|\nabla\psi_{1}\rangle=0$ and find for the mass correction term that

[TABLE]

This mass-correction function is identical with that used in Ref. Mátyus (2018a), and thus, for a single, isolated electronic state, the known expression of the second-order non-adiabatic Hamiltonian is recovered.

VI Summary and conclusions

Molecular wave functions are often approximated on the subspace $P{\mathfrak{H}}$ of the full electron-nucleus ${\mathfrak{H}}$ Hilbert space, where $P$ is the electronic subspace which governs the motion of the atomic nuclei.

We have shown that a complete neglect of the complementary electronic subspace $(1-P){\mathfrak{H}}$ introduces an ${\mathcal{O}}(\varepsilon)$ error in the Hamiltonian and also in the molecular spectrum ($\varepsilon$ is the square root of the electron-to-nucleus mass ratio). We improved upon this ${\mathcal{O}}(\varepsilon)$ approximation, by using a near-identity unitary transform of $P$, ${\widehat{\Pi}}={\mathrm{e}}^{\mathrm{i}\varepsilon{\widehat{S}}}P{\mathrm{e}}^{-\mathrm{i}\varepsilon{\widehat{S}}}$. Terms of the self-adjoint transformation operator ${\widehat{S}}={\widehat{A}}_{1}+\varepsilon{\widehat{A}}_{2}+\varepsilon^{2}{\widehat{A}}_{3}+\ldots$ were determined for increasing orders of $\varepsilon$ by making the coupling, and hence the error of the molecular energy, $\varepsilon$ times smaller at every order. The resulting transformation operators include the momentum operator $\widehat{p}$ of the atomic nuclei, thereby the transformed electronic space ${\widehat{\Pi}}$, which makes the coupling lower order, depends not only on the nuclear positions $R$ but also on the nuclear momenta $\widehat{p}$. The transformed electronic states adjusted by $\widehat{p}$ up to order $\widehat{p}^{n}$, achieve a block-diagonalization of ${\widehat{H}}$ up to terms of order $\varepsilon^{n+1}$. From the transformed, ${\mathcal{O}}(\varepsilon^{n+1})$ block-diagonal Hamiltonian, we obtained effective $n$th-order Hamiltonians for the quantum nuclear motion. Explicit expressions were derived up to the third-order corrections for a multi-dimensional electronic subspace.

In particular, the second-order non-adiabatic Hamiltonian contains correction terms quadratic in the nuclear momenta, which may be small near the bottom of the electronic band, but for highly excited states they can easily dominate the diagonal correction. These kinetic energy correction terms can be identified as a coordinate-dependent correction to the nuclear mass in the nuclear kinetic energy operator. These earlier neglected ‘mass-correction terms’ perturbatively account for the effect of the electronic states not included in the selected, explicitly coupled electronic band. For a single electronic state the multi-state expressions simplifies to the known, second-order Hamiltonian including the mass-correction function.

This perturbative decoupling can be used for isolated (groups of) electronic states and we believe that at least the second-order, multi-state expression will soon gain practical applications in rovibronic and quantum scattering computations. Examples for potential applications include the electronically excited manifold of molecular hydrogen—the first steps towards these applications are reported in Ref. FeM —, the predissociation dynamics of H${}_{3}^{+}$, in which the interaction of the electronic ground and excited states is thought to play a role, and also the H+H2 reactive scattering system.

VII Appendix

Contact transform of an operator.

Transformation of an operator ${\widehat{Y}}$ with ${\mathrm{e}}^{\mathrm{i}\varepsilon{\widehat{S}}}$ can be expanded in terms of increasing powers of $\varepsilon$ as

[TABLE]

Commutator operations with diagonal and off-diagonal operators.

[TABLE]

For the example of the diagonal $H_{\text{e}}$ and the off-diagonal ${\widehat{A}}_{i}$, we have collected the following identities (relevant for the calculations in the manuscript):

[TABLE]

Commutator expressions.

[TABLE]

**Acknowledgment

**EM thanks financial support of the Swiss National Science Foundation through a PROMYS Grant (no. IZ11Z0_166525).