Density of Avoided Crossings and Diabatic Representation

Abstract

Electronic structure theory describes the properties of solids using Bloch states that correspond to highly symmetrical nuclear configurations. However, nuclear thermal motion destroys translation symmetry. Here, we describe two approaches relevant to the time evolution of electronic states in the presence of thermal fluctuations. On the one hand, the direct solution of the time-dependent Schrodinger equation for a tight-binding model reveals the diabatic nature of time evolution. On the other hand, because of random nuclear configurations, the electronic Hamiltonian falls into the class of random matrices, which have universal features in their energy spectra. In the end, we discuss combining two approaches to obtain new insights into the influence of thermal fluctuations on electronic states.

1. Introduction

Bloch’s theorem and translational symmetry are two of the main bricks in the foundation of solid state physics. Their use in crystal lattices makes an implicit assumption, not emphasized in the texts we are aware of. Bloch’s theorem is mathematically correct of course for a fixed, periodic potential, but real nuclei move under the influence of thermal fluctuations, so that they have random configuration without translational symmetry at any moment in time. This is clearest in the coherent state representation of the lattice, which assigns a mean position and momentum to each atom [1].

The success of solid state theory suggests there should be a physical justification for using electronic states corresponding to fixed highly symmetrical positions. Recently, it has been found that for graphene $\pi$-bands, thermal fluctuations turn out to be too fast for the adiabatic Born–Oppenheimer approximation (ABO) to be valid, so that the other limit—the diabatic limit (DBO)—holds [2].

If we define the Born–Oppenheimer approximation as supposing that the electron wavefunction returns to its starting form if the nuclei return to their starting positions after an arbitrary journey, this leaves the possibility of adiabatic and diabatic motion of the electron along the journey. We allow the caveat of possible Berry phases within the Born–Oppenheimer realm.

During thermal fluctuations, atomic positions have a random component for any arbitrary moment in time. The electronic Hamiltonian, which is a function of atomic positions, acquires this random character. Properties of simple prototype random matrix models are well studied in random matrix theory; some of them are universal and can be extended to more general cases. Here, we are interested in ensembles of parameter-dependent random matrices. In addition to level statistics, statistics of singularities in their spectrum including degeneracies and avoided crossings is known. We believe these results can be used to characterize electronic time evolution in the presence of phonons.

Indeed, the ABO breaks down when degeneracies occur; for the case of avoided crossings, Landau and Zener found that the probability of adiabatic dynamics depends on level spacing and sensitivity to perturbations. Random matrix theory combined with Landau–Zener theory makes it possible to characterize the effects of thermal fluctuations on electronic states.

The paper is organized as follows. In Section 2.1, we describe parameter-dependent random matrices; in Section 2.2 and Section 2.3, we discuss statistical properties of singularities in their spectrum, namely, conical intersections and avoided crossings. In Section 3, we describe results on graphene thermal fluctuations that suggest diabatic representation. In Section 4, we sketch future directions for describing electronic dynamics in the presence of phonons using random matrix theory.

2. Singularities in the Spectra of Random Matrices

2.1. Parameter-Dependent Random Matrices

Random matrices provide a pathway to universal behavior for complex quantum systems [3]. One example of such a system is the chaotic dynamics of a kicked top [4,5]. In this paper, we consider only Gaussian ensembles, which are denoted by the Dyson index $\beta$ [3], where $\beta =1$ for GOE, $\beta =2$ for GUE, and $\beta =4$ for GSE. Corresponding matrices $H_{\beta }$ can be represented as a combination of symmetrical S and antisymmetrical A matrices, whose elements are independent Gaussian distributed random variables with mean zero and variance $(1\pm {\delta }_{ij})/\beta$:

$$H_1=S,H_2=S+iA,H_4=S{e}_0+A^1{e}_1+A^2{e}_2+A^3{e}_3$$

(1)

where $e_i$ are 2 × 2 matrix representations of quaternion algebra bases:

$$e_0=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),e_1=\left(\begin{array}{cc}0& i\\ i& 0\end{array}\right),e_2=\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right),e_3=\left(\begin{array}{cc}i& 0\\ 0& -i\end{array}\right)$$

(2)

and the direct product is assumed for $A^n{e}_n$.

We study the cases when a Hamiltonian $\widehat{H}$ depends on a set of parameters, which, in turn, depend on time. The idea of parameter-dependent random matrices first appears in [6] to quantify the dissipation rate of a driven complex quantum system. We use the following parametrization $H(x_1,x_2,\cdots ,x_n)$ suggested by Wilkinson and Austin [7]

$$H(x_1,x_2,\cdots )=cos{x}_1{H}_1+sin{x}_1{H}_2+cos{x}_2{H}_3+sin{x}_2{H}_4+\cdots .$$

(3)

where matrices $H_1$, $H_2$, ⋯ belong to the same ensemble and $\mathbf{x}=(x_1,x_2,\cdots ,x_n)$. These matrices may also belong to different ensembles to study symmetry-breaking perturbations as Fritz Haake suggested [8,9,10].

