Quantum Chaos in the Dynamics of Molecules

Abstract

Quantum chaos is reviewed from the viewpoint of “what is molecule?”, particularly placing emphasis on their dynamics. Molecules are composed of heavy nuclei and light electrons, and thereby the very basic molecular theory due to Born and Oppenheimer gives a view that quantum electronic states provide potential functions working on nuclei, which in turn are often treated classically or semiclassically. Therefore, the classic study of chaos in molecular science began with those nuclear dynamics particularly about the vibrational energy randomization within a molecule. Statistical laws in probabilities and rates of chemical reactions even for small molecules of several atoms are among the chemical phenomena requiring the notion of chaos. Particularly the dynamics behind unimolecular decomposition are referred to as Intra-molecular Vibrational energy Redistribution (IVR). Semiclassical mechanics is also one of the main research fields of quantum chaos. We herein demonstrate chaos that appears only in semiclassical and full quantum dynamics. A fundamental phenomenon possibly giving birth to quantum chaos is “bifurcation and merging” of quantum wavepackets, rather than “stretching and folding” of the baker’s transformation and the horseshoe map as a geometrical foundation of classical chaos. Such wavepacket bifurcation and merging are indeed experimentally measurable as we showed before in the series of studies on real-time probing of nonadiabatic chemical reactions. After tracking these aspects of molecular chaos, we will explore quantum chaos found in nonadiabatic electron wavepacket dynamics, which emerges in the realm far beyond the Born-Oppenheimer paradigm. In this class of chaos, we propose a notion of Intra-molecular Nonadiabatic Electronic Energy Redistribution (INEER), which is a consequence of the chaotic fluxes of electrons and energy within a molecule.

1. Introductory Remarks: What Makes Molecules Special in Quantum Science and Chaos

Molecules consist of heavy nuclei and light electrons, which give rise to quite interesting and characteristic properties. All of the living bodies and biological systems on the globe are made up with molecules and are subject to the beautiful chemical laws. In the biological activities most of those molecular dynamics and properties should appear to be deterministic in the molecular levels, whereas the general dynamics of molecules often work in stochastic and statistical manners. The origin of such statistical properties arising even from small molecular systems are widely ascribed to classical [1] and quantum [2,3,4,5,6,7,8,9,10,11,12] chaos. This in turn suggests that one may find novel chaos in the dynamics of molecules.

Light electrons in a molecule move very fast and yet are accompanied with long wavelength in their matter waves, while the nuclear motion is much slower, but the relevant wavelengths are relatively far shorter due to the heavy masses. Thus, the two intrinsic hierarchical structures emerge in time-scale and length-scale within the dynamics of a molecule. Consequently the general picture of a molecule are: (i) Nuclei move on mean-field potential functions that are generate by the electrons. (ii) The nuclear motion can bear classical and semiclassical nature due to the short wavelength, while electrons are utterly quantum. Therefore, in contrast to classical chaos of the nuclear motion, it is often claimed that there is no chaos in electronic states because of discreteness in the energy levels with large gaps among them. Physics of the mapping of classical chaos onto quantum systems is often referred to as quantum chaology [13]. As a matter of fact most of the chaos studies in molecular science have been focused on the vibrational spectroscopy and nonlinear mixing of the vibrational modes and the theoretical foundation of statistical laws in chemical reaction dynamics. All these are about the chaoticity arising from the nuclear motion alone.

However, the theoretical separation of the electronic and nuclear motions due to the fast and slow modes, which is known as the Born-Oppenheimer separation, can sometimes break down through large kinematic interaction between them, which cause a mixing of the relevant electronic states. This interaction is generally and vaguely referred to as nonadiabatic interaction [14,15,16,17,18,19,20]. This interaction is extremely important in molecular science, since the electronic states and the relevant molecular properties can change instantaneously and drastically when a molecule passes through those nonadiabatic interaction regions. Moreover, we will show that nonadiabatic mixing among many electronic states involved in the manifold of densely quasi-degenerate electronic states can result in strong chaos, which can manifest as a diffusion motion of electronic states in the Hilbert space. This article will hence discuss chaos in what we call nonadiabatic electron dynamics [20,21,22,23,24,25,26,27].

Chaos theory for Hamilton mechanics was initiated (resumed from the view point of Poincaré) around 1950s [28]. Many beautiful theoretical models, elegant mathematics of geometry along with sophisticated numerical studies on ordinary differential equations have been studied since [1]. As for molecular chaos, there are very many pieces of theoretical studies in the literature, yet, due to the complexity of molecular nature most of the theoretical developments and applications are limited to small dimensional systems, typically two-dimensional in the early days. On the other side, the present author and his group have long been studying the chemical dynamics ranging from nuclear classical and semiclassical dynamics (see ref. [29,30] and more cited there) to full quantum electron dynamics [25,26,27]. Such molecular studies on chemical dynamics often guide us to chaotic phenomena that appear only in molecules. Conversely, tracking the path of our chaos studies, along with some data added herewith, would give partly an idea about what the dynamics of molecule are. We therefore apologize in advance that this review article covers only a very biased aspect of chaos studies. We hence do not cover the topics of the main-stream theories and well-established materials for quantum chaos such as level statistics, the random matrix theory, chaos in external fields and dissipative systems, and so on. Extensive reviews about general quantum chaos theory [11,12] are now widely available in numerous monographs [2,3,4,5,6,7,8,9,10] and review articles published elsewhere. The special issue of this journal Entropy, “Quantum Chaos—Dedicated to Professor Giulio Casati on the Occasion of His 80th Birthday”, is certainly among the newest ones.

We begin this review with the very basic structure of molecular science in Section 2, along with the so-called Born-Oppenheimer approximation, which separates nuclear and electronic dynamics in molecules [31].

In Section 3, we briefly touch upon the aims and implications of chaos studies in molecular science, particularly for chemical dynamics. For this purpose only, we take a couple of examples from previous studies on classical chaos.

Section 4 resumes with the full quantum chaos using the modified Hénon-Heiles potential as a model system of two-dimensional vibrational motions. We track the quantum wavepackets running within the quasi-separatrix and construct some of the associated eigenfunctions. A class of dynamical tunneling that exists only in the quasi-separatrix will be presented. The long-time energy spectra and the dependence of the chaotic spectra on the magnitude of the Planck constant are briefly examined.

Multidimensional semiclassical mechanics is indispensable in the treatment of the dynamics of heavy nuclei in molecules. We review our developed semiclassics [29,30] in Section 5, which we call the method of Action Decomposed Function (ADF) [32,33]. With ADF we analyze the mechanism of quantization not only of forming the energy peaks but of nullification of the off-resonance components. As a consequence of our theoretical and numerical studies, we proceed to the idea of amplitude free quantization without use of the annoying amplitude factors such as the van Vleck determinant [34] and those inherent to the semiclassical kernels (Greens functions) [35,36,37]. Semiclassical tunneling is also touched upon.

Then we study the characteristic molecular science, in which the nuclear dynamics kinematically couples with the electronic states, where the electronic-nuclear separation is broken down. In Section 6, we first concentrate on the nonadiabatic nuclear dynamics, which are subject to coupled equations of motion on multiple potential energy surfaces. The electronic states are still represented in time-independent wavefunctions. We show that the phenomenon of “bifurcation and merging (remixing) of quantum wavepackets”, a fundamental mechanism of quantum chaos, is indeed observable experimentally.

The final part of this review is devoted to the nonadiabatic electron dynamics in Section 7. Electronic wavepacket description is particularly useful in the stage where the adiabatic electronic states behind are embedded in a quasi-degenerate manifold and thereby nonadiabatic mixing among them is very intense and continual. Those strong and enduring electronic state mixing result in quantum chaos for the molecular electronic states.

This review concludes with some remarks in Section 8, stressing that the theory of quantum chaos and philosophy behind is indeed fundamental and useful also in practical science like chemistry.

2. Brief Introduction to the Theoretical Framework of Molecules: Born-Oppenheimer Approximation to Separate Electronic and Nuclear Motions

We first outline briefly the very basic frame of quantum description of molecules suggesting how chaos studies of molecules proceed.

2.1. The Born-Oppenheimer (BO) Approximation

The nonrelativistic molecular Hamiltonian is given as

$$\begin{array}{cc}\hfill H(\mathbf{r},\mathbf{R})& =T_N+H^{el}(\mathbf{r};\mathbf{R})\hfill \\ & =\frac{1}{2}\sum _k\frac{{\widehat{P}}_k^2}{M_k}+H^{\left(el\right)}(\mathbf{r};\mathbf{R})\hfill \end{array}$$

(1)

and

$$i\hslash \frac{\partial }{\partial t}\Psi (\mathbf{r},\mathbf{R},t)=H(\mathbf{r},\mathbf{R})$$

(2)

is the electronic Hamiltonian defined as

$$\begin{array}{cc}\hfill H^{\left(el\right)}(\mathbf{r};\mathbf{R})& =T_e+V_c(\mathbf{r};\mathbf{R})\hfill \\ & =\frac{1}{2m}\sum _j{\widehat{p}}_j^2+V_c(\mathbf{r};\mathbf{R})\hfill \end{array}$$

(3)

with $\mathbf{r}$ and $\mathbf{R}$ representing the electronic and nuclear coordinates, respectively, while ${\widehat{p}}_j$ and ${\widehat{P}}_k$ denote the operators for their conjugate momenta of the jth and kth component of $\mathbf{r}$ (specifically written as $r_j$) and $\mathbf{R}$ ($R_k$), respectively. We here consider only the Coulombic interactions $V_c(\mathbf{r};\mathbf{R})$ among electrons and nuclei, that is

$$V_c(\mathbf{r};\mathbf{R})=\sum _{a<b}\frac{e^2}{\left|{\mathbf{r}}_a-{\mathbf{r}}_b\right|}-\sum _a\sum _A\frac{Z_A{e}^2}{\left|{\mathbf{r}}_a-{\mathbf{R}}_A\right|}+\sum _{A<B}\frac{Z_A{Z}_B{e}^2}{\left|{\mathbf{R}}_A-{\mathbf{R}}_B\right|},$$

(4)

where ${\mathbf{r}}_a$ and ${\mathbf{R}}_A$ are the ath electron and Ath nuclei, respectively, and $Z_A$ indicate the nuclear charge on the Ath nuclei.

Due to the great difference in the time scales between electrons and nuclei, their motion is assumed to be separated at each time, as the so-called the Born-Oppenheimer (BO) approximation represents the total wavefunction $\Psi (\mathbf{r},\mathbf{R},t)$ to be

$$\Psi (\mathbf{r},\mathbf{R},t)\simeq {\chi }_I(\mathbf{R},t){\Phi }_I(\mathbf{r};\mathbf{R}),$$

(5)

where ${\Phi }_I(\mathbf{r};\mathbf{R})$ and ${\chi }_I(\mathbf{R},t)$ are the electronic and nuclear wavefunctions, respectively. Note that electrons are assumed to move so fast that it can adjust themselves as a stationary wave at any nuclear position $\mathbf{R}$, and thereby the electronic wavefunction ${\Phi }_I(\mathbf{r};\mathbf{R})$ is assumed to be free of time coordinate t. This in turn requires that ${\Phi }_I(\mathbf{r};\mathbf{R})$ is to be attained at each $\mathbf{R}$ such that

$$H^{el}(\mathbf{r};\mathbf{R}){\Phi }_I(\mathbf{r};\mathbf{R})=V_I\left(\mathbf{R}\right){\Phi }_I(\mathbf{r};\mathbf{R}).$$

(6)

Notice that $\mathbf{R}$ in ${\Phi }_I(\mathbf{r};\mathbf{R})$ and in Equation (6) is thereby treated as a parameter. The electronic energy $V_I\left(\mathbf{R}\right)$ as a function of nuclear coordinates serves as a potential energy function ($I=1,2,\dots$) for the nuclear wavefunction in such a manner that

$$i\hslash \frac{\partial }{\partial t}{\chi }_I(\mathbf{R},t)=\left[\sum _k\frac{{\widehat{P}}_k^2}{2M_k}\chi +V_I\left(\mathbf{R}\right)\right]{\chi }_I(\mathbf{R},t).$$

(7)

The time independent figure $V_I\left(\mathbf{R}\right)$ is often called the potential energy hyper-surface (PES) in the chemical literature [31].

The lowest-order correction to the Born–Oppenheimer energy is found to be in the order of

$${\left(m/M\right)}^{6/4}.$$

(8)

This analysis [38] clearly demonstrates that the Born–Oppenheimer approximation is far better than would be expected from a simple mass ratio $m/M$. For example, if $m/M\sim {10}^{-4}$, which is typical in molecular systems, the error term is ${\left(m/M\right)}^{1.5}\sim {10}^{-6}$.

2.2. The Born-Huang Expansion

To approach the exact solution of Equation (2) from the BO approximation, we need to expand $\Psi (\mathbf{r},\mathbf{R},t)$ in the stationary electronic basis functions $\left\{{\Phi }_I(\mathbf{r};\mathbf{R})\right\}$ at each nuclear position $\mathbf{R}$ as

$$\Psi (\mathbf{r},\mathbf{R},t)=\sum _I{\chi }_I(\mathbf{R},t){\Phi }_I(\mathbf{r};\mathbf{R}),$$

(9)

where unknown time-dependent coefficients ${\chi }_I(\mathbf{R},t)$ are supposed to describe the nuclear wavepackets. This is referred to as the Born-Huang expansion [39]. The electronic basis functions $\left\{\right|{\Phi }_I\rangle \}$ are assumed to be orthonormal at each nuclear configuration as

$$\langle {\Phi }_I\left(\mathbf{R}\right)|{\Phi }_J\left(\mathbf{R}\right)\rangle ={\delta }_{IJ}.$$

(10)

The nuclear wavefunctions are subject to the following coupled equations

$$i\hslash \frac{\partial }{\partial t}{\chi }_I={\displaystyle \frac{1}{2}}\sum _k\frac{{\widehat{P}}_k^2}{M_k}{\chi }_I+\sum _J{H}_{IJ}^{el}{\chi }_J-i\hslash \sum _k\sum _J\frac{X_{IJ}^k}{M_k}{\widehat{P}}_k{\chi }_J-{\displaystyle \frac{{\hslash }^2}{2}}\sum _k\sum _J\frac{Y_{IJ}^k}{M_k}{\chi }_J,$$

(11)

where

$$X_{IJ}^k=〈\left.{\Phi }_I\right|{\displaystyle \frac{\partial {\Phi }_J}{\partial R_k}}〉,Y_{IJ}^k=〈{\Phi }_I|{\displaystyle \frac{{\partial }^2{\Phi }_J}{\partial {R_k}^2}}〉,$$

(12)

and

$$H_{IJ}^{el}\left(\mathbf{R}\right)=〈{\Phi }_I\left|H^{el}\right|{\Phi }_J〉.$$

(13)

Since $X_{IJ}^k$ and $Y_{IJ}^k$ serve as the mixing elements among the adiabatic electronic states, they are collectively called nonadiabatic interactions. Tough difficulties in an attempt to solve the coupled equations of motion (11) are due to the many dimensionality of $\mathbf{R}$, the oscillatory nature (due to the short wavelength) of $\left\{{\chi }_I\right\}$, and many channels to be taken into account. Besides, some of $\left\{{\chi }_I\right\}$ can extend to asymptotic regions in reaction (scattering) dynamics.

Note that the above Born-Huang expansion can be performed with use of other basis functions that do not satisfy Equation (6). In such a case $H_{IJ}^{el}\left(\mathbf{R}\right)$ in Equation (13), have significant amount of nonzero off-diagonal elements. $\left\{{\Phi }_I(\mathbf{r};\mathbf{R})\right\}$ satisfying Equation (6) are referred to as the adiabatic wavefunctions, while others are to diabatic wavefunctions collectively. There are many ways to determine the diabatic functions [16,17,18,20], a typical set being those that minimize $\left|X_{IJ}^k\right|$. We will use the adiabatic representation in what follow unless otherwise noted.

3. Implication of Chaos in Chemical Dynamics

3.1. Determinicity versus Stochasticity in Molecules

It was Arrhenius who established the well-known thermodynamical reaction rate k from initial molecules to products in the form

$$k\propto exp(-E_a/RT).$$

(14)

This formula is indeed revolutionary in that the rate or speed has been derived based on the Maxwell-Boltzmann distribution for equilibrium states. The physical insight of this rate expression was uncovered by Eyring [40,41,42]. He studied a potential energy surface for simple atomic exchange reaction

$$\mathrm{A}+\mathrm{BC}\to {\left[\mathrm{ABC}\right]}^{‡}\to \mathrm{AB}+\mathrm{C}$$

(15)

where ${\left[\mathrm{ABC}\right]}^{‡}$ is a transient state, called the transition state, at the bottleneck area on the way of reaction. He assumed a statistical distribution for the vibrational modes transversal to the translational one and associated the partition functions with them. He thus clarified the quantum mechanical meaning of the barrier height $E_a$ along with the preexponential factors to k (Polanyi and Wigner also stressed the geometrical roles of the potential energy surface $V_I\left(\mathbf{R}\right)$ of Equation (6)) [40,41,42]. It is now known that the transition state theory is approximately valid even in binary collisions in vacuum. Even a quantum mechanical extension is one of the important subjects [43,44]. On the other hand, rather simple reactions for small systems can now well be treated in full quantum mechanics in a completely deterministic manner. (I do never mean that the full quantum scattering calculation for chemical reactions is an easy task. See refs. [45,46] for examples.) Yet the fundamental question is aimed at the theoretical ground of the statistical assumption made in the transition state.

While being subject to statistical laws on one hand, the dynamics of molecular motions are generally deterministic on the other hand. Let us see such a rather intricate example. See Figure 1, which is an ab initio (first principle) dynamics of 5-methyltropolone (5MTR). Ushiyama found that the proton shift from on O${}_2$ in the panel (a) to the position next to O${}_1$ (panel (c)) induces the rotation of the methyl group -C${}_8$H${}_a$H${}_b$H${}_c$at the bottom. The reciprocal motion of proton turns to rotational motion of the methyl group through such a long distance. Careful inspection at the panel (a) and (c) shows that this determinicity is supported by the shift of double bonds (see ref. [47] for the electronic state mechanism). This example denies the simple minded view that a molecule can be always regarded as a set of points of nuclear point masses connected by “springs” among them. We therefore need a careful thought when we study the chaos of molecular vibrational motions.

