How electronic dynamics with Pauli exclusion produces Fermi-Dirac statistics

A perfect theory of electronic dynamics must relax towards the correct equilibrium distribution of population at long times. However, satisfying this condition is not trivial. The two most common approaches for simulating non-equilibrium electronic dynamics in chemistry, Ehrenfest method and the surface hopping approach, do not satisfy this condition exactly. In the case of Ehrenfest dynamics, the violation of detailed balance is severe¹ and the dynamics can only be taken seriously at short times, although corrections can be made.² Surface hopping reaches equilibrium properly under limiting conditions.³ Both these methods are usually formulated in a basis of distinguishable (Boltzmannian) electronic states. Because the basis must be limited in practice, predictions in the Boltzmannian picture depend on a choice of adiabatic or quasi-diabatic states. States with fractional orbital occupations are expensive and difficult to treat in a basis of collective electronic states. To address these problems, we are pursuing relaxation dynamics based on the reduced electron density matrices.

Fermi blocking of electronic kinetics has been known since Landau and we are not the first people to derive master equations with Fermi-blocking. Other authors have proposed ad hoc semi-classical correction factors to account for Pauli blocking in Ehrenfest dynamics,⁴ which developed out of related semi-classical corrections.^5,6 In the one-electron picture, non-equilibrium Green’s function theories^7–12 and Boltzmann transport equations^13,14 also satisfy detailed balance. Master equations which exchange both electrons and energy between a system and a bath of two-level fermions have been recently derived, and they also approach a Fermi-Dirac (FD) distribution.^15,16 For spectroscopy, the most important bath is undoubtedly the bosonic vibrations of the atoms¹⁷ which is usually treated as a harmonic bath. Following the master equation formalism to treat the harmonic bath, other authors have also derived kinetic equations with blocking factors.^18–22 Relative to these previous papers, the advancement of the present work is that no mean-field approximation or assumption of separability on the electronic state is made in our derivation. Extensions of our technique can be used to produce systematically improvable treatments of a correlated electronic system with dissipation and are compatible with correlated models of electronic dynamics.^23,24 This feature is critical to understanding dissipation in multiple excited states and for molecules with strongly correlated ground states.

Ordinary system-bath perturbation theory with a bath of bosons does not tend towards FD populations; the equations do not consider exchange statistics so they approach a canonical equilibrium instead. The goal of this paper is to show how the derivation of ordinary system-bath perturbation theory naturally leads to a FD equilibrium if expectation values are taken correctly with respect to a mixed-state, Fermionic vacuum. The resulting kinetic equations resemble Redfield theory²⁵ but also depend on the hole density. The key to this approach is that the Fermionic expectation value is taken using the extended normal ordering technique of Mukherjee and Kutzelnigg (MK)²⁶ which does not require a wavefunction and only depends on reduced density matrices (RDMs). This leads to a non-linear equation of motion for the RDM. The exact equation of motion (EOM) for many-fermion RDMs, the Bogoliubov-Born-Green-Kirkwood-Yvon hierarchy (BBGKY),²⁷ is also non-linear. Non-linearity is also a requirement to obey Fermionic detailed balance. We show how a Pauli-blocking effect slows down relaxation and dephasing of electrons. In the Appendix and the supplementary material,²⁸ we outline a systematic treatment of correlated dissipation. We also use our new expressions to show that the rate of relaxation for a particle-hole excitation is not constant because of the Pauli-blocking effect.

The application of powerful technique of Mukherjee and Kutzelnigg to finite-temperature quantum mechanics is natural, and they hinted at this application in their classic paper. Recently, Hirata and He invoked the technique to correct perturbation theory in the thermodynamic limit.²⁹ Unfortunately, the generalized normal ordering technique is challenging to derive, understand, and apply. We must refer the reader to the original papers for details. A related generalized Wick’s theorem³⁰ concept was recently developed for Green’s function-based theories.³¹ For linear response away from a determinant, the blocking effect is irrelevant³² and so the literature is basically correct. A paper using these new formulas in a useful all-electron dynamics code is an immediate follow up. Electron transport methods based on the density matrix will also benefit from the main results of this paper.^33–36 Electronic statistics can be easily enforced as long as kinetic rates depend on the density matrix.

In the non-interacting case, the population of an electronic state at long times should tend to the FD distribution

where β = (k_bT)⁻¹, ϵ_i is the energy of a one-electron state, and μ is the chemical potential satisfying $∑ i n ̄ i =N$⁠, the total number of electrons. The bar indicates that this is an equilibrium quantity, where the number of particles exchanged between any two states is equal so the corresponding populations do not change (the condition of detailed balance).