This parametrization has a number of advantages. Distribution of matrix elements is the same for all $\mathbf{x}$, and matrices H and $\partial H/\partial x_j$ are independent, therefore leading to rather simple analytical averaging within the degenerate perturbation theory; the extension to the cases when $\partial H/\partial X_1$ and $\partial H/\partial X_2$ have non-zero correlation was made by Wilkinson [11,12]. Equation (3) provides an ergodic property [13]—averaging over energy levels leads to the same level spacing distribution as averaging over parameter space for a given pair of levels—enabling, for example, geometrical considerations to determine scaling laws for the density of avoided crossings [11].

An alternative to Equation (3) is the linear parametrization

$$H(x_1,x_2,\cdots )=H+H_1{x}_1+H_2{x}_2+\cdots .$$

(4)

However, now singularities in the spectrum happen only in the vicinity of $\mathbf{x}=0$, and the distribution of eigenvalue curvatures ${\partial }^2{ϵ}_i/\partial x_j^2$ tends to zero for $x_j\to \infty$.

2.2. Geometrical Properties of Conical Intersections

For a Hamiltonian $\widehat{H}$ to have a doubly degenerate energy level, it is necessary to satisfy two independent conditions. We reproduce here the argument of Teller [14]. Let us assume that all but two of the $\widehat{H}$ eigenstates have been found. Expressing $\widehat{H}$ in the basis formed by these eigenstates and two arbitrary states i and $j,$ forming altogether a complete set of functions, the condition for the degeneracy is

$$H_{ij}=0,H_{ii}=H_{jj}.$$

(5)

For ensemble $\beta$, it leads to a system of $\beta +1$ linear equations. In two-dimensional parameter space, the levels are degenerate in isolated points called conical intersections.

The geometrical properties of conical intersections have been established in a number of papers. Longuett-Higgins found that if a wavefunction changes sign when transported adiabatically round a given loop in parameter space, then the state must become degenerate with another one at some point within the loop [15]. This topological property is one consequence of the Berry phase [16].

The distribution of conical intersections in the parameter space is characterized by density ${\mathcal{D}}_n^{\left(\mathrm{ci}\right)}$, which is the number of points, in which energy levels n and $n+1$ are degenerate, per unit area. Wilkinson and Austin proposed the idea for its analytical derivation [11]. For a small element $dA$ of the parameter space, the probability $dp$ to find a conical intersection is ${\mathcal{D}}_n^{\left(\mathrm{ci}\right)}dA$. At the same time, $dp$ is determined by the distribution $P\left[R\right]$ of the distances R from a random point to the nearest conical intersection

$$dp=P\left[R\right]dR={\mathcal{D}}_n^{\left(\mathrm{ci}\right)}dA.$$

(6)

Hence, the problem is reduced to analytical calculation of $P\left[R\right]$. It is possible to express R by using the degenerate perturbation theory. Averaging over the corresponding ensemble leads to well-known integrals [17,18,19]. ${\mathcal{D}}_n^{\left(\mathrm{ci}\right)}$ depends on the density of states ${\rho }_n=\rho \left(E_n\right)$ and the variance of the off-diagonal matrix elements ${\sigma }^2$ and has the general form $C_{\beta }{\left(\rho \sigma \right)}^{\beta +1}$:

$${\mathcal{D}}_n^{\left(\mathrm{ci}\right)}=\frac{\pi }{3}{\rho }_n^2{\sigma }^2,{\mathcal{D}}_n^{\left(\mathrm{ci}\right)}=\frac{2\sqrt{\pi }}{3}{\rho }_n^3{\sigma }^3,{\mathcal{D}}_n^{\left(\mathrm{ci}\right)}=\frac{16\sqrt{2}{\pi }^{3/2}}{45}{\rho }_n^5{\sigma }^5$$

(7)

for GOE [11], GUE, and GSE [20], respectively. We note that prefactors $C_{\beta }$ were found only after corresponding scaling had been revealed in numerical calculations [21,22]. Such an order of discoveries is typical in this field, for example, the exact distribution for energy level curvatures for GUE was found only after numerical experiments [23,24,25,26].

By putting the Wigner semicircle distribution [3] in Equation (7) and integrating the resulting expression over energy E, one gets the density of degeneracies summed over all pairs of levels. Multiplied by the volume of the parameter space, this gives the total number of degeneracies $\mathcal{N}$. For $M\times M$ GOE matrices parametrized as Equation (3), Wilkinson and Austin found [11]

$$\mathcal{N}=\frac{\pi }{2}{M}^2.$$

(8)

For GOE-parametrized matrices, the potential energy surfaces $E_n(x_1,x_2)$ and $E_{n+1}(x_1,x_2)$ form a cone near conical intersections. The corresponding energy level lines are ellipses. Wilkinson and Austin [11] calculated the distribution of eccentricities and cone slopes. We illustrate them in Figure 1, which is taken from their paper. We note that ellipses are randomly oriented and mostly elongated in accordance with the distribution

$$P\left[e\right]={\left(\frac{2e}{2-e^2}\right)}^3$$