Another example of statistical reaction rate theory was given by Marcus [48], which augmented the successors’ formalism due to Rice, Ramsperger, and Kassel, thereby now called the RRKM theory for unimolecular dissociation, which predicts the rate of slow unimolecular dissociation of a molecule energized high enough to surmount the potential barrier. This process is schematically represented as

$${\left[\mathrm{MA}\right]}^{*}\to \mathrm{M}+\mathrm{A},$$

(16)

where an energized molecule ${\left[\mathrm{MA}\right]}^{*}$ may have been produced by collision and/or photo-excitation, and so on. From the view point of quantum scattering theory, ${\left[\mathrm{MA}\right]}^{*}$ can be regarded as a core (temporally) excited resonance, and the unimolecular dissociation may be viewed as a half-collision process of the Feshbach resonance [49,50]. In case where those decaying vibrational energy levels are densely packed having high density of states, a state in the vibrational-state manifold wanders among many nonlinearly coupled vibrational states with a long life-time. It is intuitively clear that statistical prediction should work well in such a situation, since the corresponding microcanonical ensemble is expected to be highly mixing. The years of 1970s have been devoted to prove that the ergodicity and mixing in the relevant classical phase space of ${\left[\mathrm{MA}\right]}^{*}$ is indeed realized by classical chaos. Such stochastic energy redistribution within a molecule is now generally termed Intramolecular Vibrational energy Redistribution (IVR), which represents one of the most important concepts in chemical dynamics and will be briefly touched upon below. In short, statistical mechanics works surprisingly well even in such small molecules composed only of several atoms, not of the Avogadro number. Those dynamics deviating from the RRKM theory are now placed under special focus for their possible peculiar phase-space structures, which are often termed as non-RRKM dynamics.

The key quantity that brings about stochastic and statistical nature into small molecular dynamical systems is the density of state (DOS), not the number of particles (atoms and molecules). In case where the DOS is low, the quantum dynamics should be subject to the Schrödinger equation and naturally deterministic. The higher is DOS, the harder become the relevant coupled equations to solve, and statistical treatments turn out to be more feasible. Due to the heavy masses of nuclei, the vibrational and rotational degrees of freedom of a molecule tend to have very high DOS, which can readily amount to the order 10${}^{10}$ eV${}^{-1}$ even for several atomic systems. This is not as large as the Avogadro number but sufficiently large to allow for statistical treatment. On the other hand, the DOS for the typical low-lying electronic states is less than 1.0 eV${}^{-1}$, and quantum mechanics is almost only one faithful approach. For this reason, the theory of stationary molecular electronic states is widely called quantum chemistry. For these low DOS and discrete states, it is a natural assertion that there is no chaos. Yet, in the end of this review, we will penetrate the electronic state manifold of densely quasi-degenerate states, ending up with finding chaos.

3.2. Intramolecular Vibrational Energy Redistribution (IVR)

Given the statistical assumption of molecular vibrational energy distribution, the RRKM theory actually estimates the rate of unimolecular dissociation based on statistical inference at the bottle neck position on the way to products in phase space [48], using the transition state [40,41,42]. Yet, chemical dynamics has made a further step to the physical processes of intramolecular vibrational (along with rotational) energy flow (redistribution) within a molecule.

Such experimental studies were made possible by the advent of short pulse laser, which has unified two then separated subfields in molecular science, namely, spectroscopy of molecular structures and chemical reaction dynamics. As such the so-called stimulated emission pumping (SEP) spectroscopy applied to acetylene molecule H-C≡C-H (linear in shape having a triple bond) has given a huge impact [51,52,53,54]. Just briefly is outlined the essence of the experiment below. Preparing beforehand an electronically excited state of H-Ċ=Ċ-H (bent structure and longer C-C distances, double bond), they de-excited it to the high vibrational excited state on the electronic ground state by induced emission. Varying the time-width of the pumping laser, they observed very interesting vibrational emission spectra in the various hierarchical levels of resolution. Reflecting the highly vibrational excited states of polyatomic molecule, the attained spectra were extremely complicated, manifesting the chaotic feature behind. Surprisingly nevertheless, Yamanouchi succeeded to assign the spectral origins of most of the relevant peaks [52]. Among others, the most striking spectral feature is the so-called nested clump structure: Each of clump-like spectral peaks of a rather broad width is found to be composed of finer peaks attained when a longer width pulse is applied. Each of these fine peaks in turn shows another clump structure and is observed to consist of the even finer peaks with the pulse of further long width. They ascribed the nested structures to the stepwise intramolecular energy redistribution triggered by the mode change from the vibrational structure of H-Ċ=Ċ-H to that H-C≡C-H. For example, the elongated Ċ=Ċ bond releases much vibrational energy to the shorter C≡C in the first stage, which are transferred to the C-H stretching modes, which in turn coupled with the bending motion of ∠HCC in due time ordering. More complicated energy redistribution can follow such as a very large amplitude motion of the hydrogen atoms, but we do not pursue them here. Since the spectra of the longer time-width can track the long physical process of energy redistribution. The related experiments and precise analyses have vividly shown the real-time and actual energy flow within a molecule. The notion of chaos has thus physical reality in chemical dynamics.

The energy transfer among the different vibrational modes are induced by high anharmonicity of the relevant potential functions, leading to nonlinear canonical equations, and are also due to the overlap among the Fermi resonances [55] or the Feshbach resonance [50]. (see ref. [5] for the overlap criterion of chaos.) Therefore, the relevant dynamics can readily result in chaos leading to ergodicity [56]. Many classical mechanical studies have been devoted to IVR. Yet, only recently become available extensive reviews on quantum ergodicity and intramolecular energy flow due to Leitner [57] and quantum ergodicity transition in relation to IVR in a phase space perspective by Karmakar and Keshavemurthy [58].

Incidentally, energy transfer, relaxation, and redistribution of intra- and inter-molecular systems in condensed molecular systems [55,59] and biological molecules [60,61] are also among the basic processes of many of chemical phenomena, irrespective of the extent of possible contribution by quantum chaos. Currently, molecular dynamics studies are among the main streams in these fields, but those are beyond the scope of this review.

3.3. Onset of Statistical Properties: Liquid-Like Clusters

Beck and Steve Berry are the first, to the best of out knowledge, who found the critical roles of classical chaos in the statistical behavior of “melting and phase transition” appearing in structural isomerization dynamics of small atomic clusters [62]. Berry and their coworkers have achieved many fundamental works about cluster dynamics. See ref. [63] as a representative paper. We also refer to the work by Chandler et al. [64] as a pioneering study on stochastic dynamics of molecular isomerization.

Chaos not only makes the relevant phase space ergodic and mixing but also gives interesting qualitative features to the dynamics of molecules. Such properties can be typically seen in the structural isomerization dynamics of atomic clusters. Take an example from Ar${}_7$ cluster on the Lennard-Jones potential, which has 4 isomers as depicted in Figure 2. (Permutation symmetry counts more than 5000 isomers.) Each isomer is supported by its own potential energy basin. In a low energy, each isomer can stay within its potential basin keeping its structure as a solid, while at some high energy the isomerization begins to take place by surmounting the potential barriers, and further in considerably high energy frequent isomerization keeps taking place, which can be viewed as extensive large scale vibrational motions over the entire molecule or as a liquid-phase analog [63]. Even higher energy induces evaporation of one atom and two atomic molecule [65,66]. This entire process is hence regarded as a molecular level prototype of the first-order phase transition [62,63,67].

In the liquid-like phase, chaotic dynamics manifest themselves in the following phenomena. (1) Each classical path wanders among the potential basins for individual isomers in a random fashion (see Figure 2). The time-series with respect to the names of the isomers such as {PBP→COCT→IST→COCT→PBOP→SKEW$\to \dots$} is stochastic in the Newtonian dynamics, yet a path stays for somewhat long time in each basin [68,69]. (2) Memory loss of the classical paths occurs. Despite mutually different potential barriers among the individual potential basins, they behave as though they forget which basin they came from and to which basin they are going to visit next. They seem not to care about the presence of the so-called transition state, which is the lowest energy saddle between two basins [68,69]. (3) The time scale for a trajectory to penetrate into a mixing area after getting into a potential basin is so short that trajectories readily lose the memory in the mixing area in phase space. We call this situation inter-basin mixing [70]. (4) Microcanonical temperature can be defined with which the lifetimes of the isomers are subject to an Arrhenius relation. This in turn suggests the existence of the free energy despite the nonequilibrium dynamics [71].

Behind the above chaotic dynamics we can imagine geometrical structures in phase space, referred to as the reaction tubes. Suppose that a trajectory enters a PBP basin from the neighboring COCT by passing through a divided surface (many dimensional). Then there should exist a large bundle of similar trajectories flowing from the COCT to this PBP that includes this particular trajectory. See Figure 3, which illustrates a part of the unit piece of reaction tube. Into this PBP basin schematically represented as a circle, for instance, many other reaction tubes come in from other neighboring potential basins such as those of SKEW and IST. Soon after the tubes come in, each branches frequently as schematically shown in the figure. The successive branching actually reduces the diameter exponentially as time. Indeed one can measure the exponent of the decay in the form $exp\left(-\beta t\right)$ [70,72]. This is a geometrical background for classical paths to lose their memory from which basin they came in. Those extremely thin tubes entangle with other tubes of different histories, and they repeatedly merge with one another to thicker tubes and each eventually gets out of this basin to the neighboring ones (take the time reversal dynamics of the branching). This is how they lose the future “memory” of which basin they are going to visit next [70,72]. By construction it is obvious that the initial fat trunk of the reaction tube depicted in Figure 3 is already a mixture of finer tubes of different histories, each of those finer tubes being composed of even finer ones. Thus, the cross-section of the tube should exhibit very complicated fractal structures, along with globally periodic orbits of many different symbolic dynamics and different periods. In summary, repeated branching and remerging of the reaction tubes exist behind the memory loss.

The original idea of reaction tube is due to De Leon and their colleague, which they called the cylindrical manifold in the study of two dimensional two potential basin dynamics [73,74]. “Bifurcation and remerging (reconnection)” of the reaction tubes can be regarded as a basic mechanism of complex isomerization dynamics besides the baker’s transformation. There are interesting studies about the fine structure of the reaction tubes relevant to the stochastic dynamics. However we skip this aspect to quantum chaos.

4. Quantum Dynamics in the Quasi-Separatrix of Two Dimensional Molecular Vibration

In the early 1980s the studies of quantum chaos in terms of quantum wavepacket propagation were made rather extensively in the community of chemical physics [75,76,77,78,]. In particular, Feit and Fleck [78] are among the first, to the best of our knowledge, who numerically integrated the Schrödinger equation for wavepackets on the Hénon-Heiles Hamiltonian [27]. The scientific insights attained thus far with the accurate and the efficient methodologies developed for those studies are indispensable to date.

We herein study pure quantum chaos in terms of wavepackets and eigenfunctions confined in the separatrix on the modified Hénon-Heiles potential. This potential function is regarded as a two dimensional potential energy surface (PES) on which for molecular vibrational motion to develop.

4.1. Quantum Wavepackets Embedded in the Quasi-Separatrix

The Hamiltonian we adopt is

$$H=\frac{1}{2m_x}{p}_x^2+\frac{1}{2m_y}{p}_y^2+\frac{1}{2}(x^2+y^2)+x^2\left(a{y}^2+y\right)+\frac{1}{3}{y}^3(by-1),$$

(17)

where the masses are chosen as $m_x=1.0087$ and $m_y=1.0$ and $a=0.6$, $b=0.2$with the Planck constant $\hslash =0.005$. The parameters have been set so as to close the dissociation channel inherent to the original Hénon-Heiles potential. The potential functions are drawn in the upper panels of Figure 4, while the Poincaré surface of section at $x=0$ with $p_x>0$ for a selected energy is drawn at the left of the lower panel. One of the regular vibrational modes lying in a torus, marked with “lib”, is drawn in the $\left(x,y\right)$-space right to the Poincaré surface of section. An appropriate energy is chosen so that the quasi-separatrix has a small to medium size so as to coexist with the set of tori of significant scales. These altogether represent weak chaos.

We first track a wavepacket propagating in time within the major quasi-separatrix. To do so, we extended the Poincaré surface of section to a quantum version as follows [79,80]. Given a two-dimensional wavefunction $\psi \left(x,y,t\right),$ we first take a momentum representation only in the x-coordinate such that

$$\tilde{\psi }\left(p_x,y,t\right)={\int }_{-\infty }^{\infty }dxexp\left(-\frac{i}{\hslash }{p}_{x}x\right)\psi \left(x,y,t\right).$$

(18)

Then Fourier back it using only the positive $p_x$ such that

$${\psi }^{+}\left(x,y,t\right)={\int }_0^{\infty }d{p}_{x}exp\left(\frac{i}{\hslash }{p}_{x}x\right)\tilde{\psi }\left(p_x,y,t\right)$$

(19)

and put $x=0$ there for making a cross section at $x=0$. Bring it then into the Husimi representation in such a way

$$\begin{array}{cc}\hfill {\rho }_H\left(q_y,p_y,t\right)& ={(2\pi \hslash )}^{-1}{(\pi \hslash )}^{-1/4}\hfill \\ & \times {\left|{\int }_{-\infty }^{\infty }dyexp\left(-\frac{1}{2\hslash }{\left(y-p_y\right)}^2+\frac{i}{\hslash }{p}_{y}y\right){\psi }^{+}(0,y,t)\right|}^2\hfill \end{array}$$

(20)

We then plot ${\rho }_H\left(q_y,p_y,t\right)$ in the $\left(q_y,p_y\right)$-plane corresponding to the classical counterpart. This is rather a faithful extension of the classical Poincaré section.

An example of the quantum Poincaré section is exhibited in Figure 5, the upper row indicates ${\left|\psi \left(x,y,t\right)\right|}^2$ at selected times, while the lower row shows the corresponding quantum Poincaré section. It is clear in the quantum Poincaré section that the present packet is running on the quasi-separatrix. It is noteworthy that in upper panels of Figure 5a,c, the wavepacket stays around the unstable fixed points for a longer time, thereby leaving two significant strips in the wavefunction in the corresponding ($x,y$) space. These dense strips are an example of the so-called quantum scar first found by Heller [80] (see also E. J. Heller, page 547 in ref. [4]).

4.2. Energy Spectra and Eigenfunctions in the Quasi-Sparatrix

The next question is whether the eigenfunctions can be built in the quasi-separatrix without decaying to the stable tori. We next explore this aspect.

4.2.1. Energy Screening of Quantum Wavepackets

Prior to the numerical presentation we need to show how only (or at least mostly) the separatrix eigenfunctions are extracted from the very many eigenfunctions. This is a difficult task if we resort to the straightforward method of diagonalization of the Hamiltonian matrix, which may be expanded in arbitrary basis functions. Fortunately, those wavepackets we have shown above in Figure 5 are certainly running on the quasi-separatrix, and it would be an appropriate idea to quantize those packets. Our way to do so was as follows [81].

As a standard way to obtain eigenfunctions at any chosen energy E from a wave packet, one may want to calculate the time-energy Fourier transform [83,84]

$$\psi (x,E)={(2\pi \hslash )}^{-1}{\int }_{-\infty }^{+\infty }{e}^{+\frac{i}{\hslash }Et}\psi (x,t)dt.$$

(21)

However, this approach needs to pick extremely accurate $E^{\prime }s$ beforehand, usually by means of the Fast Fourier Transform of a correlation functions. (Incidentally see ref. [85] to extract very accurate spectrum from the FFT beyond a given grid size.) This is because the integral in Equation (21) simply vanishes in use of an inaccurate E. Moreover one may need to know in advance whether E should indeed belong to one of the quasi-separatrix eigenfunctions.

An efficient alternative to calculate eigenfunctions directly without the prior knowledge on the energy spectrum is the screening method [81]. Suppose we seek for an eigenfunction in an arbitrary narrow energy range $[{ϵ}_1,{ϵ}_2]$. A wavepacket residing in the interval can be generated as

$$\begin{array}{cc}\hfill \left|\Phi \left[{ϵ}_1,{ϵ}_2\right]\right.〉& ={\int }_{{ϵ}_1}^{{ϵ}_2}dE\delta (E-\mathcal{H})|\varphi \left(0\right)\rangle \hfill \\ & =\frac{1}{2\pi }{\int }_{-\infty }^{\infty }dt\frac{e^{i{ϵ}_{2}t}-e^{i{ϵ}_{1}t}}{it}|\varphi \left(t\right)\rangle ,\hfill \end{array}$$

(22)

with the help of

$$\delta (E-\widehat{H})={(2\pi \hslash )}^{-1}{\int }_{-\infty }^{+\infty }dt{e}^{\frac{i}{\hslash }(E-\widehat{H})t},$$

(23)

where $|\varphi \left(t\right)\rangle =e^{-i\mathcal{H}t}|\varphi \left(0\right)\rangle ,$ and $\left|\varphi \right(0)\rangle$ is an initial wavepacket supposed to run on the quasi-separatrix in the present purpose. If there is only one eigenstate in $[{ϵ}_1,{ϵ}_2]$, this wavefunction $\left|\Phi \right[{ϵ}_1,{ϵ}_2]\rangle$ should converge to to the exact one, aside from the norm. Its energy, if needed, is to be readily estimated as an expectation value. If the target energy range does not contain any eigenfunction, $\left|\Phi \right[{ϵ}_1,{ϵ}_2]\rangle$ converges to the null. When several eigenstates are included in the range, one can step to a further elaboration, which is not difficult at all [82,86]. Incidentally, notice the faster convergence in Equation (22) than in Equation (22) due to the presence of t in the denominator.

4.2.2. Chaotic Eigenfunctions and Spectra