It is known that a kinetic equation can be easily modified by “Pauli-blocking” factors (1 − n_i), to reach detailed balance with the FD distribution. Consider a set of non-negative Markovian transition rates between states K_ij satisfying $∑ j K i j =1$ and canonical detailed balance e^{β(ϵ_j−ϵ_i)} = K_ij/K_ji. The blocking factors close any channel which would disobey the Pauli exclusion principle. Inserting FD $n ̄ i ′ s$ into a kinetic equation of this form,

one sees the number of particles exchanged between any two states is equal, meaning the FD distribution is stationary. We have introduced η_i, the diagonal of the one-hole density matrix. Equation (2) is the form a kinetic equation must take to satisfy Fermionic detailed balance. Exclusion factors are an old concept attributed to Landau,³⁷ as an ad hoc correction. They are also important for dephasing rates.³⁸ In the remainder of this paper, we derive an expression with perturbation theory which reduces to Eq. (2).

We begin our derivation with the following Hamiltonian:

where $B ˆ ν$ is a coupling-weighted sum of boson operators attached to ν, defined below. To show that FD distribution is the exact equilibrium, the electrons must be non-interacting, but the correction we derive applies to the interacting case as well. We use a notation for second quantization which follows Ref. 26, and summation over repeated indices is implied. In the open-systems approach, the Hamiltonian is put into this form after diagonalizing the electronic part from some localized basis $|i〉= ∑ j U ˆ i ν |ν〉$ with electronic coupling. In the original basis, which is the atomic basis in our work, the electronic part of the system-bath coupling is usually taken to be diagonal:

but leads to off-diagonal coupling: $V i , j ν = ∑ ν U i , ν U j , ν †$ in the system energy eigenbasis used to propagate the dynamics. We discuss how Eq. (3) is prepared atomistically in Sec. IV. The equilibrium boson expectation values taken in this work are denoted 〈〉_eq, otherwise 〈〉 denotes a Fermionic expectation value. No restriction is made on the nature of the Fermionic state which is used to perform the expectation value. In particular, the formulas are valid for a general mixed Fermionic state which is completely defined by its RDMs and no wavefunction is invoked or implied. We also define a projection operator $P$⁠, which replaces bosonic operators with their expectation values at equilibrium leaving the fermion operators untouched, e.g., $P a a a b b † b= a a a b 〈 b † b 〉 eq$⁠, and its complement operator $Q= I ˆ −P$⁠.

The one-electron RDM (1-RDM) is defined as $γ b a =〈 a a a b 〉=〈 a b a 〉$⁠. We assume that an initial 1-RDM we would like to propagate is known, perhaps on the basis of a stationary electronic structure calculation. In the notation of this work which follows that of Ref. 26, γ is not an operator. Rather it is the expectation value of a one-particle operator. 1-RDMs are not equivalent to the tight-binding (TB) density matrices assumed in a master equation. However, Fermionic operators do satisfy Heisenberg’s EOM like TB density matrices. The steps used to produce master equations for TB operators generically apply to one-body Fermionic operators as well. Thus, until Fermionic expectation values are taken, our derivation is essentially the same as the textbook derivations of Redfield theory that can be found in Breuer³⁹ or Nitzan.⁴⁰ Our final expressions also reduce to the commonly used second-order TB theory in the limit that Fermi-blocking factors are removed.

Working in the interaction picture with respect to the first two terms of $H ˆ$⁠, where $L ( t ′ ) =−i [ V ˆ ( t ′ ) , ⋅ ]$⁠, second-order time convolutionless perturbation theory (TCL)^39,41,42 produces this EOM for the tensor product of $a b a$ with an equilibrium bosonic state

Taking the Fermionic expectation value of Eq. (5) will eventually produce an EOM for γ because $∂ ∂ t γ b a = ∂ ∂ t 〈 a b a 〉$⁠, but the equilibrium boson expectation values in $P$ must also be taken.⁴³ The first term on the right hand side of Eq. (5) is zero because it is sandwiched by $P$⁠, and the expectation value of odd-numbered boson operators is zero at equilibrium. The last term is an inhomogeneity which accounts for any initial correlations between electrons and bosons. Our previous work explored this term in detail.⁴¹ The effects on observables in that work were shown to be small relative to the Markov approximation, and other papers documented how inhomogeneous effects decay.^39,44 Because we are interested in the long-time detailed balance limit, we are well motivated to neglect the inhomogeneity.⁴⁵ This leaves only the second term, reflecting the time-dependent correlation between the system and the bath, which we denote as $d d t ( a b a ) ( 2 )$⁠,