(9)

coming from the joint probability distribution of slopes of the cone s and eccentricity e:

$$P[t,d]=\frac{d}{256{\sigma }^6}{e}^{-t/8{\sigma }^2},$$

(10)

where d and t unambiguously determine e and s [11].

One might expect that, similar to the spacings between adjacent energy levels, conical intersections between the same pair of levels would tend to avoid each other. It turns out that this is not the case for the GOE-parametrized model, Equation (3); there, degeneracies can be considered as randomly distributed points [11]. However, as one can see on Figure 1, their eccentricities and orientations are correlated.

2.3. Geometrical Properties of Avoided Crossings

Consider the case when the Hamiltonian depends on only one parameter x. This is one of a wide range of systems with an arbitrary large number of time-dependent parameters. The noncrossing rule follows from Equation (5) [27]; when x is varied, two energy levels never cross. Still, they can approach each other with an arbitrary small spacing. A local minimum in the distance between adjacent energy levels is called an avoided crossing. We refer to the difference between energy levels at the minimum as a gap. The distribution of crossing gaps is known for the GOE [28,29,30], GUE, and GSE [29,30]. Surprisingly, the ratio of mean crossing to mean spacing is not universal and varies from 0.42 to 0.52 for different realizations of the same ensemble $\beta$ [30]. Furthermore, known as level repulsion, the noncrossing rule has been generalized to non-hermitian Hamiltonians that describe dissipative systems [31].

Energy profiles near an avoided crossing have the hyperbolic form

$$E_{n,n+1}=B(x-x_0)\pm \frac{1}{2}\sqrt{{ϵ}^2+A^2{(x-x_0)}^2}$$

(11)

which is a solution for the Hamiltonian

$$H\left(x\right)=\left(\begin{array}{cc}(B+A/2)(x-x_0)& ϵ/2\\ ϵ/2& (B-A/2)(x-x_0)\end{array}\right).$$

(12)

When the Hamiltonian depends on a set of time-dependent parameters, such avoided crossings appear in the vicinity of degeneracies (the conical intersection geometry). Using this property together with ergodicity, Wilkinson and Austin have shown in a simple manner that the number of avoided crossings with gaps less than $\mathsf{\Delta }$ is proportional to $\mathsf{\Delta }$ in the limit $\mathsf{\Delta }\to 0$ [11].

The density of avoided crossings ${\mathcal{D}}_n^{\left(\mathrm{ac}\right)}$ is the number of avoided crossings between n and $n+1$ energy levels per unit length as we traverse a line by varying x. The density ${\mathcal{D}}_{\beta }^{\left(\mathrm{ac}\right)}(A,B,ϵ)$ of avoided crossings with gap $ϵ$ and parameters A and B of the hyperbola Equation (11) was found for $\beta$-parameterized models [6,20]

$${\mathcal{D}}_{\beta }^{\left(\mathrm{ac}\right)}(A,B,ϵ)=C_{\beta }^{\left(\mathrm{ac}\right)}\frac{{\rho }^{\beta +1}}{{\sigma }^{\beta +1}}{A}^{\beta +1}{ϵ}^{\beta -1}exp\left(-\frac{\beta A^2}{8{\sigma }^2}\right)exp\left(-\frac{\beta B^2}{2{\sigma }^2}\right)\frac{\sqrt{\beta }}{\sqrt{2\pi }\sigma }dAdBdϵ$$

(13)

with the prefactors [20]

$$C_1^{\left(\mathrm{ac}\right)}=\frac{\pi }{24},C_2^{\left(\mathrm{ac}\right)}=\frac{{\pi }^{3/2}}{12},C_4^{\left(\mathrm{ac}\right)}=\frac{8{\pi }^{7/2}}{135\sqrt{2}}.$$

(14)

3. Thermal Fluctuations in Solids

3.1. Born–Oppenheimer Approximation

All information about the properties of a system containing $N_n$ nuclei and $N_e$ electrons is contained in the eigenstates of the many-body Hamiltonian

$$\widehat{H}=\sum _{I=1}^{N_n}\frac{{\mathbf{P}}_I^2}{2M_I}+\sum _{i=1}^{N_e}\frac{{\mathbf{p}}_i^2}{2m}+\sum _{1\le I<J}^{N_n}\frac{Z_I{Z}_J{e}^2}{|{\mathbf{R}}_I-{\mathbf{R}}_J|}+\sum _{1\le i<j}^{N_e}\frac{e^2}{|{\mathbf{r}}_i-{\mathbf{r}}_j|}-\sum _{I=1}^{N_n}\sum _{i=1}^{N_e}\frac{Z_I{e}^2}{|{\mathbf{r}}_i-{\mathbf{R}}_I|},$$

(15)

where m is the mass of an electron, e is the charge of an electron, ${\mathbf{r}}_i$ and ${\mathbf{p}}_i$ are the position and the momentum of the i-th electron, respectively, $M_I$ is the mass of I-th nucleus, $Z_I$ is the charge of I-th nucleus, and ${\mathbf{R}}_I$ and ${\mathbf{P}}_I$ are the position and the momentum of I-th nucleus, respectively. The Hamiltonian describes kinetic energy of electrons and nuclei, Coulomb interaction between nuclei, Coulomb interaction between electrons, and finally Coulomb interaction between electrons and nuclei.