Figure 6 represents an almost degenerate pair of eigenfunctions belonging to the quasi-separatrix wavefunctions. The left (right) panel represents an odd (even) function with respect to the mirror line at the $x=0$, with energy $E=0.152094$ ($E=0.152500$). The two energies are close to the one for the classical Poincaré section of Figure 4. As seen in the nodal structures, these are actually highly excited states. The classical motions at the two major unstable fixed pints are separated by symmetry as a left-handed (right-handed) librating motions (suggested by panel (a) and (c) in the upper row of Figure 5). However, the quantum wavepackets mix together by floating from and into those fixed point areas. Looking back at the structure of the quasi-separatrix in Figure 4, it is noticed that two wide areas around the fixed points (chaotic seas) are connected through thin canals, thus forming constrictions (see at $p_y=0$ in the section for instance). Therefore the wavepackets must pass through those thin channels to move from one wide chaotic sea to the other to become a global eigenfunction. In other words, those canals serve as “dynamical barriers”, which suppress the mobility of the wavepackets. The above pair of eigenfunctions in Figure 6 are generated through such dynamical barriers, and the slight difference in their energies and the odd-even symmetry separation are all the consequences of that kind of blockade or quasi-tunneling.

Davis and Heller found a quantum mechanical tunneling between two (or more) tori, which are separated by separatrix(es). This phenomenon is referred to as the dynamical tunneling through dynamical barrier [87,88,89]. In contrast, the pair of eigenstates observed in Figure 6 are not literally dynamical tunneling, since the canal-passing is allowed both quantum- and classical-mechanically. However, the view of the Poincaré section in Figure 5 suggests that the chaotic seas around the major fixed points reserve two major wavepackets to stay long, which are weakly separated by the thin canals thereby serving as dynamical barriers effectively. We therefore referred to this tunneling-like motion as dynamical tunneling of the second kind [80]. This phenomenon seems to be peculiar to quantum chaos, and yet we are not aware whether it has been found in actual molecular spectra.

4.2.3. Dependence on the Magnitude of the Planck Constant

We are curious about the possible effects of the magnitude of the Planck constant on quantum chaos. To save the space, we refer this aspect to our original paper [79]. Yet we hereby briefly touch upon the smoothing of chaos in the larger value of ℏ as demonstrated in terms of the Poincaré section. The above two panels in Figure 7 come from a separatrix eigenfunction with $\hslash =0.005$ ($[E_1,E_2]=[0.150.0.153]$), while the lower panel displays the section $\hslash =0.010$ ($[E_1,E_2]=[0.150.0.156]$). Such quantum smoothing is not at all surprising, but the corresponding energy spectra are found not as simple (not shown here) chematic [78]. As a chemist, I suspect that the biological worlds on the globe may not be robust against a small change of the Planck constant.

4.3. Time-Dependent Spectrum of Very Long-Time Dynamics

The classical trajectories within the thin quasi-separatrix are generally not periodic nor quasi-periodic by definition (meaning a zero measure for the exception), and moreover they wander around circulating near one tori to another one from time to time, each torus of which represents different mode of vibration (see Figure 4). Therefore those trajectories seem to continue to undergo different vibrational motions mostly in an intermittent way [90]. Therefore it can take extremely long time to extract the associated vibrational quantum spectra from the correlation functions based on those trajectories. The resultant spectrum can be time dependent, meaning that it is dependent on which timing we begin to calculate the correlation function and on when we finish. Such an example of time-dependent spectra are seen in [79,80] This is a serious issue in vibrational spectroscopy for complicated molecules of chaotic behavior. (Recall the nested clump structure in the SEP spectrum, mentioned in Section 3.2) The structural isomerization dynamics in the liquid-like phase discussed in the previous section (Section 3.3) is actually a large amplitude vibrational motion of the extremely long time-scale associated with chaos. It is still left unresolved how to take vibrational spectra for such wandering motions from one potential basin to another, spending some long time at each basin. Proteins in molten globule state must have the similar nature to the liquid-like clusters.

5. Quantization of Chaos with Semiclassical Wavepackets

The mechanism of energy quantization is a landmark lying in between classical and quantum mechanics. It is therefore quite natural for semiclassical mechanics to play a major role there. Besides, the chemical dynamics quite often demands semiclassical treatment for the nuclear motion apart from chaos. There are many beautiful studies in the literature on the EBK semiclassical quantization for regular systems (see ref. [91,92] for representative examples). As for quantization of chaotic dynamics the periodic orbit theory, or the trace formula, due to Gutzwiller [35,36,37] is the central work [2,4,5,6,7,8,9,10,12]. One of the most intriguing aspects of the trace formula is its relation to the Riemann and Selberg zeta functions. However, reviewing it is by far beyond the scope of this article [93].

5.1. The Gutzwiller Periodic Orbit Theory

The periodic orbit theory [35,36,37] treats the semiclassical estimate of the density of states with the Feynman kernel K as

$$D\left(E\right)=\int dq\int dtK\left(q,q,t\right)exp\left(\frac{i}{\hslash }Et\right).$$

(24)

The theory applies the (first-order) semiclassical approximation to the Feynman kernel or its energy representation (the energy Green function) to the density of states, thereby extracting the essential roles of phase-space periodic orbits [94]. The successive applications of the stationary phase approximations give rise to the famous expression

$$D\left(E\right)\simeq \frac{1}{{\left(2\pi \hslash \right)}^N}\int \int dqdp\delta \left(E-H_{cl}\left(q,p\right)\right)+\frac{1}{\pi \hslash }\sum _{\gamma }\sum _{k=1}^{\infty }\frac{T_{\gamma }cos\left(k\left(\frac{W_{\gamma }}{\hslash }-\frac{\pi }{2}{\lambda }_{\gamma }\right)\right)}{\sqrt{\left|Det\left(M_{\gamma }^k-I\right)\right|}},$$

(25)

where the first term of the right hand side is the Thomas–Fermi (classical) density of states with the classical Hamiltonian $H_{cl}\left(q,p\right)$ arising from the periodic orbits of zero length and the second oscillatory terms represent the quantum interference among the finite periodic orbits (specified by $\gamma$) in terms of the period $T_{\gamma }$, the action integral $W_{\gamma }$, and the Monodromy matrix $M_{\gamma }^k$ in the transversal direction of phase-space periodic motion. These terms should be summed up with respect to the number of circulations k over all the possible periodic orbits. However, it is very difficult particularly in many-dimensional systems to find the necessary periodic orbits (see ref. [95] for instance) and count up their contributions accurately. Moreover, it is well known that the periodic orbits proliferate exponentially as the length of the period $T_{\gamma }$. Another intrinsic difficulty in the path summation appears in the denominator in the second term of $D\left(E\right)$. This is a direct reflection that most of the semiclassical expressions bear the serious theoretical difficulty related to divergence of the monodromy matrix, which explodes at the caustics or the turning points. The divergence is particularly prominent in chaotic systems. It is therefore often claimed that the series of periodic orbit sum in Equation (25) may not be absolutely convergent, which sometimes leaves the trace formula to be obscure in practice. Furthermore, a naive question may be cast at whether the oscillatory terms in the trace formula can indeed cancel the Thomas–Fermi density at any energy of not-to-be quantized (off-resonant energy).

It seems that there are two major opposite directions to cope with the fundamental difficulty of the trace formula. One is an attempt to recount the periodic orbit contributions in Equation (25) in such a way to attain a conditional convergence [96]. Furthermore, the convergence of the trace formula can be better considered in the complex energy plane [97]. Semiclassical mechanics of higher order ℏ is obviously a canonical candidate to improve the trace formula [98]. Deep mathematics are involved in these approaches.

The other one is to take an easier, simpler and more tractable approach, mainly based on the semiclassical propagation of wavepackets. Our approach below is along this line. Besides, semiclassical wavepacket theory is one of the main methodologies in chemical reaction dynamics. As such, we here present our own approach, which is intended to step further beyond semiclassical approximation to the exact one. We then attain some insights about the mechanism of energy quantization, that is, the explicit roles of not only constructive but destructive interferences in forming discrete energy levels. We also see that the energy levels in chaos can be determined without the tedious pre-exponential factor (or the Monodromy matrix) of the semiclassical kernel.

5.2. Wavepacket Semiclassics with Action Decomposed Function

Many semiclassical theories have been proposed in the literature [7,14,15,34,94,99,100,101,102,103,104], either starting from the Feynman kernel, WKB theory or others. Yet, we consider the following theory to clarify the dynamical hierarchical structure in the range extending from classical and semiclassical to full-quantum dynamics.

5.2.1. Short Time Dynamics of the Maslov Type Wavefunction

We represent a wavefunction $\Psi \left(q,t\right)$ in the so-called Maslov type [102], such that

$$\Psi \left(q,t\right)=F(q,t)exp\left(\frac{i}{\hslash }{S}_{cl}\left(q,t\right)\right),$$

(26)

for a short time interval $\left[t,t+\Delta t\right]$ on a multidimensional coordinate q in configuration space, where $S_{cl}$ is assumed to satisfy the classical Hamilton–Jacobi (HJ) equation

$$\frac{\partial S_{cl}}{\partial t}+\frac{1}{2m}{\left(\nabla S_{cl}\right)}^2+V=0,$$

(27)

which defines relevant classical paths during the interval, with m being the mass. The classical factors are thus factored out for this short interval and thereby the purely quantum factors are thus supposed to be included only in the function $F(q,t)$[32,33]. $F(q,t)$ is thus referred to as Action Decomposed Function (ADF) [29,30]. Then in this interval, the Schrödinger equation for $\Psi \left(q,t\right)$ is transformed to a linear equation of motion for the complex valued amplitude function $F(q,t)$ as

$$\frac{\partial F(q,t)}{\partial t}=\left(-p·\nabla -\frac{1}{2}(\nabla ·p)\right)F(q,t)+\frac{i\hslash }{2}{\nabla }^{2}F(q,t),$$

(28)

where p is a momentum at $\left(q,t\right),$ which is $q\left(t\right)$, as

$$p=\nabla S_{cl}\left(q,t\right).$$

(29)

We use the mass weighted coordinates so that all the masses (m) are scaled to unity $m=1$, and thereby p is set to be numerically equivalent to the corresponding velocity v. In transforming the Euler picture to the Lagrange by defining

$$\frac{D}{Dt}=\frac{\partial }{\partial t}+v·\nabla$$

(30)

Equation (28) reads

$$\frac{D}{Dt}F(q,t)=\left[-\frac{1}{2}(\nabla ·p)+\frac{i\hslash }{2}{\nabla }^2\right]F(q,t).$$

(31)

$F(q-q_{cl}\left(t\right),t)$ is hence carried along a classical path $q_{cl}\left(t\right)$ in this interval. Then the Trotter decomposition [94] valid for a very short time–interval $\Delta t$ gives

$$\begin{array}{c}F\left(q-q_{cl}(t+\Delta t),t+\Delta t\right)\hfill \\ \simeq exp\left[\frac{i\hslash }{2}\Delta t{\nabla }^2\right]exp\left[-\frac{1}{2}(\nabla ·p)\Delta t\right]F\left(q-q_{cl}\left(t\right),t\right)\hfill \\ \simeq exp\left[-\frac{1}{2}(\nabla ·p)\Delta t\right]exp\left[\frac{i\hslash }{2}\Delta t{\nabla }^2\right]F\left(q-q_{cl}\left(t\right),t\right)\hfill \end{array}$$

(32)

and represents the quantum effects in a product between the momentum gradient ($exp\left[-\frac{1}{2}(\nabla ·p)\Delta t\right]$) and the pure quantum diffusion ($exp\left[\frac{i\hslash }{2}\Delta t{\nabla }^2\right]$). The theory is exact except for the Trotter decomposition, provided that the classical action $S_{cl}\left(q,t\right)$ exists.

5.2.2. Semiclassical Approximation

Equation (32) is first reduced to a semiclassical level of approximation by considering the momentum gradient part only, which results in

$$\begin{array}{c}F\left(q-q_{cl}(t+\Delta t),t+\Delta t\right)\hfill \\ \simeq exp\left[-\frac{1}{2}(\nabla ·p)\Delta t\right]F\left(q-q_{cl}\left(t\right),t\right)\hfill \end{array}$$

(33)

Notice that the Planck constant does not yet appear in this level. The time evolution of it is approximately given by estimating the momentum gradient term as follows [29]

$$\begin{array}{c}{\left.F(q-q(t+\Delta t),t+\Delta t)\right|}_{q=q_{cl}(t+\Delta t)}\hfill \\ ={\left.exp\left[-\frac{1}{2}{\int }_t^{t+\Delta t}dt\sum _k\frac{\partial p_k}{\partial q_k}\right]F\left(q-q_{cl}\left(t\right),t\right)\right|}_{q=q_{cl}\left(t\right)}\hfill \\ ={\left.{\left(\frac{\sigma \left(t\right)}{\sigma (t+\Delta t)}\right)}^{1/2}F\left(q-q_{cl}\left(t\right),t\right)\right|}_{q=q_{cl}\left(t\right)}\hfill \end{array}$$

(34)

or

$${\left.\sigma {\left(t+\Delta t\right)}^{1/2}F(q-q_{cl}\left(t+\Delta t\right),t+\Delta t)\right|}_{q=q_{cl}\left(t+\Delta t\right)}={\left.\sigma {\left(t\right)}^{1/2}F(q-q_{cl}\left(t\right),t)\right|}_{q=q_{cl}\left(t\right)},$$

(35)

where the deviation determinant $\sigma \left(t\right)$ is expressed as

$$\sigma \left(t\right)=\prod _{i=1}^N\wedge \left(q_{cl}^i\left(t\right)-q_{cl}\left(t\right)\right),$$

(36)

which is an N-dimensional orientable tiny volume surrounding the point $q_{cl}\left(t\right)$ in configuration space, with $q_{cl}^i\left(t\right)$ ($i=1,\dots ,N$) being the configuration space point of an ith classical trajectories nearby the reference path $q_{cl}\left(t\right)$. See Figure 8 for geometrical definition of $\sigma \left(t\right)$. Since $\sigma \left(t\right)$ arises from the orientable volume, the square root of it gives rise to the part of the so-called Maslov phase [105]. (see ref. [106] for another geometrical meaning of the Maslov phase.)

The above results can be represented in a symmetric manner by defining “square root” of the volume element [107]. Suppose that an arbitrary function $f\left(q_t\right)$, where $q_t$ is the coordinate system in configuration space, is mapped from the original point $q_0$ at time 0 by the classical flow $q_0\to$$q_t$. An integral also mapped is represented as

$$I=\int f\left(q{}_t,t\right)d{q}_t=\int f\left(q{}_t,t\right)d{q}_t^{1/2}d{q}_t^{1/2*}$$

(37)

with

$$d{q}_t^{1/2}\equiv exp\left[\frac{i\pi M\left(q_0,q_t\right)}{2}\right]{\left|d{q}_t\right|}^{1/2},$$

(38)

where $d{q}_t^{1/2*}$ is intended to be the complex conjugate to $d{q}_t^{1/2}$, and $M\left(q_0,q_t\right)$ is the Maslov index counting the number of the change of sign of $\sigma \left(t\right)$ of Equation (36). The semiclassical conservation rule is then written up to the phase as [107]

$$F\left(q{}_t,t\right)d{q}_t^{1/2}=F\left(q{}_0,0\right)d{q}_0^{1/2}.$$

(39)

This is rather similar to the law of the Poincaré-Cartan absolute invariance in phase space [108].

The momentum gradient ${\sum }_k\partial p_k/\partial q_k$ can diverge, when it happens to hold $\sigma \left(t\right)=0$ in Equation (35), where a classical path come to a focal point either in configuration or momentum space. This is the origin of the difficulty of the theories of semiclassical level. Yet, such a singularity can be removed by the quantum diffusion below [29].

5.2.3. Gaussian Wavepacket Approximation

We next proceed to take account of the quantum diffusion term. In doing so, we here survey the dynamics of a Gaussian wavepacket among other possible choices of $F(q,t)$. Heller introduced the dynamics of a Gaussian function as early as in 1975 by assuming that the potential function at the Gaussian center is truncated up to the quadratic terms [109].

We insert a Gaussian of the following form into the general scheme of ADF Equation [29],

$$G(q-q_{cl}\left(t\right),t)=g\left(t\right)exp\left[-\frac{1}{c\left(t\right)+id\left(t\right)}{(q-q_{cl}\left(t\right))}^2\right],$$

(40)

with both $c\left(t\right)$ and $d\left(t\right)$ being real-valued. Then some (rigorous) manipulations give

$$c\left(t+\Delta t\right)=c\left(t\right){\left(\frac{\sigma \left(t+\Delta t\right)}{\sigma \left(t\right)}\right)}^2$$

(41)

and

$$d\left(t+\Delta t\right)=d\left(t\right){\left(\frac{\sigma \left(t+\Delta t\right)}{\sigma \left(t\right)}\right)}^2+\frac{2\hslash \Delta t}{m}.$$

(42)

or shows they satisfy the following differential equation, respectively,

$$\dot{c}\left(t\right)=2\frac{\dot{\sigma }\left(t\right)}{\sigma \left(t\right)}c\left(t\right).$$

(43)

and

$$\dot{d}\left(t\right)=2\frac{\dot{\sigma }\left(t\right)}{\sigma \left(t\right)}d\left(t\right)+\frac{2\hslash }{m}.$$

(44)

Note that the Planck constant appears only in this level of dynamical hierarchy. This normalized Gaussian form are propagated in time quite accurately, under the dynamical circumstances where the Gaussian may remain Gaussian. Yet, the theory needs additional procedures when the dynamics under study faces difficult situations like wavepacket bifurcation [110].

Further, one can find an interesting internal structure in the Gaussian packet. Define $\zeta \left(t\right)$ then such that

$$d\left(t\right)=\zeta \left(t\right)\sigma \left(t\right),$$

(45)

and insert this back into Equation (44) to obtain

$$\sigma \left(t\right)\dot{\zeta }\left(t\right)-\dot{\sigma }\left(t\right)\zeta \left(t\right)=\frac{2\hslash }{m},$$

(46)