At this point, we move back to the Schrödinger picture. We collect the time-dependent factors in each term (e^iω_cdt′ and e^iω_eft, where ħω_ab = ϵ_a − ϵ_b) and interaction strength scalars $V a , b ν$ together with the boson expectation value to focus on the Fermionic expectation value. An EOM for γ is obtained by taking the Fermionic expectation value of both sides of this equation. On the left hand side, $d d t e i ω a b t a b a =i ω a b e i ω a b t a b a + e i ω a b t d d t a b a$⁠. Dividing the common factor and taking the expectation value, $i ω a b 〈 a b a 〉=i ω a b γ b a$ produces the first term in Eq. (10).

Ordinarily, in order to take an expectation value of products of operators like on the right hand side of Eq. (6), one needs a description of an electronic state as excitations relative to a determinant wavefunction. Instead, we use MK’s extended normal ordering which advantageously depends only on RDMs. This double commutator also occurs as Eq. (48) in MK’s paper.²⁶ To obtain our expression from theirs, we first delete terms which contain normal-ordered operators (whose expectation values are therefore zero) and also make a simple substitution (⁠$δ b a = η b a + γ b a$⁠, with η the hole density matrix and γ the particle density matrix) so that $δ b a γ d c − δ d c γ b a = η b a γ d c − η d c γ b a$⁠, leaving four terms:

This expression has two interesting consequences. First, it shows that the expectation value of this double commutator does not depend on the higher-order RDMs. Even if the reference had strong, two-particle correlations or it was a double-excitation that does not affect this term.⁴⁶ This result supports the possibility of an accurate dissipative time-dependent Hartree-Fock equation, since the dissipation of the one-body density is not coupled to the two-body density through the bath.⁴⁷ Second, each factor in this EOM has a naturally derived Pauli-blocking factor.

The remaining steps in our derivation are equivalent to the ordinary derivation of a Redfield equation. If η is replaced with the identity matrix (blocking is shut off), the result is equivalent to an ordinary Redfield equation, and so readers can consult with one of the excellent textbooks.⁴⁰ The first assumption is a Markov approximation (⁠$lim t → ∞$⁠), in which the bath remains at boson equilibrium. This places a constraint on the boson correlation function,

which leads to the canonical detailed-balance relation in the Redfield equation, but for this equation furnishes canonical K’s like Eq. (2). Putting the result of the time integration together with the Fermionic expectation value for each term, we obtain

We note an interesting duality between the electron and the hole in this expression. On the right hand side, the third and fourth terms are particle relaxing but their blocking factors resemble hole-dephasing terms (they couple to only one index of η), whereas the first and second terms are particle dephasing but resemble hole-relaxing terms.

A secular approximation must be made to reach a kinetic equation which preserves trace, by deleting contributions which couple the diagonal and off-diagonal elements of γ. Besides the desired trace preservation, the secular approximation can be motivated with a rotating wave argument, which is identical to the Redfield case described in several previous publications.^48,49 Unit-less Kronecker deltas enforce this decoupling S_a,b,c,d = δ_abδ_cd + δ_bcδ_ad(1 − δ_abδ_cd). The resulting electronic EOM is the key result of our paper,

which is similar to the form of a secular Redfield equation but with blocking factors. We propagate this non-linear equation using a straightforward Runge-Kutta-Fehlberg adaptive integrator.⁵⁰ The hole density matrix η is not propagated separately; it is simply set equal to its definition in terms of γ. The kinetic equation rigorously preserves the trace of the 1-RDM and equilibrates to FD distribution. In the case that the density starts without coherence, η is diagonal and the equation takes the form of Eq. (2). We will now discuss the importance of these blocking factors in the evolution of populations and coherences using calculations on model systems.

First, we will demonstrate that propagation of Eq. (10) reproduces FD statistics. To produce a model with the form of Eq. (3), we consider a Hückel model of atomic sites. Each site has an energy E_ν = 1.0 a.u. and is coupled to its nearest neighbors with a coupling strength $V ˆ ν , ν e$⁠. The energy of each site fluctuates due to the system-bath coupling and the spectrum of these fluctuations is determined by a spectral density for this site $J ν ( ω ) = ∑ α ħ 2 ω α 2 g ν , α 2 δ ( ω − ω α )$ of the Drude-Lorentz form. We have previously shown a method to obtain these spectral densities atomistically from limited amounts of molecular dynamics, but for the purposes of this paper, we assume a strong, high-frequency bath⁵¹ to test positivity and rapidly achieve relaxation. After assigning a spectral density and coupling in the real-space ν basis, the electronic Hamiltonian is diagonalized to recover a Hamiltonian of form (3). Some of the dynamics are performed at a very high temperature so that the exact reproduction of FD-statistics is made clear. However, blocking has important effects even at room temperature, which we demonstrate with dynamics results at 300 K. In all cases, we half-fill the chain with electrons, and the chemical potential of the system is constant per Eq. (3).