Eigenstates of Equation (15) can be found analytically only for the hydrogen atom. Numerical solutions are possible for systems containing several atoms and electrons. If there are M points in grids representing the positions of electrons and nuclei, $M^{3(N_n+N_e)}$ real numbers would be required to store one state of the system. The exponential complexity quickly hits the memory limits of modern machines, so approximations are necessary.

The Born–Oppenheimer approximation [32,33] is the first simplification of the problem. It separates nuclear and electronic degrees of freedom. The separation is justified by the fact that nuclei are much heavier than electrons, and as a result electrons adjust to nuclear motion.

Within the Born–Oppenheimer approximation, atomic positions ${\mathbf{R}}_I$ are fixed to give the electronic Hamiltonian

$${\widehat{H}}_e({\mathbf{R}}_1,\cdots ,{\mathbf{R}}_{N_n})=\sum _{i=1}^{N_e}\frac{{\mathbf{p}}_i^2}{2m}+\sum _{1\le I<J}^{N_n}\frac{Z_I{Z}_J{e}^2}{|{\mathbf{R}}_I-{\mathbf{R}}_J|}+\sum _{1\le i<j}^{N_e}\frac{e^2}{|{\mathbf{r}}_i-{\mathbf{r}}_j|}-\sum _{I=1}^{N_n}\sum _{i=1}^{N_e}\frac{Z_I{e}^2}{|{\mathbf{r}}_i-{\mathbf{R}}_I|}.$$

(16)

The electronic Hamiltonian is then diagonalized to obtain eigenvalues

$${\widehat{H}}_{el}({\mathbf{R}}_1,\cdots ,{\mathbf{R}}_{N_n})|{\psi }_n\rangle =E_n({\mathbf{R}}_1,\cdots ,{\mathbf{R}}_{N_n})|{\psi }_n\rangle ,$$

(17)

which are functions of nuclear positions. These functions are then regarded as potentials in which nuclei move

$$H_n=\sum _{I=1}^{N_n}\frac{{\mathbf{P}}_I^2}{2M_I}+E_n({\mathbf{R}}_1,\cdots ,{\mathbf{R}}_{N_n}).$$

(18)

This Hamiltonian includes only nuclear degrees of freedom. It can be used to obtain both quantum and classical solutions. In the classical treatment, the motion of a nuclei is approximated as

$$M_I{\ddot{\mathbf{R}}}_I=-\frac{\partial E_n({\mathbf{R}}_1,\cdots ,{\mathbf{R}}_{N_n})}{{\mathbf{R}}_I}.$$

(19)

In many cases, the collection of nuclei can be represented as classical harmonic oscillators, which is justified by the Schrodinger correspondence principle [34,35].

A more restrictive simplification is the clamped nuclei approximation [36], which assumes that nuclei have fixed positions corresponding to the equilibrium structure. The original problem is then reduced to finding the ground state of the electronic Hamiltonian for a fixed set of parameters, Equation (16). This problem still has exponential complexity, since $M^{3N_e}$ real numbers are required to store an electronic state. This complexity is overcome by a single particle picture, in which an electron moves independently of all other electrons. In this picture, the interaction of a single electron with other electrons is represented as the interaction of an electron with some average external potential due to the remaining electrons. Density functional theory provides the most accurate description of electron–electron interactions in this way [37]. Another method commonly used to solve the electronic Hamiltonian is the tight-binding method [38,39].

3.2. Supercell Technique

Within the Born–Oppenheimer approximation, clamped nuclei approximation, and single-particle picture, the microscopic properties of a solid can be found by solving the Schrodinger equation $\widehat{H}|{\psi }_n\rangle =E_n|{\psi }_n\rangle$ for a single electron in the external potential set by equilibrium atomic positions $\left\{{\mathbf{R}}_I\right\}$. The problem is still too complex, since solids consist of an infinite number of atoms and electrons, making the electronic Hamiltonian a function of an infinite number of parameters. The use of translational symmetry together with periodic boundary conditions reduces the problem to a limited number of electrons inside a primitive cell, and the Hamiltonian is hence a function of a finite set of parameters—vectors that define the primitive unit cell and atomic coordinates inside the cell. Furthermore, Bloch’s theorem gives eigenstates of the resulting Hamiltonian. Since atomic positions are fixed in this description, this solution provides only the static properties of a solid. A central property is the total energy of the system as a function of its structure. Minimizing it, one can calculate equilibrium lattice constants, the bulk modulus, and the equation of state $E\left(V\right)$ or $P\left(V\right)$, which is now a standard methodology in computational materials science [37].