which is a Wronskian relation between $\sigma \left(t\right)$ and $\zeta \left(t\right)$. This suggests that something is rotating to maintain the packet in the quantum realm [29]. However, we step aside from this subject to be back to the energy quantization problem [110].

5.3. On the Quantization Condition: Creation and Elimination of the Spectral Peaks by the Phases

It is quite often that only the to-be-quantized (resonant) energy peaks are under focus without paying much attention to to-be-nullified part of $D\left(E\right)$ at off-resonant energy areas. We therefore reconsider here the basic mechanisms of evolution of the resonant peaks and those of suppression of the off-resonant counterparts.

5.3.1. Prior Quantization Conditions Based on the Periodic Orbits

To reconsider the mechanism of quantization from the view point of wavepacket dynamics, we get back to a little more elementary stage of the theory. Consider a spectral density function

$$P\left(E\right)={\int }_{PO}\frac{d{q}_{0}d{p}_0}{{\hslash }^N}{P}_{tr}\left(E_{cl}(q_0,p_0);E\right),$$

(47)

for an N-dimensional system, where the pre-spectral density is defined as

$$\begin{array}{c}{P}_{tr}\left(E_{cl}\left(q_0,p_0\right);E\right)\hfill \\ =\underset{L\to \infty }{lim}\frac{1}{L}{\int }_0^{L}dtR\left(q_t,q_0,t\right)exp\left[\frac{i}{\hslash }S\left(q_t,q_0,t\right)-i\frac{\pi }{2}M\left(q_0,q_t\right)\right]exp\left(\frac{i}{\hslash }Et\right)\hfill \end{array}$$

(48)

which is a correlation function starting from a classical point $\left(q_0,p_0\right)$ on a manifold of energy $E_{cl}$. $M(q_0,q_t)$ is the Maslov index. $R(q_t,q_0,t)$ represents an overlap integral between a wavepacket at $t=0$ and that of later time t. We implicitly assume that the wavepacket adopted here has been taken from the above ADF framework, Equation (26), yet we do not specify the functional form of it in this stage and it does not have to be a Gaussian. The Fourier transform in $P_{tr}$ eventually gives rise to the energy spectrum in Equation (47). In the following analysis, we assume that (1) only the periodic orbits are to be considered in $P_{tr}\left(E_{cl}(q_0,p_0);E\right)$ as the symbol ${\int }_{PO}$ indicates, which is a natural consequence of the stationary-phase condition, (2) the phases in Equation (47) are just made up with the standard action integral and the associated Maslov phases alone, and (3) the amplitude factor $R(q_t,q_0,t)$ is always finite and smooth along the trajectories [111,112].

5.3.2. Peaks Arising from Individual Orbits and an Extended Quantization Condition

In the time-integration of $P\left(E\right)$ in Equation (48), one can apply the stationary phase approximation (SPA) to obtain

$$\frac{\partial }{\partial t}\left[S\left(q,q_0;t\right)+Et\right]=-E_{cl}+E=0,$$

(49)

which indicates that classical trajectories with $E_{cl}$ should make a major contribution to $P_{tr}\left(E_{cl}(q_0,p_0);E\right)$, which is already well-known [35,36,37,113,114]. And $P_{tr}\left(E_{cl}(q_0,p_0);E\right)$ can be viewed as a function of both E and $E_{cl}(q_0,p_0)$.

We first survey the oscillatory property of the integrand of $P_{tr}\left(E_{cl}(q_0,p_0);E\right)$. With the classical action integral

$$S(q_t,q_0,t)={\int }_{q_0}^{q_t}p\left(q;E_{cl}\left(q_0,p_0\right)\right)dq-E_{cl}\left(q_0,p_0\right)t,$$

(50)

we define

$$\Lambda \left(t;q_0,p_0\right)=R(q_t,q_0,t)exp\left[\frac{i}{\hslash }{\left.{\int }_{q_0}^{q_t}p\left(q;E_{cl}\left(q_0,p_0\right)\right)dq\right|}_{\mathrm{PO}}-i\frac{\pi }{2}M\left(q_0,q_t\right)\right],$$

(51)

to rewrite $P_{tr}\left(E_{cl}(q_0,p_0);E\right)$ in a simplified fashion as

$$P_{tr}\left(E_{cl}(q_0,p_0);E\right)=\int dt\Lambda \left(t;q_0,p_0\right)exp\left[\frac{i}{\hslash }\left(E-E_{cl}\left(q_0,p_0\right)\right)t\right].$$

(52)

Note that both ${\int }_{q_0}^{q_t}p\left(q;E_{cl}\left(q_0,p_0\right)\right)dq$ and $\frac{\pi }{2}M\left(q_0,q_t\right)$ have the following time dependence on the periodic orbits

$$\left(\mathrm{constant}\right)\times t+\left(\mathrm{oscillatory}\mathrm{part}\right).$$

(53)

Here the constant part is calculated out as an average of ${\int }_{q_0}^{q_t}p\left(q;E_{cl}\left(q_0,p_0\right)\right)dq$ and $\frac{\pi }{2}M\left(q_0,q_t\right)$ in $\Lambda \left(t;q_0,p_0\right)$ of Equation (51) for the interval of, say, $\left[q_0,q_T\right]$. Let the period of the periodic orbit under consideration be $T\left(q_0,p_0\right)$. Then we define a quantity $X\left(q_0,p_0\right)$ such that

$$\frac{1}{T\left(q_0,p_0\right)}\left[\frac{1}{\hslash }{\displaystyle \underset{q_0}{\overset{q_T}{\oint }}}p\left(q;E_{cl}\left(q_0,p_0\right)\right)dq-\frac{\pi }{2}M\left(q_0,q_T\right)\right]=X\left(q_0,p_0\right),$$

(54)

and extract $X\left(q_0,p_0\right)\times t$ from the phase part of (51) (that is Equation (53)) in such a way that

$$\frac{1}{\hslash }{\int }_{q_0}^{q_t}p\left(q;E_{cl}\left(q_0,p_0\right)\right)dq-\frac{\pi }{2}M\left(q_0,q_t\right)-X\left(q_0,p_0\right)t.$$

(55)

Then this function should appears to be an exactly periodic function. Furthermore, since the preexponential factor $R(q_t,q_0,t)$ is also assumed to be periodic of $T\left(q_0,p_0\right)$,

$$\Lambda \left(t;q_0,p_0\right)exp\left[-iX\left(q_0,p_0\right)t\right]$$

(56)

is periodic with the period $T\left(q_0,p_0\right)$ as a whole. Then one can expand the function $\Lambda \left(t;q_0,p_0\right)$ in a Fourier series such that

$$\Lambda \left(t;q_0,p_0\right)=\sum _{n=-\infty }^{\infty }{c}_{n}exp\left(-i\frac{2n\pi }{T\left(q_0,p_0\right)}t+iX\left(q_0,p_0\right)t\right),$$

(57)

where

$$c_n=\frac{1}{T\left(q_0,p_0\right)}{\int }_0^{T}dt\Lambda \left(t;q_0,p_0\right)exp\left[i\frac{2n\pi }{T\left(q_0,p_0\right)}t-iX\left(q_0,p_0\right)t\right].$$

(58)

Here we take a limit in the time integral for $P_{tr}\left(E_{cl}(q_0,p_0);E\right)$ in Equation (48) up to sufficiently long time L and rewrite it with the help of Equation (57) as

$$\begin{array}{c}{P}_{tr}\left(E_{cl}\left(q_0,p_0\right);E\right)\hfill \\ \equiv \frac{1}{L}{\int }_0^{L}dt\Lambda \left(t;q_0,p_0\right)exp\left[\frac{i}{\hslash }\left(E-E_{cl}\left(q_0,p_0\right)\right)t\right]\hfill \\ =\frac{1}{L}\sum _{-\infty }^{\infty }{c}_n\frac{-i\hslash T\left(q_0,p_0\right)}{-2n\pi \hslash +\left(E-E_{cl}\left(q_0,p_0\right)+X\hslash \right)T\left(q_0,p_0\right)}\hfill \\ \times \left(exp\left[\frac{iL}{\hslash T\left(q_0,p_0\right)}\left(-2n\pi \hslash +\left(E-E_{cl}\left(q_0,p_0\right)+X\hslash \right)T\left(q_0,p_0\right)\right)\right]-1\right)\hfill \end{array}$$

(59)

The correlation function $R(q_t,q_0,t)$ is assumed to be so smoothly varying (or much less oscillatory than the phase factor) that it can be practically factored out from the Fourier expansion with respect to the phases.

It is now clear from Equation (59) that only the denominator in this expression satisfying the condition

$$\left(E-E_{cl}\left(q_0,p_0\right)+X\hslash \right)T\left(q_0,p_0\right)=2n\pi \hslash$$

(60)

will lead to a peak in $P_{tr}\left(E_{cl}(q_0,p_0);E\right)$ of the height of $c_n$. The substitution of $X\left(q_0,p_0\right)$ with that of Equation (54) gives rise to a condition

$$\left(E-E_{cl}\left(q_0,p_0\right)\right)T\left(q_0,p_0\right)+\left[{\left.{\int }_{q_0}^{q_t}p\left(q;E_{cl}\left(q_0,p_0\right)\right)dq\right|}_{\mathrm{PO}}-\frac{\pi }{2}M\left(q_0,q_T\right)\hslash \right]=2n\pi \hslash .$$

(61)

This is a quantization (resonance) condition for $P_{tr}\left(E_{cl}(q_0,p_0);E\right)$ (but not for $P\left(E\right)$), which we call the prior quantization condition [111,115]. Not only E but $(q_0,p_0)$ has to fulfill this condition to make a peak in $P_{tr}\left(E_{cl}(q_0,p_0);E\right)$. Although $P_{tr}\left(E_{cl}(q_0,p_0);E\right)$ has far more peaks than $P\left(E\right)$ itself, many of these peaks are to be eliminated by destructive interference upon integration over the trajectory space, which will be discussed below.

Yet, the true peaks should be contained in the set of peaks fulfilling Equation (61). Note also that the additional condition of Equation (49) nullifies

$$\left(E-E_{cl}\left(q_0,p_0\right)\right)T\left(q_0,p_0\right)$$

(62)

in Equation (61), thereby reducing it directly to the EBK-like condition.

5.3.3. Destructive Interference Suppressing the Spurious Energy Peaks: Another Essential Role of the Phases for Quantization

The off-resonant trajectories should suppress the values of $P\left(E\right)$ to the order of O$\left(\hslash \right)$ due to the Riemann-Lebesgue lemma. However, there can exist many quasi-resonant trajectories, which would form spurious peaks in $P_{tr}\left(E_{cl}(q_0,p_0);E\right)$. These peaks are anticipated to be cancelled away by destructive interference with many other anonymous trajectories. As such there are two basic case in which to suppress virtual peaks as follows. (We refer to the original papers, ref. [111,112], for more details along with numerical presentations.)

Case in which $\mathbf{E}$ Is Slightly Shifted from a True Eigenvalue (Near-Resonance Peaks)

Suppose there is a peak satisfying

$$\left[{\int }_{q_0}^{q_t}p\left(q;E_{cl}\left(q_0,p_0\right)\right)dq-\frac{\pi }{2}M\left(q_0,q_T\right)\hslash \right]=2n\pi \hslash$$

(63)

at an energy $E_{cl}\left(q_0,p_0\right)=E_x$. It is expected theoretically that the contribution from nearby orbits having $E_{cl}^{\prime }(q_0^{\prime },p_0^{\prime })$ should make $P_{tr}\left(E_{cl}^{\prime }(q_0^{\prime },p_0^{\prime });E_x\right)$ very small to the order of O$\left(\hslash \right)$. However, if $E_{cl}^{\prime }(q_0^{\prime },p_0^{\prime })$ happens to be sufficiently close to the $E_x=E_{cl}\left(q_0,p_0\right)$, the peak is split into “a pair of very high spurious peaks” that contain the peak of $P_{tr}\left(E_{cl}^{\prime }(q_0^{\prime },p_0^{\prime });E_x\right)$ in between. These peaks constitute a continuous width near the true peak. (For numerical and graphical examples, see [111].)

Case of Generation of Harmonics in $\mathbf{E}$

There is yet another mechanism to generate spurious peaks. A trajectory with $E_{cl}(q_0,p_0)$ happens to satisfy the resonance condition of Equation (63) with some “quantum number” n and the energy $E_x$. Then, it is possible that the prior quantization condition of Equation (61) may be satisfied by another set of $\left(n^{\prime },E_x^{\prime }\right)$. They must be simply the harmonics of the correct peak of $\left(n,E_x\right)$ [111]. These harmonics also should be eliminated by off-resonant trajectories in terms of the appropriate destructive interference.

Having presented the specific roles of off-resonant trajectories to suppress the spurious peaks, we are concerned whether the trace formula implements such a function in it.

5.4. Amplitude-Free Energy Quantization of Classical Chaos

One of the theoretical and numerical difficulties of the trace formula manifests itself explicitly in the denominator of Equation (25) of the second term, which comes from the amplitude factor of the energy Greens function. In the above analysis of the energy spectrum, from Equations (48)–(61), no significant role was found for the amplitude factor $R(q_t,q_0,t)$ to play. We hence try to pursue a possible energy quantization disregarding the absolute value of the amplitude factor, but taking the Maslov phase into account.

To explicitly monitor the sensitivity of the amplitude factor to the variation of the preexponential factor, we first make the following preliminary calculations [116], in which the amplitude factor is rescaled intentionally as follows.

$${}^{}\begin{array}{cc}\hfill C_S^{ADF}\left(t\right)& =\int d{q}_0{F}^{*}\left(q_t,0\right)F\left(q_0,0\right){\left|\frac{\partial q_t}{\partial q_0}\right|}^{1-s}\hfill \\ & \times exp\left[\frac{i}{\hslash }{S}_1\left(q_t,q_0,t\right)+\frac{i}{\hslash }{p}_0·\left(q_0-q_t\right)-\frac{i\pi M}{2}\right]\hfill \end{array}$$

(64)

where $C_s^{ADF}\left(t\right)$ is a correlation function using a primitive Gaussian function, and an artificial scaling parameter s has been introduced. The value $s=0.5$ reduces $C_s^{ADF}$ to the original semiclassical correlation function, while there is no physical justification to choose other values. Yet, it numerically turns out that values up to at least $s=1.5$ give correct energy positions in the resultant spectra [116]. In particular the choice of $s=1.0$ implies that one can abandon the semiclassical amplitude factor as ${\left|\frac{\partial q_t}{\partial q_0}\right|}^0$, aside from the Maslov phase, to obtain the spectral positions.

After somewhat extensive analysis and numerical studies, we have confirmed that the spectrum positions are accurately reproduced through an amplitude correlation function [117,118], of the form

$$C_{p_0}^{III}\equiv \int d{q}_0{F}^{*}\left(q_t\right)F\left(q_0\right)exp\left(-\frac{i}{\hslash }{p}_0{q}_t\right)exp\left[\frac{i}{\hslash }{S}_2(q_t,p_0,t)-\frac{i\pi M}{2}\right],$$

(65)

where $p_0$ is an initial momentum to be scanned to cover the entire spectrum, and $F\left(q_t\right)$ is a (frozen) Gaussian function to be carried from the initial position $q_0$ [115,116]. We only need the way of calculation of the Maslov index as exemplified in Figure 8. It turns out that the amplitude-free correlation function $C_{p_0}^{III}$ works quite well even in chaotic regime, provided that the semiclassical dynamics under study is essentially dominated by the phases. An example of chaotic spectrum in chaotic region for the Hamiltonian

$$H=\frac{1}{2m_x}{p}_x^2+\frac{1}{2m_y}{p}_y^2+\frac{1}{2}(x^2+y^2)+x^2\left(0.6y^2+y\right)+\frac{1}{3}{y}^3(0.2y-1)+0.1x$$

(66)

with $m_x=1.0087$ and $m_y=1.0$ is shown in Figure 9, where we have adopted $\hslash =0.005$ and $E_{cl}\simeq 0.157.$ (Note that for a large ℏ, the quality of semiclassical quantization is accordingly deteriorated and should not be used.)

The partial success of the amplitude-free quantization without paying full attention to the semiclassical amplitude factors, except for the Maslov phase, suggests that one does not necessarily have to faithfully track the every fine complexity of the Lagrange manifold behind classical chaos in order to perform semiclassical mechanics, partly due to the quantum smoothing.

5.5. Chaos Mediated by Dynamical Tunneling

Tunneling effect is one of the key phenomena to characterize chemical reactions not only in very low temperature circumstances but even in excited states. It is now widely accepted that complex semiclassical theory is very effective in the study of quantum tunneling [119]. However, classical chaos generally makes the structures of the Lagrange manifold and associated caustics very complicated in multidimensional complex space. Besides, the choice of the stationary paths there is not necessarily achieved as in the Feynman kernel over the real paths. Shudo and Ikeda have developed a pioneering work to lay out and cope with the intricate yet beautiful tunneling paths in the complex domain [120,121]. We here present a simple-minded idea on semiclassical tunneling by introducing the parity of motion, which splits phase space into many mutually exclusive nonclassical sub-spaces [122,123].

5.5.1. Semiclassical Tunneling Theory

We first need to formulate a tunneling theory in a semiclassical context [122,123]. Let us begin with Feynman’s path integral representation of the kernel

$$K(x,\Delta t;y)={\left(\frac{m}{2\pi i\hslash \Delta t}\right)}^{N/2}exp\left(\frac{i\Delta t}{\hslash }\left[\frac{m}{2}{\left(\frac{x-y}{\Delta t}\right)}^2-V\left(\frac{x+y}{2}\right)\right]\right)$$

(67)

or in a more limited alternative

$$K(q+\Delta q,\Delta t;q)={\left(\frac{m}{2\pi i\hslash \Delta t}\right)}^{N/2}exp\left(\frac{i\Delta t}{\hslash }\left[\sum _k\frac{m_k}{2}{\left(\frac{\Delta q_k}{\Delta t}\right)}^2-V\left(q\right)\right]\right),$$

(68)