The electronic dynamics resulting from Eq. (10) preserve the trace of the electronic density matrix and produce a FD distribution at long times. Simulations with several different coherent and incoherent initial conditions and a very strong bath all equilibrate to a FD distribution (Fig. 1). Ordinarily, solutions of the secular Redfield equation are exponential with time constants determined by the eigenvalues of K_ij, and hence, transient absorption spectra are fit to a superposition of decaying exponentials. However, our trajectories do not take the same form. Because of Fermi-blocking, particle-hole excitations and especially multiple particle-hole excitations do not relax with simple exponential kinetics. In Figure 2, we show the result of the simulation of an 8-site chain at T = 300 K, with otherwise unchanged parameters. We define a ground-state probability based on the overlap of the T = 0 K ground state and the density during the dynamics: $P g ( t ) = ∏ i ( 1 − | n i ( t ) − Θ ( ϵ i − μ ) | )$⁠, where Θ is the step function. We initialize the density in a particle-hole excited state by unfilling HOMO−2, filling HOMO+2, and monitor P_g(t) to observe the time dependence of the relaxation rates. For single particle-hole states, the rate of relaxation begins rapid as electrons fill in holes and then adopts a slower “blocked” and essentially exponential relaxation rate dominated by LUMO → HOMO relaxation at long times. They have direct implications for the kinetics of excited state relaxation as well as the interpretation of transient absorption spectra. However, the effective rate only varies by roughly a factor of ten because of the blocking effect, so it would be difficult to separate from other possible rate fluctuations in ultrafast spectroscopic experiments. In intrinsic semiconductors, this “blocked” effect is indirectly observed through the broadening of the apparent bandgap and is known as the Burstein-Moss shift.^52,53

Off-diagonal elements of a density matrix are called coherences and their decay is called dephasing. As has been recently reported by other authors, the dephasing rate of fermions is considerably different from what ordinary Redfield theory would predict.³⁸ In distinguishable models of electronic relaxation, coherence tends to be an ambiguous quantity because it is dependent on the chosen basis. In both the energy eigenbasis and the position eigenbasis, coherence has important and measurable implications. In the former, coherence should be zero in the thermodynamic limit. When they are nonzero, the energy-basis coherences oscillate. Our kinetic model does completely decay these coherences. If we initialize electrons on the left of our chain to reach detailed balance, they naturally flow to the right (Fig. 3). On the other hand, if only populations are allowed to relax, the populations of energy eigenstates rapidly reach a perfect FD distribution, but the undecayed coherences lead to unphysical oscillatory currents in the system. We emphasize that it is crucial to have an atomistic and improvable model for coherence decay.

In the position eigenbasis, coherences determine the measurable quantum delocalization of particles. Real-space coherence also limits the efficiency of electronic structure theory, since the cost of calculating exchange energy is completely determined by how slowly the real-space density matrix decays away from the diagonal. With our new theory, we can predict how dynamics and finite temperature effects can lead to the long-range decay of the off-diagonal elements of the electronic density matrix (Fig. 4). This locality can potentially be exploited to make electronic structure models more efficient.

Although kinetic equation (2) is non-linear in the density matrix, in our experience, it faithfully maps any density matrix to FD equilibrium at long times as it should, since FD density is a known solution of Eq. (10). We have not proven that FD is the only possible solution, but we strongly suspect it is based on repeated simulations. Fermionic relaxation rates must depend non-linearly on the density matrix to be correct at long times. Consequently, relaxation rates for Fermionic states are not constant as a function of time and tend to slow down over the course of non-radiative relaxation. Hole blocking makes rates of electronic relaxation and dephasing slower than the corresponding rates for distinguishable particles. In a Boltzmannian “exciton” model, this effect is extremely unclear.

These are second-order equations; therefore, the rates of transition they predict will become inaccurate in the strong coupling limit. The main feature of these equations is that they are systematically improvable. One can easily “Fermi-block” surface hopping rates and perform surface hopping for individual electrons. The results of this paper suggest that with Fermi blocking surface hopping for electrons would, for all intents and purposes, tend to FD. However, that electronic surface hopping would beg for a systematically improvable path to dephasing rates, like its Boltzmann counterpart.