In real crystals, of course, atoms are constantly fluctuating around equilibrium positions. At an arbitrary moment of time, atomic positions are random and there is no translational symmetry. Bloch’s theorem is no longer applicable. The problem becomes intractable since the Hamiltonian is again a function of an infinite number of parameters. In some cases, the complexity can still be harnessed by using the supercell technique [40,41]. In this technique, a supercell is constructed of N repeating primitive cells and the atoms inside it are allowed to move freely. Periodic boundary conditions are then imposed on this structure. It therefore becomes an infinite solid, where the supercell plays the role of a primitive cell. Translational symmetry is restored and Bloch’s theorem can be applied. The atomic motion inside the supercell allows access to dynamical properties such as phonons or relaxation times. Strictly speaking, the thus calculated dynamical properties would be different from those of the real solid; they would rather represent the properties of a large “molecule” consisting of N primitive cells whose motion is repeated throughout the space. However, upon an increase in N, the calculated properties would more closely resemble those of the real solid; in other words, by increasing N, it is possible to achieve numerical convergence.

Surprisingly, in many cases, the size of the supercell N necessary for numerical convergence is not exceedingly high in terms of numerical calculations. Physically, this means that in such cases interactions inside a solid quickly decay over distance. For example, for $\alpha$-quartz, phononic dispersion can be accurately calculated by displacing atoms in a supercell consisting of only 27 unit cells. The reason is that force constants between atoms in the unit cells separated by more than three lattice spacings are weak [42]. In general, of course, the value of N necessary for convergence depends both on the properties and the solid, and it is possible that numerical convergence may not be achieved due to the limited computational resources of modern computers.

The dynamical properties of solids calculated using the supercell technique also rely on the assumption that electronic dynamics has a Born–Oppenheimer adiabatic character. The electronic Hamiltonian for the supercell is set by positions of the atoms inside it. Supposing the atoms are moving classically, the Hamiltonian becomes time dependent. Therefore, to describe electronic dynamics $|\psi \left(t\right)\rangle$, it is necessary to solve the time-dependent Schrodinger equation starting from the n-th eigenstate

$$\begin{array}{cc}\hfill i\hslash \frac{\partial |\psi \left(t\right)\rangle }{\partial t}& =\widehat{H}\left(t\right)|\varphi \left(t\right)\rangle ,\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill |\psi \left(0\right)\rangle & =|{\psi }_n\left(0\right)\rangle .\hfill \end{array}$$

(21)

The time-dependent wavefunction $|\psi \left(t\right)\rangle$ often used in supercell calculations corresponds to the n-th eigenstate of the Hamiltonian $\widehat{H}\left(t\right)$ at time t

$$\begin{array}{cc}\hfill \widehat{H}\left(t\right)|{\psi }_n\left(t\right)\rangle & =E_n\left(t\right)|{\psi }_n\left(t\right)\rangle ,\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill |\psi \left(t\right)\rangle & =|{\psi }_n\left(t\right)\rangle .\hfill \end{array}$$

(23)

It is different from the wavefunction that one obtains by solving the time-dependent Schrödinger equation, but according to the adiabatic theorem, they coincide if the Hamiltonian changes slowly enough (the adiabatic Born–Oppenheimer approximation (ABO)).

The adiabatic assumption can be checked by comparing the time-dependent solution with the adiabatic one. Alternatively, the adiabatic assumption can be indirectly validated by the fact that calculated dynamical observables coincide with experimental values. In most applications, there is no need for direct validation, since, for solids, theoretical predictions turn out to correctly describe the experimental behavior.

The case of graphene is unique in this regard. Graphene is a two-dimensional honeycomb structure of carbon atoms. Each carbon atom has four valence orbitals, three of which hybridize to form strong in-plane $\sigma$-bonds between adjacent atoms. The fourth electron is in a p${}_z$ orbital perpendicular graphene’s plane. Out-of-plane p${}_z$ electrons form $\pi$-bonds that are significantly weaker than $\sigma$-bonds, and therefore their impact on atomic motion is weak. Although atomic motion is primarily determined by $\sigma$-bonds, electrons in $\sigma$-bonds do not contribute to electronic transport. Therefore, the fact that numerical calculations reproduce experimental values for atomic vibrations confirms the validity of the adiabatic assumption for $\sigma$-bonds only. One can check the validity of the adiabatic theorem for $\pi$-electrons by directly comparing time-dependent electronic wavefunctions to the adiabatic ones. This comparison was performed by Mohanty and Heller [2], which we describe below.

3.3. Simulating Graphene Thermal Fluctuations

To model the thermal fluctuations of graphene, the authors constructed a 4 × 4 supercell with 32 carbon atoms. To check that qualitative conclusions hold for larger sizes, they repeated calculations for a 5 × 10 supercell. Within a supercell, atoms were allowed to move freely under the harmonic $\sigma$-bonded force field, with atomic motion taken as a superposition of graphene phonons that satisfied the periodic boundary conditions of a given supercell. In turn, graphene phononic frequencies were computed by considering only force constants between adjacent atoms, and the corresponding value was taken from infrared spectroscopy of the C=C aromatic bond. The amplitude of each phonon corresponded to a temperature $T=300$ K. Eventually, atomic positions were predetermined for each moment of time according to