where $\Delta q$ is not necessarily as short as $\Delta t$, and $V\left(q\right)$ is an arbitrary potential function. In the path integration, the Euler variational principle

$$\frac{d}{dt}\left(\frac{\partial L}{\partial {\dot{q}}_k}\right)-\frac{\partial L}{\partial q_k}=0,$$

(69)

with the Lagrangian

$$\begin{array}{cc}\hfill L& =\sum _k\frac{m_k}{2}{\left(\frac{\Delta q_k}{\Delta t}\right)}^2-V\left(q\right)\hfill \\ & =\sum _k\frac{m_k}{2}{\left({\dot{q}}_k\right)}^2-V\left(q\right)\mathrm{as}\Delta q\to 0\mathrm{with}\Delta t\to 0\hfill \end{array}$$

(70)

brings about the stationarity in the integral [94]. It is seen though that the kernel of Equation (68) is invariant to a set of parameters (constant) ${\sigma }_k$ that are introduced in a manner that

$$K(q+\Delta q,\Delta t;q)={\left(\frac{m}{2\pi i\hslash \Delta t}\right)}^{N/2}exp\left(\frac{i\Delta t}{\hslash }\left[\sum _k\frac{m_k}{2}{\left({\sigma }_k\frac{\Delta q_k}{\Delta t}\right)}^2-V\left(q\right)\right]\right)$$

(71)

as long as

$${\sigma }_k^2=1.$$

(72)

Define accordingly the corresponding Lagrangian $L\left(\sigma \right)$ as

$$L\left(\sigma \right)=\sum _k\frac{m_k}{2}{\left({\sigma }_k{\dot{q}}_k\right)}^2-V\left(q\right)$$

(73)

which should satisfy the stationary condition

$$\frac{d}{dt}\left(\frac{\partial L\left(\sigma \right)}{\partial {\dot{q}}_k}\right)-\frac{\partial L\left(\sigma \right)}{\partial q_k}=0.$$

(74)

An additional set of the Newtonian equations then results

$$m_k{\ddot{q}}_k=-{\sigma }_k\frac{\partial V\left(q\right)}{\partial q_k}.$$

(75)

Nothing happens when all ${\sigma }_k=1$. If all the ${\sigma }_k$ is $-1$, the paths emerged from Equation (75) are equivalent to the so-called instanton path [94]. We refer to ${\sigma }_k$ as the parity of motion. A variety of the set $\left\{{\sigma }_k\right\}$ gives the various combinations of non-tunneling and tunneling paths in the real time. Thus, the sets of $\left\{{\sigma }_k\right\}$ defines the individual phase spaces, which we call sheets. No communication exists among them as long as the parity sets are frozen.

We can move to the Hamiltonian framework by defining

$$p_k=\sqrt{{\sigma }_k}\dot{q_k},$$

(76)

where $p_k$ is pure imaginary for the negative parity. The corresponding canonical equations of motion

$${\sigma }_k\dot{\overline{p_k}}=-\frac{\partial H\left(\sigma \right)}{\partial q_k}\mathrm{and}{\sigma }_k\dot{q_k}=\frac{\partial H\left(\sigma \right)}{\partial \overline{p_k}},$$

(77)

is defined with the modified Hamiltonian

$$H\left(\sigma \right)=\sum _k\frac{{\sigma }_k}{2m_k}{\overline{p_k}}^2+V\left(q\right).$$

(78)

Setting ${\sigma }_k=-1$ is effectively equivalent to flipping the mass $m_k$ from positive to negative value.

To see the invariance property inherent to the Hamilton system, the following coordinates are convenient, that is

$$Q_k=\sqrt{{\sigma }_k}{q}_k\mathrm{and}{P}_k=\sqrt{{\sigma }_k}{\overline{p}}_k,$$

(79)

in which the new coordinates $Q_k$’s in configuration space are also pure imaginary when their parities are negative. Then Equation (77) are transformed to

$${\dot{P}}_k=-\frac{\partial H\left(\sigma \right)}{\partial Q_k}\mathrm{and}{\dot{Q}}_k=\frac{\partial H\left(\sigma \right)}{\partial P_k}.$$

(80)

Due to the Poincaré-Cartan theorem of absolute invariance [108], the symplectic structure the 2-form

$${\omega }^2=\sum _{k}d{P}_k\wedge d{Q}_k$$

(81)

is conserved and the Stokes theorem allows to define the path-action $S_{cl}$, that is,

$$S_{cl}=\int \sum _k{P}_{k}d{Q}_k-H\left(Q,P\right)dt.$$

(82)

In order to adopt these paths in semiclassical mechanics, the Q-coordinates should be rotated back to q-coordinates keeping P in such a way that

$$\left(Q,P\right)\to \left(q,P\right).$$

(83)

Then we have an action integral

$$\begin{array}{cc}\hfill S_{cl}(Q,P)& \to S_{trl}(q,P)=\int \sum _{k}Pd{q}_k-H\left(\sigma \right)dt\hfill \\ & =\sum _k\int \sqrt{{\sigma }_k}{\overline{p}}_{k}d{q}_k-H\left(\sigma \right)dt\hfill \end{array}$$

(84)

and a semiclassical wavefunction for tunneling is formally written as

$$\psi =Aexp\left(\frac{i}{\hslash }{S}_{\mathrm{tnl}}\right),$$

(85)

where the pre-exponential factor A represents an amplitude. Since $S_{\mathrm{tnl}}$ in Equation (84), is generally complex-valued, we rewrite $S_{\mathrm{tnl}}$ in the standard form as

$$S_{\mathrm{tnl}}=S_R+i{S}_I.$$

(86)

Insertion of Equation (86) into Equation (85) gives

$$\psi =Aexp\left(-\frac{1}{\hslash }\left|S_I\right|\right)exp\left(\frac{i}{\hslash }{S}_R\right),$$

(87)

where the phase convention to $S_{\mathrm{tnl}}$ is taken so that $S_I$ becomes positive. Thus, the damping of the norm due to tunneling is estimated as

$$exp\left(-2\frac{\left|\mathrm{Im}{S}_{\mathrm{tnl}}\right|}{\hslash }\right),$$

(88)

which gives a natural definition of the relevant tunneling probability.

The tunneling paths induce additional phases when they are incorporated into the semiclassical Feynman kernel such that [124]

$$K(q_t,t;q_s,s:\sigma )=\sum _{\alpha }{\left(\frac{1}{-2\pi i\hslash }\right)}^{\frac{N}{2}}\frac{1}{\sqrt{\left|\frac{\partial q_t}{\partial {\overline{p}}_s}\right|}}exp\left(\frac{i}{\hslash }{S}_{\mathrm{tnl}}+i\varphi \right),$$

(89)

where N is the number of dimensions and the additional phase $\varphi$ is given by

$$\varphi =-\frac{\pi }{2}{N}_f+\frac{\pi }{4}{N}_p$$

(90)

at an entrance of a tunneling path and

$$\varphi =-\frac{\pi }{2}{N}_f-\frac{\pi }{4}{N}_p$$

(91)

at an exit of a tunneling path. Here $N_p$ is the number of the negative parities while $N_f$ is the sum of the multiplicities of the zeros of $\partial q_f/\partial {\overline{p}}_i=0$ along the trajectory, which gives the Maslov index.

Connection Problem

Tunneling paths can arise only upon transferring from $\left\{{\sigma }_k=1\right\}$ to another, typically only one of ${\sigma }_k$ is negative. This is because as the number of negative parities increases the damping probability in Equation (88) also become large. Those paths on the different sheets should be connected smoothly. To ensure smoothness, we change the sign of one of the parities at a caustic in the direction normal to the caustic line. The caustics, at which the density of the paths becomes very high (actually infinite in the primitive semiclassical approximation), are defined with the so-called Jacobi field [94] in such a way that

$$det|\frac{\partial \overrightarrow{q_f}}{\partial \overrightarrow{p_i}}|=0,$$

(92)

where $\overrightarrow{p_i}$ is a real-valued (quasi-)momentum vector for a given parity set {${\sigma }_k=\pm 1$} at the initial point on each sheet and $\overrightarrow{q_f}$ is the final points of a trajectory on the same sheet. The tunneling path can get back to the Newtonian sheet of $\left\{{\sigma }_k=1\right\}$ in a similar manner [125].

The stationary phase condition in the semiclassical approximation to the Feynman kernel requires that the quantum phase $S_{HJ}$ in Equation (89) should be connected smoothly at the connection point. This requirement is expressed as

$$\frac{\partial S_{HJ}}{\partial q}\left(\mathrm{in}\mathrm{classical}\mathrm{region}\right)=\frac{\partial S_{\mathrm{tnl}}}{\partial q}\left(\mathrm{in}\mathrm{tunnel}\mathrm{region}\right),$$

(93)

which simply implies

$$\overline{p_k}=-\overline{p_k}$$

(94)

at the connection point. Obviously,

$$\overline{p_k}=0$$

(95)

is the unique solution [126]. Therefore, to connect the wavefunctions smoothly at the turning points in the relevant direction k, we may change the parity there (or parities depending on the degeneracy of turning points) to the direction of zero momentum.

Tunneling phenomena are quite important in many chemical reaction dynamics, and some applications of the present treatment are found in Refs. [124,126,127,128].

5.5.2. Statistical Redistribution of Tunneling Paths

In view of the nature of this review article, we briefly show an example of chaos that is induced by tunneling dynamics, just for curiosity. The Hamiltonian is

$$H=\frac{1}{2m_x}{p}_x^2+\frac{1}{2m_y}{p}_y^2+\frac{1}{2}(x^2+y^2)+x^{2}y-\frac{1}{3}{y}^3,$$

(96)

where the masses are chosen as $m_x=1.0087$ and $m_y=1.0$ so as to break the symmetry of $C_3$. The vibrational modes and the corresponding tori in the Poincaré surface of section are depicted in Figure 10. The bottom panel displays a long time Poincaré section arising from a single trajectory, which repeats dynamical tunneling from time to time among tori, thus giving rise to chaos effectively. Here, however, we have intensionally disregarded the damping of the amplitude induced by tunneling. Careful inspection will find the faint trace of a couple of tori [125]. See ref. [125] for details about the present tunneling dynamics.

6. Chaos Arising from Repeated Bifurcation and Merge of the Quantum Wavepackets on Nonadiabatically Coupled Potential Basins

One of the characteristic features of quantum chaos is reflected on the energy level statistics such as the patterns in the nearest neighbor level spacing distribution (NNLSD) [114,129]. The extensive level repulsions among the states are widely recognized to be the origin of the Wigner distribution in NNLSD. It is therefore natural to pay attention to the nonadiabatic interactions as one of the key factors to modulate the structure of the energy levels. As such, the theory of level dynamics [9,130,131,132] sets a focus on nonadiabatic interactions as a factor that can modulate the energy levels in an anomalous fashion. Shudo [133,134] explicitly clarified the effects of an avoided crossing on the level distributions, the avoided crossing of which is induced by a variation of a part (s) of interaction potential but not by those interactions presented in Equation (12). The nonadiabatic interactions of Equation (12) are basically about electronic-state mixing due to the nuclear kinematic coupling on the occasion of the variation of nuclear configuration (molecular geometry) and are one of the most critical interactions in chemical dynamics [14,16,17,18,19,,20]. Therefore, the beautiful methods in the level dynamics does not seem to apply to the dynamics of molecules.

Besides, we are quite often concerned with the wavepacket dynamics of electrons rather than the electronic energy levels at a given nuclear configuration. In other words, the essential characteristics of molecular nonadiabatic dynamics is seen when the relevant quantum wavepackets bifurcate and coherently merge (remix) many times at the molecular geometry of significant nonadiabatic interactions. It is quite natural therefore to suspect that the wavepacket bifurcation and merging should be the dynamical origin of purely quantum chaos in place of “stretching and folding” in the baker’s transformation for classical chaos. We hence survey the nuclear wavepacket dynamics of bifurcation and remixing in this section first and then proceed to the nonadiabatic electron wavepacket dynamics in the next section.

6.1. Experimental Observation of Bifurcation and Merge of the Quantum Wavepackets

We here show that “bifurcation and merging” of quantum wavepackets can be observed with the modern real-time experiments. Among the cutting-age experimental methodologies to track the femtosecond-scale nuclear motions in chemical dynamics, the time- and angle-resolved photoelectron spectroscopy (TRPES) has been very well established since 1999 [135,136,137], when the group of this author happened to start the first complete ab initio full quantum mechanical simulation to track nuclear wavepackets in TRPES [138,139,140]. Among others, we showed for the first time that the instance of bifurcation and merging of nuclear wavepackets in passing across the nonadiabatic regions can be directly observed [141,142]. That is, such bifurcation and merging are of physical reality. Now other novel experimental methodologies are available.

For readers who might be interested in the recent experimental progress since then, which are related to observation nonadiabatic dynamics, we refer to the following works (we apologize for not a complete list). T. Suzuki is the first who used photoelectron imaging technique to observe the dynamics passing through an intersystem crossing of pyradine [143,144]. Wörner and his group observed the passage of NO${}_2$ molecule across the conical intersection by means of the high harmonic spectroscopy [145] and time-resolve photoelectron spectroscopy [146]. Further novel methods have been introduced, such as attosecond stimulated X-ray Raman spectroscopy [147], time-resolved fluorescent spectroscopy [148], ultrafast electron diffraction technique [149], and so on. We also found that the wavepacket bifurcation can be directly observed not only by pump-probe technique by means of pulse lasers but with induced photoemission spectroscopy from molecules placed under a CW laser [150].

The basic experimental setting for the time-resolved photoelectron spectroscopy is as follows: (1) First excite a ground state wavepacket to an electronically excited state by a pump pulse-laser. This laser is of a sufficiently short width so that the nuclear wavepacket thus excited begins to run spontaneously, which eventually passes across the nonadiabatic region(s). Then, (2) shine another laser (probe laser) to ionize the excited state. (3) Measure the energy and angular distributions of photoelectrons as a function of the time delay (the time lag between the probe and pump), and thereby retrieve the information of the dynamics of the nuclear (chemical) dynamics. This TRPES technique is quite successful in that there is no closed channel to be probed.

The theoretical framework of TRPES is outlined below [138,139,140]. Briefly, the total wavefunction is expanded in terms of the three relevant electronic states: ${\Phi }_1$ (with a potential function $V_1$), ${\Phi }_2$ (with $V_2$), and ${\Phi }_{\mathbf{k}}^{(-)}$ (ion state), as

$$\Psi (\mathbf{r},R,t)={\chi }_1(R,t){\Phi }_1(\mathbf{r};R)+{\chi }_2(R,t){\Phi }_2(\mathbf{r};R)+\int d\mathbf{k}{\chi }_{\mathbf{k}}(R,t){\Phi }_{\mathbf{k}}^{(-)}(\mathbf{r};R),$$

(97)

where ${\Phi }_{\mathbf{k}}^{(-)}(\mathbf{r};R)$ is a scattering state having the in-going boundary condition [49] with $\mathbf{k}$ labeling the wave vector of the photoelectron. The time-dependent coupled Schrödinger equations for the vibrational wavepackets ${\chi }_1$, ${\chi }_2$, and ${\chi }_{\mathbf{k}}$ are to be solved in a coupled Schrödinger equations (note that ${\chi }_{\mathbf{k}}$ is continuous with respect to $\mathbf{k}$) under all the matrix elements being evaluated in full quantum mechanics [138,139,140]. For instance, we need the matrix elements with respect to the electronic continuum representing the electronic ionization process. This is not an easy task. Once they are all solved through, the photoelectron spectra are calculated from the vibrational wavepackets after the probe interaction as

$$P\left({\epsilon }_k\right)=k\int dRd\widehat{k}{\left|{\chi }_{\mathbf{k}}\left(R\right)\right|}^2.$$

(98)

Figure 11 shows the photoelectron spectra of NaI molecule exhibiting the wavepacket bifurcation and merging, taken for the pump-probe delays of 675, 700, 725, and 750 fs. Snapshots of the absolute square of the wavepackets on the diabatic potentials are shown in the upper column. The red and blue ones represent ${\left|{\chi }_1^{\left(dia\right)}(R,t)\right|}^2$ on $V_1^{\left(dia\right)}\left(R\right)$ and ${\left|{\chi }_2^{\left(dia\right)}(R,t)\right|}^2$ on $V_2^{\left(dia\right)}\left(R\right)$, respectively. These wavepackets are passing across the crossing point (nonadiabatic region) from the right to the left. The photoelectron signals, the kinetic energy distribution of photoelectrons, at each time are correspondingly draw in the lower column. It is clearly seen that the peaks of photoelectron signals are modulated in real time from one to two and three peaks, highlighting the wavepacket bifurcation to those components. These dynamics can be comprehended primarily with the Condon principle: The wavepackets running on the lower potential energy make the higher photoelectron kinetic energy signals. The wavepacket bifurcation and merging and can thus be physically observed in real time.

6.2. Chaotic Eigenfunctions in Nonadiabatically Coupled Potential Functions

We next survey how the infinite repetition of bifurcation and merge of the quantum wavepackets results in the specific signature in eigenfunctions. Fujisaki systematically investigated this aspect by placing a focus on the transition from regular to chaotic states of the wavefunctions as functions of energy and the magnitude of nonadiabatic coupling elements [151,152]. The system used was the so-called Heller model, which was defined in their pioneering work [153]. It is defined as (see also Figure 12)

$$\begin{array}{cc}\hfill \mathcal{H}& =\left(\begin{array}{cc}T+V_A& J\\ J& T+V_B\end{array}\right)\hfill \\ \hfill T& =\frac{1}{2}\left(p_x^2+p_y^2\right),V_A=\frac{1}{2}\left({\omega }_x^2{x}^2+{\omega }_y^2{y}^2\right)+{ϵ}_A,V_B=\frac{1}{2}\left({\omega }_x^2{\xi }^2+{\omega }_y^2{\eta }^2\right)+{ϵ}_B\hfill \end{array}$$

(99)