We have found that the MK normal ordering is an especially useful technique for dynamics, allowing us to work around the impurity of electronic states at finite temperature and away from equilibrium. We used it to provide a promising EOM for electrons and simulate non-equilibrium vibronic relaxation. As in our previous work, the perturbation theory could also be used to develop electronic equations of motion that treat electron correlation effects (see the Appendix). This paper opens the door to answering whether interacting equilibrium differs meaningfully from the FD equilibrium in real materials, and how that equilibrium is approached. We are currently extending this work towards higher electronic density matrices.

It is significant that these equations approach Fermionic detailed balance naturally, as opposed to enforcing detailed balance in an ad hoc way. In a related vein, sometimes Pauli-blocking factors are sometimes fixed at their equilibrium value $η ̄ i$²¹ instead of η_i(t), which is correct at equilibrium. Other papers¹⁸ and this work derive blocking factors which are time-dependent and change with the density matrix. This leads to the interesting conclusion that even in a Markovian limit, the rate of relaxation is not constant. An interesting experimental/theoretical challenge would be to separate this effect from other non-Markovian rate changes in a real relaxation process. Applications of this work to predict rates of non-radiative relaxation in molecules are an active pursuit in our laboratory. We are also actively extending these equations to correlated electronic dynamics.

We thank the University of Notre Dame’s College of Science and Department of Chemistry and Biochemistry for generous start-up funding and Honeywell Corporation. T.S.N. gratefully acknowledges the support of the NSF Graduate Research Fellowship under Grant No. 1313583. We acknowledge valuable conversations with Liviu Nicolaescu (Notre Dame), Geoff Thompson (Iowa State), Jarrod McClean, Samuel Blau, and Thomas Markovich (Harvard).

The derivation technique employed in this work is more laborious than other approaches to electronic master equations^12,19,38 which simply invoke a Hartree-Fock Wick’s theorem. These papers obtain the correct blocking factor for products of single-particle operators, but they are limited to mean-field dynamics because they are based on an uncorrelated model of the electronic state. This limitation can be bypassed by virtue of the MK technique. To clarify the relationship with previous approaches, this appendix outlines how two beyond-mean field phenomena (electron-electron correlations and dissipation in correlated states) can be treated with extensions of the formalism developed in this paper.

To exactly calculate the energy of an electronic system, one must propagate a two-electron RDM (2-RDM):²³ $Γ r s p q =〈 a r s p q 〉$⁠. The 2-RDM has a separable component that can be captured by a determinant wavefunction $γ r p γ s q − γ r q γ s p$ and a non-separable cumulant part which is important for the correlation energy of the system $λ r s p q = Γ r s p q − ( γ r p γ s q − γ r q γ s p )$⁠. Following the derivation that leads to Eq. (7), a bath which couples to the 2-RDM will contribute a term to the EOM for the 2-RDM depending on the expectation value: $〈 [ [ a r s p q , a u t ] , a w v ] 〉$⁠. It was not previously possible to calculate this expectation value for a correlated, mixed electronic state in a systematic way without invoking some explicit wavefunction. However, this is a straightforward extension of the theory in this paper. We can show explicitly that the mean-field result would neglect contributions to $d d t 〈 a r s p q 〉$ from the cumulant of the 2-RDM. Using our automated version of the MK normal ordering written in Mathematica, we find the expectation value of the double-commutator involves not only mean-field like terms but also non-separable contributions,²⁸ 

where $A$ is an operator which antisymmetrizes upper and lower indices. Interestingly, the same decoupling of higher-order density matrices holds for 2-RDM. There are no terms which couple the 3-RDM cumulant to the 2-RDM via a single-particle bath.

We have previously shown that TCL perturbation theory makes a good second-order approach to incorporating electron-correlation effects in the propagation of the 1-RDM.⁴¹ When TCL theory is used in this way, the electron-electron interaction $V ˆ ee = V p q r s ⋅ a r s p q$ just replaces the system-bath coupling in the EOM, contributing a term to the EOM for $γ q p$ of the form,

Following this formalism, we developed corrections to configuration-interaction singles dynamics which resembled other second-order excited state theories and gave similar results for linear response spectra. We employed the Hartree-Fock normal ordering, making the assumption that the expectation value is taken with a determinant wavefunction: 〈〉 ≈ 〈Ψ_SD||Ψ_SD〉. This approximation limits the accuracy of dynamics, and has plagued EOM models for a long time. Some recent theories use correlated vacuum states^31,54,55 to avoid this approximation. To handle the correlated expectation value, these methods make factorization approximations⁵⁴ or use a complete basis for the state⁵⁵ which is intractable. The techniques developed in this paper allow us to relax this approximation in a way which is totally systematic and can be automated.