$$\mathbf{r}\left(t\right)=\sum _n\sqrt{\frac{2k_{B}T}{M{\omega }_n^2}}{\mathbf{q}}_{n}cos{\omega }_{n}t,$$

(24)

where ${\mathbf{q}}_n$ is the normal mode for the n-th phonon with frequency ${\omega }_n$. $\pi$-electrons do not directly influence this motion. This is a simplification, but it is expected that the qualitative conclusions hold since the stiff $\sigma$-bonds contribute most to the force.

At each moment of time, atomic positions determine the nearest-neighbor tight-binding Hamiltonian for $\pi$-electrons. The hopping integral between adjacent atomic sites depends on the distance $r_{ij}\left(t\right)$ between the atoms

$$H\left(t\right)=-\sum _{\langle ij\rangle }\beta exp[-\alpha r_{ij}\left(t\right)]({\widehat{c}}_i^{+}{\widehat{c}}_j+H.c.),$$

(25)

where $\langle ij\rangle$ stand for the indices of adjacent atoms.

The time-dependent Hamiltonian $\widehat{H}\left(t\right)$ is used to obtain two types of wavefunction. The first corresponds to the adiabatic solution $|{\psi }_n\left(t\right)\rangle$, Equation (22). Adiabatic diagonalization based on the instantaneous atomic positions provides the eigenenergies $E_n\left(t\right)$. The second type of wavefunction is the actual electronic wavefunction that is given by Equation (20).

To study the validity of the ABO, the true wavefunction $|\psi \left(t\right)\rangle$ is compared with the adiabatic solution $|{\psi }_n\left(t\right)\rangle$ by calculating the overlap between them

$$a_n\left(t\right)={\left|\langle {\psi }_n\left(t\right)|\psi \left(t\right)\rangle \right|}^2,$$

(26)

which is called the adiabatic correlation.

The adiabatic correlation reveals whether the adiabatic theorem, and consequently the ABO, is valid for this process. However, it does not say how strong the changes in the wavefunctions $|{\psi }_n\left(t\right)\rangle$ and $|\psi \left(t\right)\rangle$ are. For this, the overlap of the actual wavefunction with the initial wavefunction was calculated. The corresponding quantity is called the diabatic autocorrelation

$$d\left(t\right)={\left|\langle \psi \left(0\right)\left|\psi \right(t)\rangle \right|}^2.$$

(27)

3.4. Breakdown of the ABO for Thermally Fluctuating Graphene

Eigenenergies of the time-dependent Hamiltonian $E_n\left(t\right)$ for a 4 × 4 supercell, Equation (22), are shown in Figure 2a. Time is given in the units of the shortest vibrational period, which is 191 $\hslash /$Ry for the case of the 4 × 4 supercell. There are 32 eigenenergies because the supercell accommodates 32 atoms; for any moment of time, eigenvalues are not degenerate because atomic positions do not correspond to a highly symmetrical configuration including the initial moment of time $t=0$, Equation (24). Eigenenergies form clusters that are separated from each other by gaps of several eV.

The time evolution of the eigenenergies of one of the clusters is seen to fluctuate around the average values in Figure 2b. Adjacent levels constantly approach each other; however, spacings between them never go to zero. If the dynamics were adiabatic, $E_n\left(t\right)$ would describe how the expected values of the time-dependent states evolve. To reveal the character of the dynamics, adiabatic correlations and diabatic autocorrelations are calculated for the states starting from the ground state and from the $n=18$ state at $t=0$.

For the ground state, the adiabatic correlation remains one $a_1\left(t\right)=1$ during the simulation time equal to tens of the shortest vibrational period. The diabatic autocorrelation remains one $d\left(t\right)=1$ as well. Preservation of adiabatic correlation means that the solution of Equation (20) $|\psi \left(t\right)\rangle$ and the solution of Equation (22) $|{\psi }_1\left(t\right)\rangle$ coincide up to an arbitrary phase, $|\psi \left(t\right)\rangle =|{\psi }_1\left(t\right)\rangle$; therefore, the adiabatic theorem and consequently the ABO hold for time evolution of this particular state. Since the diabatic autocorrelation remains one as well, we conclude that the thermal motion of atoms has no effect on the ground state, and neither the adiabatic nor the true wavefunctions change.

The behavior of higher states that form the clusters is different. Figure 3a shows adiabatic correlations for the state $|\psi \left(t\right)\rangle$, starting from the n = 18 eigenstate at t = 0. In contrast to the ground state, the overlap probability with the corresponding adiabatic state sharply decays to zero during the first vibrational period. For the first four vibrational periods, it remains zero most of the time, occasionally sharply rising to $a_{18}\approx 0.5$ and then sharply decaying to zero again three times during the aforementioned time. In addition to this, we can also observe the overlap probability with other adiabatic states in the Figure. These states form a cluster of eigenenergies, Figure 2b. The overlap with each of these states is constantly changing, rising from zero to a maximum value and then decaying back to zero. The maximum value varies among states of the cluster and lies within the range from $a_{25}\approx 0.25$ for the $n=25$ adiabatic eigenstate to $a_{21}=0.8$ for the n = 21 adiabatic eigenstate. For certain states, for example, $n=21$, peaks in the adiabatic correlation merge, while for the others, peaks in the adiabatic correlation stand separately.