with $\xi =(x+2asin\theta )cos2\theta +(y-2acos\theta )sin2\theta$ and $\eta =-(x+2asin\theta )sin2\theta +(y-2acos\theta )cos2\theta ,$ where T is a nuclear kinetic energy (nuclear masses are assumed to be unity), $V_i$ $(i=A,B)$ two-dimensional diabatic harmonic potentials of each electronic state i, J the nonadiabatic coupling constant between the electronic states induces an interaction between the two otherwise uncoupled harmonic potentials. ${ϵ}_i$$(i=A,B)$ indicates the bottoms of each harmonic potential. The geometrical meanings of a and $\theta$ are shown in Figure 12. The angle $2\theta$ referred to as the Duschinsky angle cause a coupling between x and y. Unlike the original Heller model, we chose the frame for $V_A$ to be ‘diagonalized’ (no crossing terms like $xy$ in $V_A$). The location of the crossing manifold (usually called crossing seam) should be presented in the original paper [152].

Since the Heller model is simply an uncoupled pair of harmonic potentials at $J=0$, it can have merely regular states depending on the values on J and $\theta$. Yet, the intense mixing induced by these coupling elements lead the wavefunctions and energy spectra to exhibit chaotic patterns. As an example, Figure 13 demonstrates the spatial distribution of some (selected) very chaotic wavefunctions. The chaoticity was quantified in terms of a transition from the Poisson to Wigner distributions in the nearest neighbor level spacing distribution (NNLSD), ${\Delta }_3$ statistics, and so on [151]. Here in the graph of Figure 13, only the so-called amplitude distributions of the eigenfunctions, proposed by Berry [154], on one of the surfaces (A) are shown. The definition of the amplitude distributions is

$$P\left(\varphi \right)=\frac{1}{\sqrt{2\pi {\sigma }^2}}exp\left(-\frac{{\varphi }^2}{2{\sigma }^2}\right),$$

(100)

where $P\left(\varphi \right)$ is the frequency distribution of the amplitude of the eigenfunction that is randomly sampled at points that are energetically accessible in configuration space on each potential. About 6000 points are randomly sampled in Figure 13, and it turns out that $P\left(\varphi \right)$ is found to be fitted with a large $\sigma$, reflecting the spatially wide distribution due to chaos. Although not shown here, we found that the regular wavefunctions give much smaller $\sigma$ (far narrower distributions) [151,152]. Thus, we have seen a regular to chaotic transition of eigenfunctions (and eigenvalues as well, not shown here) in nonadiabatically coupled quantum system.

6.3. Need for Measures of Chaoticity in Quantum Wavepackets

One of the most interesting issues in quantum chaos is how to measure the extent of chaoticity. This would depend on the definition of quantum chaos [2,3,4,5,6,7,8,9,10,11,12]. Characterization of chaoticity embedded in discrete quantum energy levels has been extensively studied (see Refs. [5,9,11] for instance). For example, the measures such as NNLSD, ${\Delta }_3$ statistics, and so on, work well to capture chaos. On the other hand, we need to characterize chaoticity in the dynamics of wavefunctions of molecules. Berry’s amplitude distributions of the eigenfunctions $P\left(\varphi \right)$ in Equation (100) can be used as one of such measures by applying at each time of the time-propagation of a packet. However, this measure seems to be based on a rather intuitive expectation that the chaotic wavefunction should be always widely distributed in space. Another concern is whether it can distinguish well between regular systems and those of quantum localization and the quantum scar.

The studies on classical chaos have already established such measures like the Lyapunov exponent [155], the real-complex valuedness of the eigenvalues of the stability matrix, KS entropy, and so on [1]. In practical studies on chemical dynamics we often need a crisp measure to detect and/or characterize chaoticity in the propagation of wavefunctions. However, it is known that quantum coherence working in quantum wavepackets particularly make a significant difference from the classical one. Suppose we have two wavepackets ${\varphi }_{\mu }\left(t\right)$ and ${\varphi }_{{\mu }^{\prime }}\left(t\right)$, which are slightly different from each other in their initial conditions $\mu$ and ${\mu }^{\prime }$, the initial condition of two wavepackets for instance. In marked contrast to the exponential deviation of two classical trajectories of slightly different initial conditions, the following overlap between the two packet

$$O_{\mu {\mu }^{\prime }}\left(t\right)={\left|〈\left.{\varphi }_{\mu }\left(t\right)\right|{\varphi }_{{\mu }^{\prime }}\left(t\right)〉\right|}^2$$

(101)

is kept constant as long as they are subject to the same Hamiltonian. For the same reason, quantities like the Shannon entropy

$$S\left(t\right)=-\int dq{\left|{\varphi }_{\mu }(q,t)\right|}^{2}log{\left|{\varphi }_{\mu }(q,t)\right|}^2$$

(102)

is conserved as well. Then, a primitive idea to break the coherence is as follows [156]. Prepare a projection operator

$$P_i={\int }_{A_i}dr\left|r\right.〉\left.〈r\right|$$

(103)

where $A_i$ is a subset (local area) of the entire space and satisfies $A_i\cap A_j=\varnothing$ for $i\ne j$. In place of the overlap integral $O_{\mu {\mu }^{\prime }}\left(t\right)$

$$O_{\mu {\mu }^{\prime }}\left(t\right)={\left|〈{\varphi }_{\mu }\left(t\right)\left|\sum _i{P}_i\right|{\varphi }_{\mu }\left(t\right)〉\right|}^2$$

(104)

we consider the following modulation

$$D_{\mu {\mu }^{\prime }}\left(t\right)=\sum _i{\left|〈{\varphi }_{\mu }\left(t\right)\left|P_i\right|{\varphi }_{{\mu }^{\prime }}\left(t\right)〉\right|}^2,$$

(105)

which obviously destroys the coherence among the different regions. In chaos ${\varphi }_{\mu }\left(t\right)$ and ${\varphi }_{{\mu }^{\prime }}\left(t\right)$ deviate from each other in configuration space, and $D_{\mu {\mu }^{\prime }}\left(t\right)$ is anticipated to exhibit a quick damp from a slight difference in the initial condition [156]. A good compromise to balance the number of $A_i$ and size is necessary to sensitive yet economical detection of chaos. It has been shown numerically that the method works very well even for the nonadiabatic dynamics as in the above Heller model system [156]. However, this decoherence method costs much labor in many dimensional systems. We hence need general and tractable alternatives. We will be back to this aspect in the next section.

7. Chaos in Nonadiabatic Electron Dynamics of Molecules

In this section, we review chaotic dynamics in electronic state manifolds, those states involved in which are mixed frequently and extensively by the nuclear kinematic couplings, or as vaguely termed, nonadiabatic interactions [21,23,24,25,26,27]. (We do not review the traditional and canonical theories of nonadiabatic transitions for nuclear wavefunctions, essentially originating from the spirit of the Landau-Zener-Stuekelberg theory [14,16,17,19,20]).

7.1. The Hamiltonian Studied for Electron Dynamics

7.1.1. Electron Dynamics

The coupled Equations (11) based on the Born-Huang expansion, Equation (9), set the theoretical foundation of molecular science. Yet, various drawbacks arise coming from the time-independent representation of the electronic wavefunctions $\left\{{\Phi }_I(\mathbf{r};\mathbf{R})\right\}$. First of all, it is technically hard to prepare a set of potential functions $\left\{V_I\left(\mathbf{R}\right)\right\}$, on which the nuclear wavepackets $\left\{{\chi }_I(\mathbf{R},t)\right\}$ propagate in time. In particular, the most serious situation appears when the electronic states are densely quasi-degenerate and are strongly coupled to one another through the nonadiabatic interactions. In such a case a huge time-dependent fluctuation over the electronic states can appear. To cope with those difficulties, we resume the dynamics with time-dependent electronic wavefunctions from the beginning.

We hereby restart our basic framework from the following representation of the total Hamiltonian operator [21,157]

$$\begin{array}{c}H(\mathbf{R},\mathrm{elec})\hfill \\ \equiv \frac{1}{2}\sum _k{\left({\widehat{P}}_k-i\hslash \sum _{IJ}\left|{\Phi }_I\right.〉X_{IJ}^k\left.〈{\Phi }_J\right|\right)}^2+\sum _{IJ}\left|{\Phi }_I\right.〉H_{IJ}^{\left(\mathrm{el}\right)}\left.〈{\Phi }_J\right|\hfill \end{array}$$

(106)

where ${\widehat{P}}_k$ is the nuclear momentum operator. $\mathbf{R}$ indicates the nuclear coordinates, and we represent the electronic states in the Hilbert space, not in the configuration space $\mathbf{r}$. Note that ${\widehat{P}}_k$ are to be operated only on the nuclear wavefunctions, not on the electronic wavefunctions ${\Phi }_I(\mathbf{r};\mathbf{R})$, since such kinematic interactions have already been taken into account in $X_{IJ}^k$ of Equation (106). The sum over the electronic states suffixed with I and J may include continuum states.

No approximation has been made so far, and, we here introduce an approximation that the nuclear momentum operators ${\widehat{P}}_k$ are replaced with their classical counterparts $P_k$ as

$$\begin{array}{c}\tilde{H}(\mathbf{R},\mathbf{P},\mathrm{elec})\hfill \\ \equiv \frac{1}{2}\sum _k{\left(P_k-i\hslash \sum _{IJ}\left|{\Phi }_I\right.〉X_{IJ}^k\left.〈{\Phi }_J\right|\right)}^2+\sum _{IJ}\left|{\Phi }_I\right.〉H_{IJ}^{\left(\mathrm{el}\right)}\left.〈{\Phi }_J\right|\hfill \end{array}$$

(107)

Then the total wavefunction subject to the Hamiltonian of Equation (106) is to be reduced to an electron wavepacket state running on paths $\mathbf{R}\left(t\right)$ (but not necessarily Newtonian trajectories). More explicitly, we expand it as

$$\left|{\Phi }_{\mathbf{elec}}(t;\mathbf{R}\left(t\right))\right.〉=\sum _I{C}_I\left(t\right)\left|{\Phi }_I(\mathbf{r};\mathbf{R})\right.〉{|}_{\mathbf{R}=\mathbf{R}\left(t\right)},$$

(108)

along each $\mathbf{R}\left(t\right)$, where $\left\{{\Phi }_I(\mathbf{r};\mathbf{R})\right\}$ is a basis functions taken at each nuclear coordinate $\mathbf{R}$, and ${\Phi }_I(\mathbf{r};\mathbf{R})$ do not have to be the eigenfunctions of the electronic Hamiltonian. Equations (107) and (108) are widely referred to as mixed quantum-classical representation. The time-dependent coefficients $C_I\left(t\right)$ are to be evaluated with the coupled equations of motion for electron wavepackets in such a way that

$$i\hslash \frac{d{C}_I}{dt}=\sum _J\left[H_{IJ}^{\left(\mathrm{el}\right)}-i\hslash \sum _k{\dot{R}}_k{X}_{IJ}^k-\frac{{\hslash }^2}{4}\sum _k(Y_{IJ}^k+Y_{JI}^{k*})\right]C_J,$$

(109)

where $H_{IJ}^{\left(\mathrm{el}\right)},$$X_{IJ}^k,$ and $Y_{IJ}^k$ have been defined in Equations (12) and (13), respectively. Again, the electronic basis $\left\{{\Phi }_I(\mathbf{r};\mathbf{R})\right\}$ are not demanded to be adiabatic ones. The bra-ket notation used here indicates integration over the electronic coordinates $\mathbf{r}$ only. Furthermore, the second order derivative terms $Y_{IJ}^k$ in Equation (109) are quite often neglected in practice because they are always accompanied by the small quantity ${\hslash }^2$ [24,25], although it is not always negligible in general.

7.1.2. Nuclear Dynamics: Path Branching

The nuclear paths $\mathbf{R}\left(t\right)$ are driven by the force field provided by the electronic states, more explicitly in terms of the force matrix $F_{IJ}^k$ (see ref. [20]) defined as

$$F_{IJ}^k=-\frac{\partial H_{IJ}^{\left(\mathrm{el}\right)}}{\partial R_k}-\sum _K\left(X_{IK}^k{H}_{KJ}^{\left(\mathrm{el}\right)}-H_{IK}^{\left(\mathrm{el}\right)}{X}_{KJ}^k\right)+i\hslash \sum _l{\dot{R}}_l\left[\frac{\partial X_{IJ}^l}{\partial R_k}-\frac{\partial X_{IJ}^k}{\partial R_l}\right],$$

(110)

which arises from the variational principle or the Hamilton canonical equations of motion for $(\mathbf{R},\mathbf{P})$ [20].

The non-zero off-diagonal elements of $F_{IJ}^k$ can induce branching of the nuclear paths $\mathbf{R}\left(t\right)$ [21]. If, on the other hand, the force matrix is diagonal, that is, $F_{IJ}^k=0$ for $I\ne J$, each of $F_{II}^k$ gives a Newtonian force on each potential function I, which gives the theoretical foundation for the so-called ab initio dynamics, molecular dynamics, the first principle dynamics, and so on.

7.1.3. Nuclear Dynamics: Mean-Field Approximation (Semiclassical Ehrenfest Theory)

To save the work needed for the path branching calculations [22], we often adopt a mean-field approximation, the so-called semiclassical Ehrenfest theory (SET), in which the force matrix of Equation (110) is to be averaged over the electron wavepacket as

$$\begin{array}{c}\tilde{H}(\mathbf{R},\mathbf{P},\mathrm{elec})\hfill \\ \equiv \frac{1}{2}\sum _k{\left(P_k-i\hslash \sum _{IJ}\left|{\Phi }_I\right.〉X_{IJ}^k\left.〈{\Phi }_J\right|\right)}^2+\sum _{IJ}\left|{\Phi }_I\right.〉H_{IJ}^{\left(\mathrm{el}\right)}\left.〈{\Phi }_J\right|\hfill \end{array}$$

(111)

which gives rise to a single force for a single path. If the basis set $\left\{{\Phi }_I(\mathbf{r};\mathbf{R})\right\}$ happens to be complete, Equation (111) is reduced to the form of Hellmann–Feynman force

$$\langle F_k\rangle =-〈{\Psi }_{\mathrm{elec}}\left|\frac{\partial {\widehat{H}}^{\left(\mathrm{el}\right)}}{\partial R_k}\right|{\Psi }_{\mathrm{elec}}〉.$$

(112)

The use of the mean-field approximation is justified for short-time propagation, typically shorter than 100 fs. Indeed, it has been numerically shown [21] that the electronic-state mixing is very well (accurately) reproduced by SET for a short time, provided that the nonadiabatic coupling elements $X_{IK}^k$ and electronic matrix elements $H_{IK}^{\left(\mathrm{el}\right)}$ are both correctly evaluated. As the nuclear positions proceed to some extent, the Ehrenfest paths should be modified accordingly [22].

7.2. Longuet–Higgins (Berry) Phase and Lorentz-like Force

Before proceeding to the electronic-state chaos, we make a detour to a specific force working on the nuclei, which can cause breaking the mirror symmetry of molecular shape, namely molecular chirality.

Let us review the molecular Hamiltonian of Equation (107) now in classical electromagnetic field

$$\begin{array}{c}\tilde{\mathcal{H}}(\mathbf{R},\mathbf{P},\mathrm{elec},\mathbf{A})\equiv \hfill \\ \frac{1}{2}\sum _k{\left(P_k-\frac{Z_{k}e}{c}{A}_k-i\hslash \sum _{IJ}\left|{\Phi }_I\right.〉X_{IJ}^k\left.〈{\Phi }_J\right|\right)}^2+\sum _{IJ}\left|{\Phi }_I\right.〉{\mathcal{H}}_{IJ}^{\left(el\right)}\left({\mathbf{A}}^{ele}\right)\left.〈{\Phi }_J\right|\hfill \end{array}$$

(113)

where $A_k\left(\mathbf{R}\right)$ are the (classical) electromagnetic vector potentials on nuclei. The electronic Hamiltonian also includes the vector potential ${\mathbf{A}}^{ele}$. In this expression we notice that the nonadiabatic couplings mechanically stand parallel with the vector field $\mathbf{A}$, which suggests the qualitative roles of the nonadiabatic couplings working on the nuclear momentum may be akin to those of the vector potentials.

Indeed it is well-known that there is a mathematical similarity between the Aharonov–Bohm phase [158] and the Longuet–Higgins phase [159], the latter being well known in chemical dynamics relevant to the so-called conical intersection, which is a special case of nonadiabatic interaction. (see ref. [17] for the interesting properties of conical intersection.) These two phases have been unified as the so-called geometrical phase, often called Berry’s phase [160,161]. In light of the above parallelism in Equation (113), we suspect another similarity to exist between them: As the electromagnetic Lorentz force works on a charged particle running in a magnetic field, a novel force, which we call the Lorentz-like nonadiabatic force, appears to work on nuclei running in a field created by the nonadiabatic interaction field [21,162].

To be more specific, let us imagine a molecular system as illustrated in Figure 14.

At the center we have an atom C, say carbon atom, and around it lie three atomic or molecular elements A${}_1$, A${}_2$, A${}_3$. All the elements A${}_1$, A${}_2$, A${}_3$ in the molecule CA${}_1$A${}_2$A${}_3$ are mutually different and sit stably on a common plane. Suppose an external atom A is coming in on the same plane of an $x-y$ coordinates. (A${}_1$, A${}_2$, A${}_3$, and A should be all different species from each other.) The z-direction is set to be perpendicular to this plane and the plus or minus of this direction determines the right or left with respect to the plane. Let $\left(x_a,y_a,z_a\right)$ be a point of the incident atom A. Then the force matrix in $z_a$-direction (perpendicular to the plane) [162] is given such that

$$\begin{array}{cc}\hfill F_{IJ}^{z_a}=& -{\left(\frac{\partial H^{\left(el\right)}}{\partial R_{z_a}}\right)}_{IJ}-i\hslash \left({\dot{R}}_{x_a}\left[\frac{\partial X_{IJ}^{x_a}}{\partial R_{z_a}}-\frac{\partial X_{IJ}^{z_a}}{\partial R_{x_a}}\right]+{\dot{R}}_{y_a}\left[\frac{\partial X_{IJ}^{y_a}}{\partial R_{z_a}}-\frac{\partial X_{IJ}^{z_a}}{\partial R_{y_a}}\right]\right)\hfill \\ & -i\hslash \sum _{b\ne a}^{{\mathrm{A}}_1,{\mathrm{A}}_2,{\mathrm{A}}_3}\sum _{j_b}^{x,y,z}{\dot{R}}_{j_b}\left[\frac{\partial X_{IJ}^{j_b}}{\partial R_{z_a}}-\frac{\partial X_{IJ}^{z_a}}{\partial R_{j_b}}\right]\hfill \end{array}$$

(114)

where the suffix b specifies A${}_1$, A${}_2$, and A${}_3$, and $j_b$ indicates the $\left(x,y,z\right)$-coordinate of each element. This is a force while the state undergoes nonadiabatic transition from an Ith to the Jth one. Just for simplicity, we assume that all the atoms other than the incident atom $\mathrm{A}$ stay still, thus neglecting the possible vibrational motion. For the sake of clear presentation, we simply set to

$${\dot{R}}_{j_b}=0(j=x,y,z\mathrm{and}b\ne a).$$

(115)

Hence only the atom $\mathrm{A}$ is supposed to be moving on the plane with a nonzero value of $\left({\dot{R}}_{x_a},{\dot{R}}_{y_a}\right)$ and ${\dot{R}}_{z_a}=0$. $F_{IJ}^{z_a}$ is then reduced to

$$F_{IJ}^{z_a}=-{\left({\displaystyle \frac{\partial H^{\left(el\right)}}{\partial R_{z_a}}}\right)}_{IJ}-i\hslash \left({\dot{R}}_{x_a}\left[{\displaystyle \frac{\partial X_{IJ}^{x_a}}{\partial R_{z_a}}}-{\displaystyle \frac{\partial X_{IJ}^{z_a}}{\partial R_{x_a}}}\right]+{\dot{R}}_{y_a}\left[{\displaystyle \frac{\partial X_{IJ}^{y_a}}{\partial R_{z_a}}}-{\displaystyle \frac{\partial X_{IJ}^{z_a}}{\partial R_{y_a}}}\right]\right).$$

(116)

We here define the Lorentz-like nonadiabatic force with the second term of the right hand side of this equation as [162]

$$f_{IJ}^{z_a}\equiv -i\hslash \left({\dot{R}}_{x_a}\left[{\displaystyle \frac{\partial X_{IJ}^{x_a}}{\partial R_{z_a}}}-{\displaystyle \frac{\partial X_{IJ}^{z_a}}{\partial R_{x_a}}}\right]+{\dot{R}}_{y_a}\left[{\displaystyle \frac{\partial X_{IJ}^{y_a}}{\partial R_{z_a}}}-{\displaystyle \frac{\partial X_{IJ}^{z_a}}{\partial R_{y_a}}}\right]\right)$$

(117)

and simply write as

$$F_{IJ}^{z_a}=-{\left({\partial }_{z_a}{H}^{\left(el\right)}\right)}_{IJ}+f_{IJ}^{z_a},$$

(118)

where $-{\left({\partial }_{z_a}{H}^{\left(el\right)}\right)}_{IJ}$ is the short-hand notation of the off-diagonal Hellmann-Feynman term. (Note that the diagonal elements of the nonadiabatic interaction $X_{II}$ are zero if the function used as the basis functions are real. However, this is not the case for complex-valued wavefunctions.)

The expression in Equation (117) can be reexpressed in a more familiar manner: Let us consider a vector of the force matrix elements as

$${\overrightarrow{f}}_{IJ}^a=\left(f_{IJ}^{x_a},f_{IJ}^{y_a},f_{IJ}^{z_a}\right).$$

(119)

In terms of the vector ${\overrightarrow{X}}_{IJ}^{{}_a}\left(\mathbf{R}\right)$ of the nonadiabatic interaction fields

$${\overrightarrow{X}}_{IJ}^{{}_a}\left(\mathbf{R}\right)=\left(X_{IJ}^{x_a}\left(\mathbf{R}\right),X_{IJ}^{y_a}\left(\mathbf{R}\right),X_{IJ}^{z_a}\left(\mathbf{R}\right)\right)$$

(120)

we define a vector field

$${\overrightarrow{Z}}_{IJ}^{{}_a}\left(\mathbf{R}\right)=-i\hslash {\overrightarrow{\nabla }}^a\times {\overrightarrow{X}}_{IJ}^{{}_a}\left(\mathbf{R}\right),$$

(121)

in which the dependence on the nuclear coordinates $\mathbf{R}$ is explicitly denoted. Then ${\overrightarrow{f}}_{IJ}^a$ is expressed as

$${\overrightarrow{f}}_{IJ}^a\left(\mathbf{R}\right)={\overrightarrow{v}}^a\times {\overrightarrow{Z}}_{IJ}^{{}_a}\left(\mathbf{R}\right),$$

(122)

where ${\overrightarrow{v}}^a$ is the velocity vector of the atom $\mathrm{A}$, and ${\overrightarrow{\nabla }}^a\times$ indicates the curl on the coordinates for the atom $\mathrm{A}$ only [162].

It is immediately noticed that this expression is essentially of the same mathematical form as that behind the Fleming left-handed rule in classical electromagnetic theory. The Lorentz force distinguishes the front and back sides of a plane made up by an electron current and a magnetic field. See Figure 15. Thus, a simple relationship is found between the electromagnetic Lorentz force and the Lorentz-like nonadiabatic force in such a way that

$${\overrightarrow{A}}^a\leftrightarrow {\overrightarrow{X}}_{IJ}^{{}_a},{\overrightarrow{B}}^a\leftrightarrow {\overrightarrow{Z}}_{IJ}^{{}_a},{\dot{\mathbf{R}}}_a\times q_{a}e{\mathbf{B}}_a\leftrightarrow {\overrightarrow{v}}^a\times {\overrightarrow{Z}}_{IJ}^{{}_a},$$

(123)

where $\mathbf{B}$ is a magnetic field.

We thus recognize that the nonadiabatic force can distinguish right and left with respect to the plane of the vector field of nonadiabatic interactions, which can serve as one of the fundamental mechanisms of the origin of molecular chirality. Yet, further study is needed to quantify the actual magnitude of anisotropic force.

7.3. Chaotic Electron Dynamics in Densely Quasi-Degenerate Electronic-State Manifold; B${}_{12}$ Cluster as an Example

As a candidate among molecules possibly having electronic-state chaos, we study some of the highly excited states of a cluster consisting of 12 boron atoms, B${}_{12}$. The nonadiabatic interactions working among the excited states of this cluster cause very frequent transitions among adiabatic states [163,164,165], resulting in frequent and continual state-mixing as represented as (recall Equation (108))

$$\sum _I^{\mathrm{adiabatic}}{C}_I\left(t\right){\Phi }_I^{\mathrm{ad}}(\mathbf{r};\mathbf{R}\left(t\right)).$$

(124)

In contrast to classical chaos, there is no crisp criterion with which to define or identify chaos in quantum dynamics, to the best of our knowledge. Below we only quantify the extent of chaoticity of the relevant electron dynamics represented as in Equation (124). A special attention is paid to the random mixing among the deterministic electronic states in the state space. We then discuss the possible manifestation of those strong chaoticity in chemical and physical properties of the clusters.

7.3.1. Nearest Neighbor Level Spacing Distribution (NNLSD)

At any selected time $t,$ or at any molecular shape $\mathbf{R}$, one can diagonalize the electronic Hamiltonian matrix $H_{IJ}^{\left(\mathrm{el}\right)}$ in Equation (109). The energy levels are therefore geometry dependent. We have surveyed the NNLSD for the B${}_{12}$ cluster (see ref. [164] and its supplemental material stored in J. Chem. Phys.). It turns out that whereas the NNLSD shows a clear Wigner distribution in case of low symmetry of molecular shape, level clustering occurs to some extent when the molecule is of higher symmetry, actually C${}_{3v}$ symmetry. This is natural because the systems of higher symmetry tend to bear accidental degeneracy more often. On the other hand, the possibility (measure) for this system to take C${}_{3v}$ symmetry is actually zero in the dynamical studies of variable $\mathbf{R}\left(t\right)$. Thus, we can judge that the present system has a chaotic structure in its energy level. However, our goal is to clarify and characterize the chaotic feature possibly hidden in the individual electronic wavepackets.

7.3.2. Diffusive Dynamics in the State Space: Fractional Brown Motion

It is experimentally hard to pick and prepare a pure adiabatic excited state out of the densely quasi-degenerate electronic state. Even if one could successfully pump the ground state to one of the highly adiabatic states, the resultant state would quickly begin to undergo continual nonadiabatic mixing with other states, each of which in turn further mixes with other states. These enduring state-mixings look like a diffusion as a whole in the electronic Hilbert space. Figure 16 shows the magnitude of ${\left|C_I\left(t\right)\right|}^2$ of Equation (124) in a color code in the ($I,t$) space. The brighter color indicates the larger ${\left|C_I\left(t\right)\right|}^2$. It visibly demonstrates that diffusive dynamics has been certainly realized in the electronic-state space.

To quantify the extent of the diffusive motion of the wavepackets in the electronic state-space explicitly, we survey the scaling property of the dynamics in terms of the fractional Brownian motion (fBm) of Mandelbrot [166]. Consider

$$\sum _{J=1}{(J-I\left(0\right))}^2{\left|C_J\left(t\right)\right|}^2=D_{I\left(0\right)}{t}^{\alpha },$$

(125)

where $I\left(0\right)$ indicates the state number of the initial state (that is $C_I\left(0\right)$), and $C_J\left(t\right)$ is the mixing coefficient at time t starting from this initial $I\left(0\right)$ state. Again $C_J\left(t\right)$ have been evaluated with the SET electron wavepackets. Recall that the exponent $\alpha =1.0$ implies the ordinary diffusion and $D_{I\left(0\right)}$ is just the ordinary diffusion constant. States of the initial adiabatic state $I\left(0\right)=3,$$40,$$100,$$300,$ and 500 have been numerically examined [164]. Furthermore, it turns out that the nonadiabatic electron dynamics of $I\left(0\right)\ge$ 40 are indeed the fractional Brownian motion (fBm) in the sense of Equation (125). The exponent of the present fBm is found to be about $\alpha =0.15$, which suggests a milder stochastic motion than ordinary diffusion. Interestingly, the values of the exponent is mostly independent of the initial states chosen, whereas the diffusion constants $D_{I\left(0\right)}$ significantly depend on the initial states, with the higher states having the larger diffusion constant.

Thus, Figure 16 convinces enough that the concept of adiabatic state in the Born-Oppenheimer approximation and the associated the energy level (potential energy hypersurface) all lose sense in the present situation. Besides, the dynamical aspect of chaos is far more directly relevant to the present huge electronic-state fluctuation than the static signatures of chaos like the level statistics.

7.3.3. Shannon Entropy

We next calculate the Shannon entropy in terms of ${\left|C_I\left(t\right)\right|}^2$ taken from Equation (124), which is the probability for the electronic wavepacket to occupy the Ith adiabatic state at time t. Hence,

$$S=-\sum _I{\left|C_I\left(t\right)\right|}^{2}log{\left|C_I\left(t\right)\right|}^2$$

(126)

quantifies the information about randomness in the electronic Hilbert space. Note this entropy is different from Equation (102). Incidentally, the present system is too large to apply the spirit of decoherence as in Equation (105). Figure 17 shows S for the three wavepackets, those starting from the adiabatic state #100, #200, and #300. A very quick rise from $S=0$ is observed for all the wavepackets up to about 5 fs, and much slower uprises follow since then. The absolute magnitude of the entropy is found to be larger for the higher state, which seems natural. In case of the adiabatic dynamics, where no electronic-state mixing is induced, the entropy is kept zero since ${\left|C_I\left(t\right)\right|}^2=1.0$. We thus see that the present electron dynamics keep distributing a large amount of the population over the extensive range in the Hilbert space.

7.3.4. Lyapunov Exponents for the Loss of Electronic State Memory

In the present nonadiabatic electron wavepacket dynamics, the main origin of chaoticity is not nonlinearity in the nuclear path dynamics but the nonadiabatic electronic state-mixing. To quantify the chaoticity due to the nonadiabatic interactions directly, we compare two dynamics that start from the same initial condition and one running on an adiabatic potential and the other to run on the nonadiabatic averaged potential. Suppose we have the Ith pure adiabatic state ${\Phi }_I^{ad}(\mathbf{r};\mathbf{R}\left(t\right))$ at the nuclear position $\mathbf{R}\left(t\right)$ at time t. Let a classical path run on this potential surface $V_I\left(\mathbf{R}\right)$ for a short time $\Delta t$ with an arbitrary momentum $\mathbf{P}\left(t\right)$, reaching a point at ${\mathbf{R}}^{ad}(t+\Delta t)$ with the electronic state ${\Phi }_I^{ad}(\mathbf{r};{\mathbf{R}}^{ad}(t+\Delta t))$. Likewise we propagate a nonadiabatic electronic state starting from the same set of $\left({\Phi }_I^{ad}(\mathbf{r};\mathbf{R}\left(t\right)),\mathbf{R}\left(t\right),\mathbf{P}\left(t\right)\right)$ but the nonadiabatic interactions are taken account this time. The nuclear path then arrives at ${\mathbf{R}}^{SET}(t+\Delta t)$, which is slightly different from ${\mathbf{R}}^{ad}(t+\Delta t),$ but a significant different electronic wavepacket ${\Phi }_{\left(I\right)}^{SET}(\mathbf{r};{\mathbf{R}}^{SET}(t+\Delta t))$ is to be attained. The process is summarized as