The diabatic autocorrelation for the $n=18$ state also has a different behavior to the one of the ground state, Figure 3b. In contrast to the ground state case, it does not remain one throughout the simulation. It gradually decays to zero within 60 periods of the shortest vibration. The decay time of the diabatic correlation is several orders of magnitude larger than the decay time of the adiabatic correlation. After dropping to zero, the diabatic correlation function grows.

From the quick decay of the adiabatic correlations, it is clear that in the presence of thermal fluctuations, the adiabatic theorem and the ABO do not hold for states that form clusters. The adiabatic theorem and the ABO are still valid for the ground state that is separated from all the others by gaps of the order of several eV. Such a sharp decay of the adiabatic correlation can be attributed to the near degeneracy of the states, suffering frequent avoided crossings. As we discuss in the next section, these avoided crossings indeed have parameters that lead to the breakdown of the adiabatic theorem.

Secondly, the fact that the diabatic autocorrelation decays on a timescale several orders larger than adiabatic correlations means that during the very first oscillations, electronic wavefunctions stay intact in the sense that the atomic orbital coefficients remain nearly unchanged as they are carried around by their respective nuclei. Since adiabatic correlations drop sharply to zero in the meantime, one can conclude that even the smallest changes in atomic positions lead to enormous changes in the adiabatic wavefunctions $|{\psi }_n\left(t\right)\rangle$ that form a cluster and that are solutions of Equation (22). This is in sharp contrast to the ground state, where both adiabatic and diabatic wavefunctions did not change during the whole simulation time. From the perspective of the zero-th order degenerate perturbation theory, almost degenerate states form a complete basis set to describe changes in adiabatic electronic wavefunctions induced by thermal motion of atoms. Quick decay of adiabatic correlations and preservation of the initial wavefunction means that even slight atomic displacements lead to strong mixing between initial adiabatic states. This is a characteristic feature of states that form clusters.

Finally, the decay of the diabatic autocorrelation at much larger times than adiabatic correlations means that electrons cannot adjust to nuclear motion. The ABO ansatz is based on this assumption. Therefore, it is not optimal for description of such dynamical processes. Qualitatively, one may rather say that electronic states are preserved during dynamics. Since time evolution has diabatic character, then another ansatz, which the authors called the diabatic Born–Oppenheimer approximation, would be more suitable than the ABO for a description of such processes.

Qualitatively, the diabatic character is in accordance with experiments. Indeed, electronic structure theory and Bloch functions are built assuming atoms form perfect symmetrical structures. The success of electronic structure theory indicates that during thermal fluctuations, electronic wavefunctions are preserved by obeying a diabatic solution near avoided crossings. Therefore, there is an implicit diabatic assumption in electronic structure theory, otherwise it would be impossible to exploit translational symmetry. The diabatic representation is that the ABO solution for the highest symmetry is used no matter what the actual positions of the nuclei are.

4. Discussion and Further Work

In this section, we suggest future directions in the study of time-dependent dynamics of electrons in the presence of phonons. Mohanty and Heller [2] have established its diabatic character by running time-dependent simulations. It is interesting, however, that the lowest frequency phonons have not been studied, owing to the the finite size of the supercell. They are slow; thus, they may be again in the adiabatic realm. It would be very constructive to determine this.

We now suggest a new perspective based on avoided crossings. We outline how one can potentially come to the same conclusions by using Landau–Zener theory in tandem with random matrix theory without running time-dependent calculations.

The idea of using Landau–Zener theory in tandem with random matrix theory belongs to Wilkinson [6]. They considered a driven finite-sized system and calculated the rate of increase in energy of the system using a distribution of gap sizes and the slopes of avoided crossings, similar to Equation (13). It was noted in reference [6] that this idea could be applied in the context of the Born–Oppenheimer approximation to characterize vibronic coupling in complex molecules.

Diabatic preservation of electronic character suggests a breakdown of the adiabatic approximation. Adiabaticity breaks down when the Hamiltonian changes too rapidly. One can estimate a timescale of adiabatic motion by using Heisenberg’s uncertainty principle; for a process with the minimal spacing between energy levels $\mathsf{\Delta }E$, adiabatic evolution can take place only for time $\tau$

$$\tau \ge \hslash /\mathsf{\Delta }E.$$

(28)

From Equation (28), it follows that avoided crossings or local minima in the distance between adjacent levels determine the timescale for adiabatic motion. We have considered their statistical properties in Section 2.3.

Figure 2b shows how energy levels vary over time for graphene thermal fluctuations. Each state experiences several avoided crossings during one cycle of the shortest vibrational mode. We conclude that the presence of multiple avoided crossings is an inherent feature of this motion, and therefore should be explicitly incorporated into theoretical treatment. Specifically, avoided crossings determine a characteristic energy gap, which in turn determines a timescale for adiabatic evolution. This timescale is to be compared with the vibrational period, which is primarily set by $\sigma$-bonds and only weakly depends on p${}_z$-orbitals.

For the case [2], the shortest vibrational period is $\tau \approx 191\hslash /$Ry, so the energy gap should be at least $\mathsf{\Delta }E=1/191\approx 0.005$ Ry. As one can see in Figure 2b, avoided crossings have a much smaller gap, leaving no hope for adiabatic motion. Indeed, in Figure 3a, the adiabatic autocorrelation of the 18th state drops to zero at $t=0.25$, which is exactly the first avoided crossing one can observe in Figure 2b.

Landau and Zener gave a more accurate estimation for adiabatic timescales [43]. They considered a two-level system with the Hamiltonian in the form of Equation (12) for the case $B=0,x_0=0$, and $x\equiv t$. The probability of adiabatic time evolution by the end of the motion is

$$p=1-exp\left(-\frac{\pi {ϵ}^2}{2A\hslash }\right).$$

(29)

The estimation of adiabatic $A_{\mathrm{adia}}$ based on the exact solution $A_{\mathrm{adia}}=\frac{\pi }{2}{ϵ}^2/\hslash$ is close to the one based on Heisenberg’s uncertainty principle, $A={ϵ}^2/\hslash$. In the case of graphene thermal fluctuations, the electronic structure of the p${}_z$ orbitals determines $ϵ$, while the temperature and $\sigma$-bonds determine A.

Thus, instead of explicitly calculating the time-dependent wavefunction, one can more simply and more insightfully study the evolution of electronic character by finding avoided crossings in the time-dependent electronic spectrum, extracting A and $ϵ$ from the fit to Equation (11) and putting these parameters into Equation (29).

The procedure we have just outlined still requires performing numerical calculations to obtain a time-dependent electronic spectrum. Further analytical insights at much less effort might have been obtained using the theoretical results from Section 2.3. As a result of the randomness of thermal fluctuations, one can model the electronic Hamiltonian by a random matrix. This theory allows determining the distribution of $ϵ,A$, and B for avoided crossings, Equation (13). In tandem with Landau–Zener theory, it is possible to understand the diabatic character of time evolution even without running simulations. In order to make random matrix theory applicable, it is necessary to adjust the density of states and the sensitivity to perturbations by rescaling energy and time; this procedure has been described by Wilkinson [22].

There is one caveat: Landau–Zener theory considers only two coupled states, whereas all states that form a cluster turn out to be coupled. There is a close and suggestive analogy with the Wigner surmise for the nearest neighbor level spacing distribution; a very accurate result was obtained by Wigner using only an ensemble of 2 × 2 random matrices instead of the full N × N case.

One can see a manifestation of the full N-state cluster coupling in Figure 3a for the first avoided crossing of the 18th state; the state makes a transition to the 21st but not to the adjacent 19th or 17th state. A general solution to the multi-level Landau–Zener problem is necessary here, which has not been obtained yet [44,45,46].

We have suggested a new way to analyze electronic dynamics in the presence of phonons, which enables to use random matrix theory and may shed new light on this problem.

5. Summary

The main purpose of this paper is to show that two seemingly distinct subjects, random matrix theory and electronic dynamics in the presence of thermal fluctuations, can be brought together through considering statistically known properties of the electronic structure spectrum. This perspective can potentially give more information than just the preservation of diabatic electronic character recently found in the time-dependent simulations of graphene thermal fluctuations.

We first introduced parameter-dependent random matrices and the reasons for selecting a particular parametrization. Then, we considered singularities in their spectra, i.e., conical intersections and avoided crossings. We have described the geometrical properties of conical intersections, along with a typical workflow for calculating their statistical distributions analytically. Qualitatively, analytical predictions were in agreement with the illustrated numerical results. Then, we quoted the statistical distribution of fitting parameters for energy levels near an avoided crossing, which is important for the perspective we suggest.

We next examined the time evolution of graphene $\pi$-band electronic states in the presence of phonons. A direct solution of time-dependent Schrodinger equation shows that a projection of the time-dependent state on an instantaneous adiabatic state drops to zero during the very first vibrational cycle, whereas overlap with the initial electronic state decays at a much longer timescale. Moreover, it tends to recover at an even larger timescale. That indicates preservation of the electronic character and the diabatic nature of time evolution.

Finally, we sketched a possible way to treat the dynamics of electrons in the presence of phonons based on random matrix theory and avoided crossings. We started with the observation that numerous avoided crossings arise in the supercell electronic bandstructure because of nuclear motion. Overlap of the time-dependent wavefunction with the instantaneous adiabatic states suggests avoided crossings are responsible for the rapid decay of the adiabatic states. A simple estimation based on the uncertainty principle confirms that avoided crossings have gaps that are too small for the adiabatic theorem to be valid. We then outlined a more sophisticated analytical treatment based on Landau–Zener theory. All in all, the suggested perspective may enable the use of random matrix theory to characterize the time evolution of electrons in the presence of phonons. This could be a powerful new perspective in this important problem.