$$\left\{\begin{array}{c}{\Phi }_I^{ad}(\mathbf{r};\mathbf{R}\left(t\right))\to {\Phi }_I^{ad}(\mathbf{r};{\mathbf{R}}^{ad}(t+\Delta t))\mathrm{by}\mathrm{adiabatic}\mathrm{path}\\ {\Phi }_I^{ad}(\mathbf{r};\mathbf{R}\left(t\right))\to {\Phi }_{\left(I\right)}^{SET}(\mathbf{r};{\mathbf{R}}^{SET}(t+\Delta t))\mathrm{by}\mathrm{nonadiabatic}\mathrm{path}.\end{array}\right.$$

(127)

The deviation of ${\Phi }_{\left(I\right)}^{SET}(\mathbf{r};{\mathbf{R}}^{SET}(t+\Delta t))$ from ${\Phi }_I^{ad}(\mathbf{r};{\mathbf{R}}^{ad}(t+\Delta t))$ should reflect a memory loss of the initial state, and we quantify this the process by regarding it as an exponential decay of

$$\begin{array}{c}{\left|〈{\Phi }_I^{ad}\left(\mathbf{r};{\mathbf{R}}^{ad}(t+\Delta t)\right)\mid {\Phi }_{\left(I\right)}^{SET}\left(\mathbf{r};{\mathbf{R}}^{SET}(t+\Delta t)\right)〉\right|}^2\hfill \\ ={\left|〈{\Phi }_I^{ad}(\mathbf{r};\mathbf{R}\left(t\right))\mid {\Phi }_I^{ad}(\mathbf{r};\mathbf{R}\left(t\right))〉\right|}^{2}exp\left(-{\xi }_I\Delta t\right)\hfill \end{array}$$

(128)

where t is chosen arbitrary. Unfortunately the integral in the left hand side of Equation (128) is not easy to estimate. However, for a sufficiently short $\Delta t$, the following approximation holds very accurately

$${\mathbf{R}}^{SET}(t+\Delta t)\simeq {\mathbf{R}}^{ad}(t+\Delta t),$$

(129)

since their deviation is far slower and smaller. In fact, the deviation is numerically found to be about 0.005 Å in 1.0 fs (not shown graphically). Thus, we may define

$${\xi }_I(t;\Delta t)=-\frac{1}{\Delta t}log{\left|〈\left.{\Phi }_I^{ad}(\mathbf{r};{\mathbf{R}}^{SET}(t+\Delta t))\right|{\Phi }_{\left(I\right)}^{SET}(\mathbf{r};{\mathbf{R}}^{SET}(t+\Delta t))〉\right|}^2,$$

(130)

and with use of Equation (124) it is estimated by

$${\xi }_I=-\frac{1}{\Delta t}log{\left|C_I(t+\Delta t)\right|}^2.$$

(131)

Figure 18 represents the process of the above memory-loss process as observed in the selected states. Panel (a) shows that the nonadiabatic electronic state starting from the adiabatic #300 one quickly loses its memory about which adiabatic state it emerged from [167]. Actually it takes only 0.5 fs for the nonadiabatic state to lose most of the memory in ${\left|C_I\left(t\right)\right|}^2.$ The relevant exponents defined in Equation (128) is shown in panel (b). In this calculation, $\Delta t$ was set to 0.5 fs and the time t was scanned. Panel (b) displays two exponents ${\xi }_I$ of Equation (131) for the nonadiabatic states starting from the #100 and #300 adiabatic states, respectively. They fluctuate from time to time with the average value around 5.0 fs${}^{-1}$. This implies that the adiabatic path can immediately lose its memory as soon as the nonadiabatic interactions are switched on, no matter where the paths run. In other words, the adiabatic states are entirely embedded in a wide chaotic sea.

7.3.5. Turbulent Electron Flow in the Cluster

To observe one of the physical consequences from the present stochastic electron dynamics, we next survey the internal electron current driven by the nonadiabatic interactions. The probability current density of electrons $\overrightarrow{j}(\mathbf{r},t)$, which is also termed electron flux, quantifies and visualizes the flow pattern. $\overrightarrow{j}(\mathbf{r},t)$ is supposed to satisfy the continuity equation

$$\frac{\partial \rho (\mathbf{r},t)}{\partial t}+\nabla \overrightarrow{j}(\mathbf{r},t)=0,$$

(132)

and in quantum mechanics $\overrightarrow{j}(\mathbf{r},t)$ is defined as

$$\overrightarrow{j}(\mathbf{r},t)=\frac{\hslash }{2i{m}_e}[{\psi }^{*}\nabla \psi -\psi \nabla {\psi }^{*}],$$

(133)

where $m_e$ and $\psi$ are the mass and wave function of the involved particles [168]. For many-electron systems it is reduced to

$$\begin{array}{cc}\hfill \overrightarrow{j}(\mathbf{r},t)& ={\left.\frac{\hslash }{2i{m}_e}\left[{\nabla }_{\mathbf{r}}\rho \left({\mathbf{r}}^{\prime },\mathbf{r}\right)-{\nabla }_{{\mathbf{r}}^{\prime }}\rho \left({\mathbf{r}}^{\prime },\mathbf{r}\right)\right]\right|}_{{\mathbf{r}}^{\prime }\to \mathbf{r}}\hfill \\ & ={\left.\frac{\hslash }{2i{m}_e}\sum _i^{\mathrm{NO}}{n}_i\left[{\lambda }_i^{*}\left({\mathbf{r}}^{\prime }\right){\nabla }_{\mathbf{r}}{\lambda }_i\left(\mathbf{r}\right)-{\nabla }_{{\mathbf{r}}^{\prime }}{\lambda }_i^{*}\left({\mathbf{r}}^{\prime }\right){\lambda }_i\left(\mathbf{r}\right)\right]\right|}_{{\mathbf{r}}^{\prime }\to \mathbf{r}},\hfill \end{array}$$

(134)

where ${\nabla }_{\mathbf{r}}$ and ${\nabla }_{{\mathbf{r}}^{\prime }}$ are the nabla with respect to $\mathbf{r}$ and ${\mathbf{r}}^{\prime }$, respectively, and $\left\{{\lambda }_i\right\}$ are the natural orbitals that are eigenfunction of the density operator $\widehat{\rho }\left(t\right)$ ($\rho ({\mathbf{r}}^{\prime },\mathbf{r})=〈{\mathbf{r}}^{\prime }\left|\widehat{\rho }\right|\mathbf{r}〉$) with occupation numbers $n_i$. Note that only the complex-valued natural orbitals can make contributions to Equation (134). Stationary-state electronic wavefunctions like most of the eigenfunctions of the electronic Hamiltonian ${\widehat{H}}^{\left(el\right)}$ are real-valued (with exceptions with respect to angular momentum eigenfunctions) and thereby give zero flux only (see [169,170] for the general electronic and nuclear flux in general molecules, refs. [171,172,173,174] for the early applications of electron flux in chemical dynamics, ref. [175] for an extensive review, and ref. [176] for one of the latest).

Figure 19 shows the electron flux induced by the electron wavepacket dynamics starting from the 74th and 75th states, individually, with no nuclear kinetic energy but released spontaneously from the initial geometry. The snapshots have been taken at two very short times at 0.6 fs and 1.0 fs. These flux vectors are viewed from the top. The packet starting from the 74th state happens to exhibit the electron current in the counterclockwise direction at 0.6 fs, while the 75th counterpart shows the current mostly in the direction clockwise. Although the cluster seems to be in C${}_{3v}$ symmetry, it is relaxed to the lower one, and consequently, these circular patters are not conserved and shortly transformed to a more complicated flow. Indeed, only shortly after 1.0 fs both fluxes begin to fall into very complicated in spatial patterns. After all, those electron currents become turbulent or random flow (not shown graphically). The most exciting chemical study is then to identify the consequences of electronic turbulence in chemical reactivity and molecular properties.

Incidentally, one can define the flux of electronic energy [177], which is of course different from the above electron flux. It is quite interesting to track real-time energy flow in molecules. We have observed the turbulent flux of the electronic energy too in the above excited-state cluster dynamics [177].

We have thus observed the clear signatures of quantum chaos in molecular nonadiabatic electron dynamics.

7.4. The Long Life-Time of Dynamical Chemical Bonding: Hyper Resonance

The above diffusive and chaotic phenomena due to the very frequent nonadiabatic interactions prevent molecular dissociation and bring about a notion of quantum resonance we have not seen before. As resonance phenomena the Pauling resonance in the theory of chemical electronic structure [178] and the Feshbach resonance [49] in quantum scattering theory [50] are most relevant to the present study. The Pauling resonance in the so-called valence bond theory is regarded as a configuration mixing (interaction) among the so-called resonance structures, each of which corresponds intuitively to a valence bonding. The Feshbach resonance, or core-excited resonance, is schematically represented as

$$M+A\to {\left(MA\right)}^{*}$$

(135)

dissociations eventually with a long lifetime, where a state of a target nuclei, atom, or molecule (denoted as M) is pumped to the higher levels by borrowing the kinetic energy of an incident particle A. The Feshbach resonances and their overlapping set the foundation of Intramolecular Vibrational energy Redistribution (IVR) in the unimolecular dissociation dynamics, as discussed in Section 3.2. Thus, the dissociation channels are temporarily closed by a temporal compound state ${\left(MA\right)}^{*}$ surviving a long life-time. Both the Pauling and Feshbach resonances can be regarded basically as events on a single adiabatic potential surface.

During the continual state-mixing in the diffusive motion in B${}_{12}$, the electronic and nuclear energies are mutually exchanged through the nonadiabatic interactions. This phenomenon can be regarded as another kind of quantum mechanical resonance. The present dynamics found in B${}_{12}$ may be termed as hyper-resonance (resonance among the resonant states), since each adiabatic state is already in a Feshbach and/or Pauling resonance state. We refer to such a very long life-time tentative chemical bonding as dynamical bonding [165].

7.5. Intra-Molecular Nonadiabatic Electronic Energy Redistribution

7.5.1. Huge Inflation of Phase-Space Volume

Despite the very high energy of those excited-state boron clusters, they do not readily break apart, since the dissociation channels of them tend to be closed by very frequent nonadiabatic transitions. This is a natural consequence of hyper-resonance. We then survey a possible mathematical ground of the long life-time in terms of a drastically simplified statistical theory based on the phase space theory, originally proposed by J. C. Light [179]. This theory has been extensively applied to and modified before in the context of cluster dynamics [180,181,182], Calvo [183,184,185,186,187], and Fujii [65,66].

We consider only the simplest dissociation of B${}_{12}$ cluster in which only one atom is evaporated from N-body cluster ($N=12$). The phase space theory under no external field predicts the dissociation rate k as

$$k=\frac{W_{N-1}\left(E\right)}{{\Omega }_N\left(E\right)},$$

(136)

where ${\Omega }_N\left(E\right)$ is the density of the states of the N-particle cluster at the total energy E, and $W_{N-1}\left(E\right)$ denotes the total flux induced by the dissociation of a single atom [179]. As for a single adiabatic potential energy surface, say $V_I\left(\mathbf{R}\right)$ at a nuclear configuration $\mathbf{R}$, ${\Omega }_N\left(E\right)$ is approximately estimated with the expression

$${\Omega }_N^{\left(I\right)}\left(E\right)={\hslash }^{-3N}\int d\mathbf{R}d\mathbf{P}\delta \left(E-H_I(\mathbf{R},\mathbf{P})\right),$$

(137)

where

$$H_I(\mathbf{R},\mathbf{P})=\frac{\mathbf{P}·\mathbf{P}}{2}+V_I\left(\mathbf{R}\right),$$

(138)

with the obvious notation of the total nuclear momentum $\mathbf{P}$. We assume that the total energy is high above the potential barrier of the dissociation channels and that there are not significant centrifugal barriers for the present non-rotating clusters. In such systems, the so-called bottle-neck volume, say $\partial {\Omega }_N^{\left(I\right)},$ of the phase-space volume ${\Omega }_N^{\left(I\right)}$ to the dissociation channel cannot be appropriately defined. In place of resorting to $\partial {\Omega }_N^{\left(I\right)}$ therefore, we use the phase space volume of the daughter cluster of $N-1$ particles, denoted as ${\Omega }_{N-1}^{\left(I\right)}(E-\Delta E)$, where $\Delta E$ is the energy for the dissociating particle to take away with the relative velocity $\overrightarrow{v}$. The flux is approximately expressed as

$$W_{N-1}^{\left(I\right)}\left(E\right)=\int {\Omega }_{N-1}^{\left(I\right)}(E-\Delta E)\overrightarrow{v}\left(\Delta E\right)·d\overrightarrow{v}.$$

(139)

(The reduced mass is set to be unity for simplicity.) Thus, the critical quantity for the statistical rate of the dissociation probability turns out to be the ratio

$$\frac{{\Omega }_{N-1}^{\left(I\right)}(E-\Delta E)}{{\Omega }_N^{\left(I\right)}\left(E\right)}$$

(140)

for a given $\Delta E$.

We further suppose that the phase space of a system under study has the mixing property and is completely statistical due to chaos. Then we may regard the quantity

$${\Omega }_N^{\left(I\right)}\left(E\right)-{\Omega }_{N-1}^{\left(I\right)}(E-\Delta E)$$

(141)

as a phase space volume for trajectories to “tentatively” remain bound as N-body cluster. Then the related entropy $s^{\left(I\right)}(E,\Delta E)$ can then be defined as

$$s^{\left(I\right)}(E,\Delta E)=-log\left(\frac{{\Omega }_N^{\left(I\right)}\left(E\right)-{\Omega }_{N-1}^{\left(I\right)}(E-\Delta E)}{\Delta \Omega }\right),$$

(142)

where $\Delta \Omega$ represents the size of a unit cell of the phase space volume, such as ${\hslash }^{3N}$, which is used to count the number of cells in the volume $\left[{\Omega }_N^{\left(I\right)}\left(E\right)-{\Omega }_{N-1}^{\left(I\right)}(E-\Delta E)\right]$.

Although the state mixing is caused by the interaction

$$-i\hslash \sum _{k=1}^{3N}{\dot{R}}_k\left(t\right)X_{IJ}^k\left(\mathbf{R}\left(t\right)\right)$$

(143)

and/or $H_{IJ}^{\left(\mathrm{el}\right)}\left(\mathbf{R}\left(t\right)\right)$ as in Equation (109), we herein assume for simplicity that the transition among the electronic states occur freely. Then the entropy confined in the individual state $s^{\left(I\right)}(E,\Delta E)$ should be represented by the sum of them as

$$\begin{array}{c}{s}^{(max)}(E,\Delta E)=\sum _I{s}^{\left(I\right)}(E,\Delta E)\hfill \\ =-\sum _{I}log\left(\frac{{\Omega }_N^{\left(I\right)}\left(E\right)-{\Omega }_{N-1}^{\left(I\right)}(E-\Delta E)}{\Delta \Omega }\right)\hfill \end{array}$$

(144)

In such an extreme case the connected phase-space volume within $\Delta E$ is effectively extended and can be approximated as

$$\prod _I\left[{\Omega }_N^{\left(I\right)}\left(E\right)-{\Omega }_{N-1}^{\left(I\right)}(E-\Delta E)\right].$$

(145)

On the other hand, the phase space volume directly connected to the dissociation channel should be represented as the simple sum of each as

$$\sum _I{\Omega }_{N-1}^{\left(I\right)}(E-\Delta E).$$

(146)

Hence the the dissociation probability can be dramatically lowered from the order of Equation (140) to

$$\frac{\sum _I{\Omega }_{N-1}^{\left(I\right)}(E-\Delta E)}{\prod _J\left[{\Omega }_N^{\left(J\right)}\left(E\right)-{\Omega }_{N-1}^{\left(J\right)}(E-\Delta E)\right]},$$

(147)

indicating that the relative phase-space volume towards the dissociation is by far smaller than the phase space volume of the original N-body cluster, and the probability for the cluster to come across the dissociation channel is extremely small. Consequently, the life-time of the clusters is expected to be far longer than in case of no nonadiabatic transitions.

7.5.2. Intra-Molecular Nonadiabatic Electronic Energy Redistribution (INEER): Nonadiabatic Interaction to Close Dissociation Channels

Nonadiabatic electronic-state chaos thus inflates the volume of phase-space of nuclear motion of clusters to a very large extent. Due to this inflation, the molecular states lose chances to reach one of their dissociation channels. That is, a dissociating molecule can be brought back to the bound state by nonadiabatic interactions, which would dissociate otherwise. Here is an example of such an event [165]. In Figure 20, a B${}_{12}$ of initial hindered geometry is prepared, which is close to a dissociation channel. (Such dissociation channels can be readily found by colliding a B${}_{11}$ cluster and a B atom) In panel (a) is shown a trajectory for the 47th adiabatic state at this geometry to run on its adiabatic potential surface (no nonadiabatic interaction). After 110 fs of running, one of the atoms leaves the cluster leaving B${}_{11}$ cluster behind. Likewise, the same 47th adiabatic state in panel (b) begins to run from the same geometry with the same momentum but with the nonadiabatic interactions explicitly included in the electron dynamics. The trajectory clearly demonstrates that an atom being about to leave the cluster is brought back again to the B${}_{12}$, thus evidencing that the dissociation channel is closed. This phenomenon is what we call intra-molecular nonadiabatic electronic energy redistribution (INEER) due to hyper-resonance [188]. This makes a contrast to intramolecular vibrational energy redistribution (IVR) caused by classical chaos on a single potential basin. INEER is a crucial consequence of the intense chaos due to the continual nonadiabatic electronic-state mixing.

8. Concluding Remarks

To review “quantum chaos” in the scope of chemical dynamics, we herein have tracked the following subjects. We began with full quantum dynamics, wavepackets and eigenfunctions localized in the quasi-separatrix in two-dimensional potential function. This study was made in search for a novel molecular vibrational mode. Indeed, in this small dimensional system lie several vibrational modes, which give birth to very long-time wandering motions in the quasi-separatrix. A pair of the eigenstates of the dynamical tunneling of the second kind have been identified.

Semiclassical theory and energy quantization thereof, which is indispensable in the study of molecular vibration were next reviewed. We took a wavepacket approach in studying semiclassical dynamics rather than the direct semiclassical reduction of the Feynman kernel. With this approach, we discussed the mechanism of energy quantization for not only making the peaks but suppressing the off-resonant components. One of the most annoying parts in semiclassical kernel and related quantities is the divergence at caustics and/or turning points, which casts an inherent difficulty to the Gutzwiller trace formula. To avoid the difficulty, we have shown that the amplitude-free energy quantization actually works quite well. Chaotic paths that appear only through quantum tunneling have been shown along with a semiclassical tunneling formalism based on the notion of the parity of motion.

One of the main aims of the present review has been to guide the readers to chaos in molecular nonadiabatic dynamics, which is induced by the kinematic couplings between nuclei and electrons. The entanglement between electrons and nuclei is one of the most important subjects in general chemical science. We have emphasized the critical roles of “bifurcation and merging” of quantum wavefunctions. The most characteristic molecular chaos is found in nonadiabatic electron dynamics in the densely quasi-degenerate electronic-state manifold, every state of which mixes together with other members of the manifold. We have shown that the dynamical properties of electronic-state chaos can be measured in terms of the quantities related to diffusion and memory loss in the Hilbert space, turbulent flow of electron current, and so on. The chaos of the electronic excited state should become more and more vital in future molecular science, which will proceed far beyond the Born-Oppenheimer paradigm of the current theoretical framework of chemistry.

Some of the characteristics found for molecular quantum chaos are summarized below:

The basic mechanism of chaos in nonadiabatic wavepacket dynamics is “bifurcation and merging”, which is a quantum mechanical synonym of “stretching and folding” in the classical baker’s (or the Smale horse-shoe) transformation. Indeed repeated occurrences of “bifurcation and merging” of quantum wavepackets result in genuine chaos in quantum mechanics not only in energy spectra but more importantly in the dynamics of wavepackets. We would like to particularly emphasize again that the bifurcation and merging of wavepackets are indeed physically observable. Interestingly, “bifurcation (branching) and merging of the reaction tubes” sets also a geometrical foundation of classical chaos in structural isomerization dynamics. Extension of the tube geometry to quantum mechanics has not yet been explored.

“Loss of memory” is one of the central concepts to characterize chaos in general. The positive Lyapunov exponent in classical systems indicates an exponential separation of nearby orbits, leading to the loss of the memory of their initial stage. In our quantum chaos, the exponential decay of the memory of an electronic state in the time-propagation against the nonadiabatic (kinematic) interactions can take place in the Hilbert space. Again, we have observed that classical chaos in the structural isomerization dynamics of liquid-like clusters is associated with a characteristic feature of memory loss, which leads to a random yet statistical appearance of the isomers in the time series.

The advent of the ultrafast pulse laser has revolutionized molecular science since the 1970s. The femtosecond laser chemistry [189] to date has made it possible to real-time track the molecular deformation and reactions, an example of which can be seen in Figure 11 of this review. We are now in the age of technical development and applications of attosecond pulse lasers, the width of which is compatible with the time scales of electronic motions in molecules. This situation has been driving us in nonadiabatic electron wavepacket dynamics in Section 7. Nevertheless, we also call for attention for the study of extremely long-time spectroscopy. As we noted in this article, the structural isomerization of not only clusters but also large molecules like proteins can wander from one shape (potential valley, potential basin) to another taking a very long time at each local structure. These many-dimensional dynamics are often accompanied by large dimensional chaotic spaces, which makes the entire dynamics even longer and more complicated. It is extremely challenging therefore to figure out spectroscopic means to capture the essential feature of the long-time dynamics. The time-dependent spectroscopy, which was proposed in Section 4 and discussed briefly in the section of semiclassical mechanics Section 5, is just one primitive idea in an attempt to cope with this matter. A parallel computation scheme with many different initial conditions should be a good candidate to overcome long-time quantum dynamics in ergodic phase space, if a method to appropriately connect the pieces of spectra taken individually is devised.

A final remark before closing. The essential nature of classical and quantum chaos may lie in mathematical beauty and intricate geometry behind the complicated-looking phenomena. However, the complexity of molecules often makes such elegant mathematical structures obscure and even less mighty. Moreover, We are afraid that sticking to the too rigorous definitions of quantum chaos would isolate itself from other natural science. Therefore, I would like to conclude this review by stressing that more studies on quantum chaos made mathematically easy and conceptually useful are demanded than ever, since there are so many unknowns left (not only) in the molecular phenomena, notions, and laws that need an understanding of quantum chaos.