Time-dependent Orbital-free Density Functional Theory: Background and Pauli kernel approximations
TTime-dependent Orbital-free Density Functional Theory: Background andPauli kernel approximations
Kaili Jiang ∗ Department of Chemistry, 73 Warren St.,Rutgers University, Newark, NJ 07102, USA
Michele Pavanello † Department of Chemistry, 73 Warren St.,Rutgers University, Newark, NJ 07102, USA andDepartment of Physics, 101 Warren St., Rutgers University, Newark, NJ 07102, USA (Dated: February 12, 2021)Time-dependent orbital-free DFT is an efficient method for calculating the dynamic prop-erties of large scale quantum systems due to the low computational cost compared to standardtime-dependent DFT. We formalize this method by mapping the real system of interactingfermions onto a fictitious system of non-interacting bosons. The dynamic Pauli potentialand associated kernel emerge as key ingredients of time-tependent orbital-free DFT. Usingthe uniform electron gas as a model system, we derive an approximate frequency-dependentPauli kernel. Pilot calculations suggest that space nonlocality is a key feature for this ker-nel. Nonlocal terms arise already in the second order expansion with respect to unitlessfrequency and reciprocal space variable ( ωq k F and q k F , respectively). Given the encourag-ing performance of the proposed kernel, we expect it will lead to more accurate orbital-freeDFT simulations of nanoscale systems out of equilibrium. Additionally, the proposed pathto formulate nonadiabatic Pauli kernels presents several avenues for further improvementswhich can be exploited in future works to improve the results. ∗ [email protected] † [email protected] a r X i v : . [ phy s i c s . c h e m - ph ] F e b I. INTRODUCTION
The simulation of large scale quantum systems (such as nanometer-scale quantum dots) whentheir electrons are out of equilibrium has been a challenging undertaking and cause for frustration[1–5]. The challenge is multifaceted. First, the simulation methods need to be predictive, andalthough time-dependent DFT (TD-DFT) is the workhorse for these types of simulations, it stilllacks broad applicability across the possible excited states characters (valence, charge transfer,Rydberg) [6–13]. Second, the computational cost is a major concern. The cubic scaling of theground state DFT algorithm and also of the real-time or linear-response TD-DFT algorithms(provided a small number of states are computed) cripple the applicability of these methods tonanoscale systems [14, 15]. The problem is further exacerbated when systems with dense spectraare considered [16–18]. Beyond DFT, an array of accurate methods is available. However, theircomputational cost is typically orders of magnitude larger [19].In this work, we tackle issues of computational feasibility by exploring an alternative path:orbital-free DFT (OF-DFT) [20, 21]. OF-DFT effectively reduces the complexity of the problemby considering only a single active orbital. This is in stark contrast with commonly adopted Kohn-Sham DFT (KS-DFT) methods in which a set of occupied orbitals equal in number to the electronsin the system needs to be considered [22, 23]. Thus, OF-DFT massively reduces the complexityof the problem, provided that accurate approximations for the density functionals involved areavailable. An important feature distinguishing OF-DFT from KS-DFT is that in addition to theexchange-correlation (XC) functional, OF-DFT also requires the knowledge of the noninteractingkinetic energy functional.Employing currently available density functional approximants, the ground state version of OF-DFT scales linearly with the system size and has been shown to be applicable to main group metalsand III-V semi-conductors [24–28]. Recently, non-local functionals such as HC [25] and LMGP [29]have achieved chemical accuracy for a wider range of systems.OF-DFT can be formulated in the time domain to approach systems out of equilibrium [30–42]. It is sometimes referred to as the time-dependent Thomas-Fermi method [31–33] and alsohydrodynamic DFT [36–39]. In this work, we will refer to it as TD-OF-DFT. Practical imple-mentations initially featured the adiabatic Thomas-Fermi (TF) [43, 44] plus von Weizs¨acker (vW)[45] approximation [30, 32] (TFW, hereafter), and later including nonadiabatic corrections in thepotential [40–42]. TD-OF-DFT has seen a wide range of applications including atoms and clustersin laser fields, electron dynamics and optical response in nano-structures, and oscillators in electricfields [32, 33, 46–50]. Despite the drastic approximations, TD-OF-DFT has been quite successfulin describing the optical spectra of metal clusters [32, 51].In this work we formulate exact conditions which are a foundational aspect of functional de-velopment [52] and focus particularly on the TD-OF-DFT Pauli kernel, its nonadiabatic behavior,spatial nonlocality and approximations needed for practical calculations.To cast our developments in the current state of the art, we should mention other quantumdynamic methods, such as time-dependent density functional tight binding (TD-DFTB) [53–56],which is capable of handling systems of much greater size than conventional TD-DFT as well assimplified versions of TD-DFT [57–60].This work is organized as follows. In section II, we first establish an exact map between the realsystem of interacting electrons, a system of noninteracting electrons and a system of noninteractingbosons. This allows us to formulate TD-OF-DFT as an exact formalism with theorems and exactproperties/conditions for the building blocks of the method. In section III we derive the real-time formalism of this approach (i.e., the Schr¨odinger-like equation), and the Pauli potential thatneeds to be approximated in actual calculations, as well as some properties of the exact Paulipotential. In section IV we focus on the linear response formalism and derive Dyson equationsrelating the response to external perturbations of the real system, the KS system and the fictitiousnoninteracting boson system. We then proceed to uncover a route to approximate the Pauli kernelrecovering needed nonadiabaticity and nonlocality. Pilot calculations are also presented to showcasethe newly developed Pauli kernel.
II. NONINTERACTING BOSON SYSTEM
The foundation of ground state Kohn-Sham (KS)-DFT and TD-DFT is the existence of aunique and invertible map between the real system of interacting fermions and a fictitious systemof noninteracting fermions (KS system, hereafter). However, other unique and invertible maps canbe found between the real system and other fictitious systems. For example, it is possible to usea system of bosons having the same charge, point-like shape and mass as the electrons. Such abosonic system is fictitious, and thus does not need to be rooted in reality. A bosonic system atzero temperature yields the following wavefunction to density relationship, n ( r ) = N | φ B ( r ) | . (1)OF-DFT implicitly takes advantage of the fictitious bosonic system because it can be formulatedin a way that reduces to Eq.(1), i.e., only utilizes a single active orbital. This can be seen byinvoking the KS-DFT total energy functional E [ n ] = T S [ n ] + E H [ n ] + E XC [ n ] + (cid:90) d r n ( r ) v ( r ) , (2)where T S is the kinetic energy of N noninteracting electrons having a density of n ( r ). E H is theHartree energy, E XC is the XC energy, and v ( r ) is the external potential. A search over electrondensities to find the minimum of the energy functional must be carried out with the constraint thatthe densities must integrate to N electrons. Thus, the following Lagrangian is typically invoked, L [ n ] = E [ n ] − µ (cid:20) (cid:90) d r n ( r ) − N (cid:21) . (3)Formally, the minimum of the above Lagrangian is also a stationary point [61] and therefore itsfunctional derivative with respect to the electron density must vanish,0 = δT S [ n ] δn ( r ) + v S ( r ) − µ, (4)which is known as the Euler equation of DFT. If the first term on the RHS of the above equation iswritten as the vW potential, v vW ( r ) = √ n ( r ) (cid:16) − ∇ (cid:112) n ( r ) (cid:17) , plus a correction, δT S [ n ] δn ( r ) = v vW ( r ) + (cid:16) δT S [ n ] δn ( r ) − v vW ( r ) (cid:17) , then because the vW potential is exact for wavefunctions of up to two electronsand for bosonic systems, the correction term is known as Pauli potential, v P ( r ), as it accounts forFermi statistics.The single-orbital wavefunction now emerges as simply φ B ( r ) = √ N (cid:112) n ( r ) and the followingSchr¨odinger-like equation derives directly from Eq.(3), (cid:18) − ∇ + v P ( r ) + v S ( r ) (cid:124) (cid:123)(cid:122) (cid:125) v B ( r ) (cid:19) φ B ( r ) = µφ B ( r ) . (5)Following Eq.(1), the many-body wavefunction of the noninteracting boson system is a Hartreeproduct: Φ B ( r , r , . . . , r N , t ) = N (cid:89) l φ B ( r l , t ) . (6)In this article we use the subscript B for the noninteracting boson system, subscript S for the KSsystem, and use no subscript for the real system.When considering systems away from equilibrium (i.e., when the density, wavefunction andpotentials become time-dependent) the DFT theorems remain largely valid [62]. Before going in tothe details of the formal proofs, we summarize in Table I the most important quantities involved TABLE I. Comparison between the interacting, KS and noninteracting boson system.Interactingsystem KS system Noninteractingboson systemdensity n ( r , t ) n ( r , t ) n ( r , t )wavefunction Ψ( r , . . . , r N , t ) √ N ! det[ { φ S ,l ( r l , t ) } ] (cid:81) Nl φ B ( r l , t )effectivepotential v ( r , t ) v S [ n ]( r , t ) = v ( r , t ) + v H [ n ]( r , t )+ v XC [ n ]( r , t ) v B [ n ]( r , t ) = v S [ n ]( r , t ) + v P [ n ]( r , t )Hamiltonian ˆ T + ˆ W + v ( r , t ) ˆ T + v S ( r , t ) ˆ T + v B ( r , t ) in the three reference systems (interacting, KS and noninteracting boson). The wavefunctionsof the three systems are different. For the interacting system it can be a very complex functionof coordinates of every electron and time; for the KS system it is a single Slater determinant;and, as mentioned, for the noninteracting boson system it is a simple Hartree product of thesame function, φ B . The interacting system features an electron-electron Coulomb repulsion termˆ W in the Hamiltonian, whereas the two noninteracting systems do not. To compensate for themissing electron-electron interaction, the KS system employs an effective potential, v S [ n ]( r , t ),which includes the external potential v ( r , t ), the Hartree potential, v H [ n ]( r , t ), and the XC potential, v XC [ n ]( r , t ). For the noninteracting boson system, there is a different effective potential, v B [ n ]( r , t ),which includes an additional Pauli potential term, v P [ n ]( r , t ), required to compensate for the neglectof the Pauli exclusion principle. In the end, all three systems yield the same time-dependent density.To provide a visual and formal framework, we introduce bijective maps connecting the threesystems and realized in practice by the v S and v P potentials. We shall indicate with Γ, Γ S and Γ B bijective maps linking every v -representable time-dependent density to a unique time-dependentmulti-electron wavefunction Ψ( r , . . . , r N , t ), every non-interacting v S -representable densities to aunique set of KS orbitals { φ l ( r l , t ) } , every non-interacting v B -representable densities to a uniquesingle orbital φ B ( r , t ), respectively (see Figure 1).To establish the exactness of these maps, we now prove analogs of two fundamental theoremsin TD-DFT. The first theorem establishes one-to-one bijective maps between time-dependent den-sities, time-dependent wavefunctions and time-dependent effective potentials for any one of thesystems, particularly for the noninteracting boson system. From Eq. (1) and Eq. (6), given a wave-function the density is uniquely determined, and with a given density, the wavefunction is uniquely FIG. 1. Three representations of an electronic system discussed in this work. The circles represent electronsand the dotted lines represent electron-electron interactions. The noninteracting fermion and boson systemsare depicted in the bottom left and right enclosures, respectively, while the interacting (real) system is atthe top. Grey circles represent bosons. All three systems have a one to one bijective map between itswavefunction, effective potential and time-dependent density, n ( r , t ) which are indicated by the symbols Γ,Γ S and Γ B . The density is the same across all representations. determined up to a phase factor. The forward one-to-one map between the time-dependent effectivepotentials and time-dependent densities is also found by pluging the potential in the correspondingHamiltonian and solving for the time-dependent Schr¨odinger-like equation. The reverse one-to-onemap between the time-dependent effective potentials and time-dependent densities is proved bythe Runge-Gross theorem [63]. The proof of this theorem does not require the Hamiltonian andthe wavefunction of the system to be in a particular form, except for the effective potential tobe Taylor expandable [64]. Therefore, it applies to the noninteracting boson system without anyeffort.Next, we are going to establish a one-to-one map between the external potential of the inter-acting system and the effective potential of the noninteracting boson system (up to a constant).This requires to prove an analog of the van Leeuwen theorem [65]: For any time-dependent density, n ( r , t ), associated with the system of interacting fermions and the external potential, v ( r , t ), andinitial interacting state, Ψ , there exist a unique potential v B ( r , t ) up to a purely time-dependentconstant, c ( t ), that will reproduce the same time-dependent density of a system of noninteractingbosons, where the initial state Φ B of this system must be chosen such that the densities and theirtime derivative of the two systems are the same at the initial time.The proof of our analog of the van Leeuwen theorem is very similar to that of the originalone. Here we are going to follow Ref. 62 and only point out the differences. We start with theHamiltonian of the two systems, the interacting fermions,ˆ H = ˆ T + ˆ W + ˆ V ( t ) , (7)with the time-dependent wavefunction Ψ( t ) and initial state Ψ , and the noninteracting bosons,ˆ H B = ˆ T + ˆ V B ( t ) , (8)with the time-dependent wavefunction Φ B ( t ) and initial state Φ B . Both systems yields the sametime-dependent density n ( r , t ). We can Taylor expand v ( r , t ) v ( r , t ) = ∞ (cid:88) k =0 k ! v k ( r )( t − t ) k , (9)where v k ( r ) = ∂ k v ( r , t ) /∂t k | t = t , and v B ( r , t ) v B ( r , t ) = ∞ (cid:88) k =0 k ! v k B ( r )( t − t ) k , (10)where v k B ( r ) = ∂ k v B ( r , t ) /∂t k | t = t . Our goal is to find a way to uniquely determine each Taylorexpansion coefficients v k B ( r ).Next, we note that the current density operator should be the same for both systems. Namely,its expectation value ˆ j ( r ) = 12 i N (cid:88) l [ ∇ l δ ( r − r l ) + δ ( r − r l ) ∇ l ] . (11)should be the same whether it is evaluated with the KS wavefunction or the bosonic wavefunction.Thus, applying Eq. (11) to the wavefunctions of the interacting, KS, and noninteracting bosonsystems we obtain the current densities for these systems, j ( r , t ) = (cid:68) Ψ( t ) (cid:12)(cid:12)(cid:12) ˆ j ( r ) (cid:12)(cid:12)(cid:12) Ψ( t ) (cid:69) , (12) j S ( r , t ) = 12 i N (cid:88) l [ φ S ∗ ,l ( r , t ) ∇ φ S ,l ( r , t ) + φ S ,l ( r , t ) ∇ φ S ∗ ,l ( r , t )] , (13) j B ( r , t ) = N i [ φ ∗ B ( r , t ) ∇ φ B ( r , t ) + φ B ( r , t ) ∇ φ ∗ B ( r , t )] . (14)If we apply the equation of motion for ˆ j ( r ) in the noninteracting boson system, we will get (SeeAppendix A for detailed derivation) ∂∂t j B ( r , t ) = − n ( r , t ) ∇ v B ( r , t ) − F kin ( r , t ) . (15)Where the kinetic force F kin is in the following form F kin ( r , t ) = N { ∇ φ ∗ B ( r , t )][ ∇ φ B ( r , t )] + 4[ ∇ φ ∗ B ( r , t )][ ∇ φ B ( r , t )] − ∇ [ φ ∗ B ( r , t ) φ B ( r , t )] } . (16)Eq.(15) is in the same form as Eq. (3.27) in Ref. 62 if we consider the interaction force F int = 0because the system is assumed to be noninteracting. Take the divergence of Eq.(15) and use thecontinuity equation, we obtain ∂ ∂t n ( r , t ) = ∇ · [ n ( r , t ) ∇ v B ( r , t )] + q B ( r , t ) , (17)where q B ( r , t ) = ∇ · F kin ( r , t ). Subtract Eq.(17) from its counter part for the interacting systemEq.(3.48) in Ref. 62, we obtain ∇ · [ n ( r , t ) ∇ γ ( r , t )] = ζ ( r , t ) , (18)where γ ( r , t ) = v ( r , t ) − v B ( r , t ) and ζ ( r , t ) = q B ( r , t ) − q ( r , t ). Note that Eq.(18) is in the exactlysame form as Eq.(3.50) in Ref. 62.The last part of the proof is to derive the equations that uniquely determine v k B ( r ) from Eq.(18).This part of the proof is exactly the same as those from Eq.(3.51) to Eq. (3.55) in Section 3.3 ofRef. 62. Now we reach the conclusion which is that v B ( r , t ) is uniquely determined by the densityand the initial states, e.g., v B ≡ v B [ n, Ψ , Φ B ]( r , t ) . (19)We formally proved that for every interacting system with a v -representable time-dependentdensity and given an initial condition, there is a unique, noninteracting boson system that yieldsthe same density. This is important because without such a map a reference noninteracting bosonsystem cannot be employed. We recall that this approach has been analyzed before by severalauthors in several regimes under different sets of approximations [1, 2, 32, 34, 36, 39, 40, 42, 66–69] as well as in the context of the exact factorization method of recent formulation [70]. In thenext section, we discuss how to determine the effective bosonic potential in practical calculationsand we also review some of its exact properties that can be useful for guiding the development ofapproximations. III. ADDITIONAL THEOREMS AND PROPERTIESA. Time-dependent Schr¨odinger-like equation
The typical setup is that the system starts in its ground state at t = t and begins to evolveunder the influence of a time-dependent external potential for t > t . Thus, the initial wavefunction, Φ B , can be found with ground-state OF-DFT by solving Eq.(5). In the following, it willbe convenient to define the time-independent Pauli potential as a functional derivative of the Paulikinetic energy. Namely, v P [ n ]( r ) = δT P [ n ] δn ( r ) , (20)where T P [ n ] = T S [ n ] − T vW S [ n ].Given the initial orbital φ B , φ B ( r ) = 1 √ N (cid:112) n ( r ) , (21)where n ( r ) is the ground-state density at t = t , in order to propagate the system after t > t , weneed to solve a time-dependent Schr¨odinger-like equation: (cid:20) − ∇ v B [ n ]( r , t ) (cid:21) φ B ( r , t ) = i ∂∂t φ B ( r , t ) , (22)with the initial condition φ B ( r , t ) = φ B ( r ) (23)and the time-dependent effective potential is given by v B [ n ]( r , t ) = v P [ n ]( r , t ) + v S [ n ]( r , t ) , (24)which defines the time-dependent Pauli potential. B. Properties of v P ( r , t ) While the exact form of the time-dependent Pauli potential is unknown for the general case, wecan still derive some properties and exact conditions. Eq.(24) indicates that the role of the time-dependent Pauli potential in TD-OF-DFT is similar to the role of time-dependent XC potential inTD-DFT because it ensures that the fictitious boson system has the same electronic dynamics asthe fermion system. In TD-DFT, the XC potential ensures that the dynamics of the noninteracting0fermion system is the same as the interacting fermion system. Therefore, for many properties ofthe time-dependent XC potential, we can find analogies for the time-dependent Pauli potential.Here we are going to list some of such properties.
1. Functional dependency
The functional dependencies of v P , v S and v B , in terms of n , Ψ , Φ S and Φ B v S ≡ v S [ n , Ψ , Φ S ]( r ) , (25) v B ≡ v B [ n , Ψ , Φ B ]( r ) , (26) v P ≡ v P [ n , Φ S , Φ B ]( r ) . (27)Eq.(27) links the KS system with the noninteracting boson system. Thus, it does not depend onthe interacting wavefunction.
2. The zero-force theorem
The total momentum of a many-body system can be expressed as [62] P ( t ) = (cid:90) d r j ( r , t ) . (28)Combine Eq.(15) and Eq.(28) and notice that F kin in the form of Eq.(16) is the divergence of astress tensor and thus the integral of F kin over the space vanishes, we obtain ∂ P B ( t ) ∂t = − (cid:90) d r n ( r , t ) ∇ v B ( r , t ) . (29)The total momentum of the noninteracting boson system and the KS system are the same atall time because the densities of the two systems are the same, thus0 = ∂ P S ( t ) ∂t − ∂ P B ( t ) ∂t = − (cid:90) d r n ( r , t ) ∇ [ v S ( r , t ) − v B ( r , t )] (30)or (cid:90) d r n ( r , t ) ∇ v P ( r , t ) = 0 (31)The implication of this theorem are that v P ( r , t ) cannot exert a total, net force on the electronicsystem. This is also the case for the XC potential [71].1
3. The one- or two-electron limit
Similar to the XC potentials in KS-DFT and TD-DFT, an exact condition for the Pauli potentialarises for one- and two-electron densities. In the event of the system having only one electron ortwo spin-compensated electrons, the KS system is given by a single orbital. This orbital is thesame as the orbital of the noninteracting boson system. Thus, v P ( r , t ) = 0 for (cid:90) d r n ( r ) ≤ . (32)Even though this is a clear exact condition, its imposition in real-life density functional approxi-mations is nearly impossible resembling the self-interaction error for the XC functional [72].
4. The relation between T P and v P Let us define the time-dependent energy of the noninteracting boson as E ( t ) = (cid:68) Φ B ( t ) (cid:12)(cid:12)(cid:12) ˆ H B ( t ) (cid:12)(cid:12)(cid:12) Φ B ( t ) (cid:69) . (33)We then plug Eq. (8) and Eq. (24) in Eq.(33) to obtain E ( t ) = T B ( t ) + T P ( t ) + E H ( t ) + E XC ( t ) + (cid:90) d r v ( r , t ) n ( r , t ) , (34)where T B ( t ) = − (cid:112) n ( r , t ) ∇ (cid:112) n ( r , t ) , (35) T P ( t ) = T S ( t ) − T B ( t ) , (36) E H ( t ) = 12 (cid:90) d r (cid:90) d r (cid:48) n ( r , t ) n ( r (cid:48) , t ) | r − r (cid:48) | , (37) E XC ( t ) = T ( t ) − T S ( t ) + W ( t ) − E H ( t ) . (38)Note that the definition of E H ( t ) and E XC ( t ) are exactly the same as those in TD-DFT [62].If we apply the Heisenberg equation of motiond O d t = (cid:42) ∂ ˆ O∂t (cid:43) + i (cid:104) [ ˆ H, ˆ O ] (cid:105) (39)to the Hamiltonians of the KS and the noninteracting boson system, we get [62]d T S ( t )d t = − (cid:90) d r ∂n ( r , t ) ∂t v S ( r , t ) , (40)2and (see Appendix B for detailed derivation)d T B ( t )d t = − (cid:90) d r ∂n ( r , t ) ∂t v B ( r , t ) . (41)Subtracting Eq.(41) from Eq.(40) we obtaind T P ( t )d t = (cid:90) d r ∂n ( r , t ) ∂t v P ( r , t ) . (42)Eq.(42) can be used to check the accuracy of numerical calculations, or as an exact constraintfor approximations of v P as done for other functionals [73]. This is also the case of its analog in v XC [74] being used in developing XC functionals [75].
5. Nonadiabaticity and causality
It is possible to define time-dependent potentials and kernels from action functionals. Here,we use an analogy of the variational principle method proposed by Vignale [5] which takes intoaccount the “causality paradox” [1].The action integral of the noninteracting boson system can be written as: A B [ n ] = A B [ n ] − (cid:90) t t d t (cid:90) d r n ( r , t ) v B ( r , t ) , (43)where A B [ n ] = (cid:90) t t d t (cid:42) Φ B [ n ]( t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) i ∂∂t − ˆ T (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Φ B [ n ]( t ) (cid:43) . (44)Imposing A B in Eq. (43) to be stationary with respect to variations of the density leads to v B ( r , t ) = δA B [ n ] δn ( r , t ) − i (cid:28) Φ B [ n ]( t ) (cid:12)(cid:12)(cid:12)(cid:12) δ Φ B [ n ]( t ) δn ( r , t ) (cid:29) . (45)Relating the action integrals of the noninteracting boson system and the KS system, we obtain A B [ n ] = A S [ n ] + A P [ n ] , (46)and carrying out a similar analysis as Eq.(43) through Eq. (45) for the KS system, plugging in theresults of both systems into Eq.(46), we obtain v P ( r , t ) = δA P [ n ] δn ( r , t ) + i (cid:28) Φ S [ n ]( t ) (cid:12)(cid:12)(cid:12)(cid:12) δ Φ S [ n ]( t ) δn ( r , t ) (cid:29) − i (cid:28) Φ B [ n ]( t ) (cid:12)(cid:12)(cid:12)(cid:12) δ Φ B [ n ]( t ) δn ( r , t ) (cid:29) . (47)The above equation highlights the fact that in order to avoid the so-called “causality paradox”,the functional derivative of A P should be augmented by two boundary terms, one stemming fromthe KS system and one from the noninteracting boson system. This is an analog to the functionalderivative of A XC in conventional TD-DFT augmented by the boundary terms from the interactingsystem and the KS system.3 IV. APPROXIMATING THE PAULI KERNEL
In this section, we will first derive appropriate Dyson equations relating response functions ofthe interacting system, the KS system and the noninteracting boson system. In a second step, weanalyze the poles of these response functions to better appreciate the nonadiabaticity of the Paulikernel. Finally, we derive an approximate Pauli kernel and test it in several pilot calculations.
A. Response functions and Dyson equations
Dyson equations relating the interacting, KS and the noninteracting boson systems are impor-tant as they give us a framework to formulate approximations for the involved kernels [76–79].In linear response theory, we expand the time-dependent density in orders of v B (where v B isthe first order potential term) and the first order change in the density, n , is given by n ( r , t ) = (cid:90) d t (cid:48) (cid:90) d r (cid:48) χ B ( r , t, r (cid:48) , t (cid:48) ) v B ( r (cid:48) , t (cid:48) ) , (48)where the density-density response function is χ B ( r , t, r (cid:48) , t (cid:48) ) = δn [ v B ]( r , t ) δv B ( r (cid:48) , t (cid:48) ) (cid:12)(cid:12)(cid:12)(cid:12) v B [ n ]( r ) . (49)Defining the time dependent Pauli kernel f P ( r , t, r (cid:48) , t (cid:48) ) = δv P [ n ]( r , t ) δn ( r (cid:48) , t (cid:48) ) (cid:12)(cid:12)(cid:12)(cid:12) n , (50)and by noticing that both χ B and f P only depend on t − t (cid:48) , they can be represented in frequencyspace using Fourier transformations χ B ( r , r (cid:48) , ω ) = (cid:90) d( t − t (cid:48) ) e iω ( t − t (cid:48) ) χ B ( r , r (cid:48) , t − t (cid:48) ) , (51)and same for f P . We can now introduce the relevant Dyson equations, f P ( r , r (cid:48) , ω ) = χ − B ( r , r (cid:48) , ω ) − χ − S ( r , r (cid:48) , ω ) , (52) f PHXC ( r , r (cid:48) , ω ) = χ − B ( r , r (cid:48) , ω ) − χ − ( r , r (cid:48) , ω ) , (53) f HXC ( r , r (cid:48) , ω ) = χ − S ( r , r (cid:48) , ω ) − χ − ( r , r (cid:48) , ω ) . (54)Eq.(52) shows the role the Pauli kernel f P plays in linking the poles of the response functionsof the noninteracting boson system and the KS system. Clearly, the response functions of the twosystems have poles at different frequency values. To further investigate this, in Figure 2 and Table4 D e n s i t y o f S t a t e s ( e V ) Ag KSBoson-TF0.0 0.5 1.0 1.5 2.0 2.5 3.0KS 0.0 0.5 1.0 1.5 2.0 2.5 3.0Boson-TF 0 5 10 15 20 25 30Energy (eV)Boson-Inversion
FIG. 2. Top panel: Density of states of a Ag rod (see depiction in inset of Figure 6). Red solid line:noninteracting boson system using the TF potential as v P . Blue dashed line: KS system. Both systemsare computed with the same XC functional (LDA). Bottom three panels show the energy value of the polesof various response functions. Blue dots: KS system; red dots: approximate noninteracting boson systemwith TF as v P ; black dots: noninteracting boson system with v P found by inversion of the KS density (seeEq.(55)). Note: the bottom panel x-axis runs from 0 to 30 eV as opposed to the other panels which runfrom 0 to 3 eV. II we compare the frequency positions of the poles of the noninteracting boson system and theKS system for a Ag nano-rod [16] computed by simply diagonalizing the KS and noninteractingboson Hamiltonians. LDA is used as the XC potential to determine v S , and TF or an exact inver-sion are used to determine v P . TF is equivalent to employing the Thomas-Fermi-von-Weizs¨ackerapproximation (TFW) for T S [ n ] while the exact inversion of the KS system is carried out via thefollowing relation v B [ n ]( r ) = ∇ (cid:112) n ( r )2 (cid:112) n ( r ) , (55)where n ( r ) is KS-DFT density.5 TABLE II. Comparison of the lowest bosonic and KS pole in eV for the Ag nano-rod and icosahedral Na cluster (see Figure 6 for a depiction of their structures and Figure 2 for additional information).System KS Boson-TF Boson-inversionAg Figure 2 and Table II shows that for the Ag system the onset of the poles of the KS responsecomes at significantly lower frequencies compared to the corresponding poles of the boson response.It is interesting to note that the approximate TF v P yields excitation energies for the boson systemthat are much closer to the KS excitations compared to the true boson system computed byinversion.This demonstrates the role of f P which is to red shift the poles of the boson response towardsthose of KS response in addition to changing their character by mixing the boson excited statesthrough the Pauli kernel. Unlike the poles of KS response which are generally close to the polesof the interacting response, the differences between poles of boson and KS response seem muchlarger, highlighting the importance of accounting for nonadiabaticity in f P . B. Approximating f P accounting for nonadiabaticity and nonlocality The KS response function for the noninteracting uniform electron gas (free electron gas, or FEG,hereafter), aka the frequency-dependent Lindhard function, plays an important role in TD-DFT forderiving approximations to the XC kernel [80–86] and is also used as the target response functionfor parametrizations of nonadiabatic Pauli potentials [40, 66]. Following a similar paradigm, wefirst derive the exact response function for the bosonic FEG (BFEG) and then we use the Dysonequation in Eq. (59) to derive approximations to the Pauli kernel.For deriving the BFEG response function, assume N noninteracting bosons confined in a d -dimensional cubic cell with the volume of L d , with a uniform external potential. Assuming spincompensation, the BFEG response function can be written as [87]: χ B ( q, ω ) = 1 L d (cid:88) k n k − n k + q ω + ε k − ε k + q + iη , (56)where n k = 1 / ( e ( ε k − µ ) /T −
1) is the Bose-Einstein average occupation of state k at temperature T and chemical potential µ . We refer the reader to chapter 4 of Ref. 87 for a thorough introductionand derivation of the Lindhard functions in 1, 2 and 3 dimensions.6Eq.(56) can be rewritten with a change of variable k → k − q : χ B ( g, ω ) = 1 L d (cid:32)(cid:88) k n k ω + ε k − ε k + q + iη + (cid:88) k n k − ω + ε k − ε k − q − iη (cid:33) , (57)where the summation is over all occupied states. In the limit of T →
0, only k = 0 statesare occupied with the occupation number N , therefore we can replace the summation with amultiplication of N : χ B ( q, ω ) = NL d (cid:18) ω − q / iη + 1 − ω − q / − iη (cid:19) = k F π (cid:18) ω − q / iη + 1 − ω − q / − iη (cid:19) , (58)where the Fermi wavevector, k F = (3 π n ) / . R e () - I m () SBTFW
FIG. 3. Comparison between the Lindhard function ( χ S for the FEG), χ B for the BFEG (see Eq.(58)) andthe Thomas-Fermi-von-Weizs¨acker (TFW) response function χ TFW (Eq. 18 of Ref. 66) with respect to q fora uniform density value of n = 0 .
004 a.u. (corresponding to a Fermi wavevector value of k F = 0 .
491 a.u.),the excitation frequency is set to ω = 0 .
125 a.u. and spectral broadening is set to η = 0 . We note that the formula for χ B of the BFEG is simpler than the one for χ S of the FEG dueto the omission of the integral over k . Figure 3 shows the comparison between χ S , χ B and theapproximate response function for TFW functional χ TFW for the FEG. χ B has the same asymptoticbehavior as χ S for both q → q → ∞ . It has a similar shape to χ TFW but the singularityis at a different position. Additionally, both χ B and χ TFW feature vary narrow peak dispersionscompared to χ S .Combining Eq.(58), the Lindhard function and the Dyson equation Eq.(52), we can generatethe exact Pauli kernel for the FEG (in reciprocal space):7 f FEG P ( q, ω ) = 3 π k F (cid:18) ω − q / iη + 1 − ω − q / − iη (cid:19) − − π qk F (cid:20) Ψ (cid:18) ω + iηqk F − q k F (cid:19) − Ψ (cid:18) ω + iηqk F + q k F (cid:19)(cid:21) − , (59)where Ψ ( z ) = z − z z + 1 z − . (60) R e ( f ) - I m ( f ) = 0= 0.125= 0.25 FIG. 4. Pauli kernel, f P , for the free electron gas. Plots of f P with respect to q are shown for different ω (ina.u.) at n = 0 .
004 a.u. ( k F = 0 .
491 a.u.) and η = 0 . R e ( f ) - I m ( f ) n = 0.001 n = 0.004 n = 0.02 n = 0.1 FIG. 5. Pauli kernel, f P , for the free electron gas. Plots of f P with respect to ω are shown for different n (in a.u.) at q = 1 a.u. and η = 0 . f P has two cusps at the same position of the cusps of the Lindhardfunction for χ S , and the imaginary part has a wide feature with the same width as the peak in theLindhard function for χ S . Figure 5 compares the f P - ω dependence for different choices of densityvalues. f P gets closer to 0 as the density increases because the system becomes more Thomas-Fermilike [88]. On the other hand, at lower density, the non-zero features of f P are more pronounced,for example the feature between ω from 0.2 a.u. to 0.8 a.u. for n = 0 .
001 a.u..
C. Extending the Pauli kernel from the uniform electron gas to inhomogeneous systems
One way to extend Eq.(59) to general systems is to replace the Fermi wave vector k F with atwo-body Fermi wave vector ξ ( r , r (cid:48) ), such as ξ ( r , r (cid:48) ) = (cid:20) k γF ( r ) + k γF ( r (cid:48) )2 (cid:21) /γ , (61)or a geometric average ξ ( r , r (cid:48) ) = k / F ( r ) k / F ( r (cid:48) ) . (62)Eq.(61) is used in functionals of the noninteracting kinetic energy ( T S [ n ]) such as CAT[89] andWGC[24] and has the advantage of ξ ( r , r (cid:48) ) always being positive as long as the density is not 0at both r and r (cid:48) . However, it leads to a quadratic scaling algorithm because it is not possible toformulate the resulting equations in terms of convolution-like integrals where fast Fourier transformcan be applied. With the choice of Eq.(62) for the Fermi wavevectors, we propose the non-adiabaticcorrection to the Pauli kernel in the following form (See Appendix C for detailed derivation): f nad P ( ω, q ) = iπ (cid:18) ξ q + 6 qξ (cid:19) ω + π (16 − π )4 ξ q ω + π (16 − π )24 k F ω . (63)Note the last term in Eq.(63) is purely local (e.g. does not have a contribution from q ).To use the kernel in practical calculations, Casida matrix elements need to be computed, suchas, K ij = (cid:90) d r (cid:90) d r (cid:48) φ ( r (cid:48) ) φ ∗ j ( r (cid:48) ) f nad P ( ω, r , r (cid:48) ) φ ∗ ( r ) φ i ( r ) , (64)we need to treat the space dependence of k F in the following way. For last term in Eq.(63) whichis local, we apply the local-density approximation or LDA. Namely, k F ≡ k F ( r ) = (3 π n ( r )) / .9The LDA results in the following commonly employed integrals (e.g., for the so-called ALDA kernelin TD-DFT), K local ij = (cid:90) d r (cid:90) d r (cid:48) φ ( r (cid:48) ) φ ∗ j ( r (cid:48) ) f nad , local P ( ω, r ) δ ( r (cid:48) − r ) φ ∗ ( r ) φ i ( r ) . (65)For the non-local terms, for example, the second to the last term in Eq.(63), we use k F = k / F ( r ) k / F ( r (cid:48) ), and K ij = (cid:90) d r (cid:48) φ ( r (cid:48) ) φ ∗ j ( r (cid:48) ) k − / F ( r (cid:48) ) F − (cid:26) π (16 − π ) q ω F { k − / F ( r ) φ ∗ ( r ) φ i ( r ) } (cid:27) , (66)where F {·} and F − {·} represent forward and inverse Fourier transform, respectively. D. Pilot calculations involving the nonadiabatic Pauli kernel
To demonstrate the effect of the nonadiabatic Pauli kernel, we present several pilot calculationsinvolving the solution of the Dyson equation in Eq.(52) linking the boson and KS response func-tions. That is, we begin with the boson response function, χ B , and through Eq.(52) employing anapproximate Pauli kernel we recover an approximate χ S which we compare to the exact one.The role of the nonadiabatic Pauli kernel is to line up the poles of χ B with the poles of χ S .According to the examples provided in Figure 2 and Table II, the Pauli kernel should red shiftthe bosonic poles by as much as 4.5 eV and 2.8 eV for the silver nano-rod and sodium cluster,respectively. Therefore, we expect f p and its matrix elements with respect to occupied-virtualorbital products to be large in size and comparable to the orbital energy differences.Therefore, it is problematic that the simplest approximation to it [i.e., the TF functional withkernel f TF p ( r , r (cid:48) ) = C T F δ ( r − r (cid:48) ) n − ( r )] is strictly positive. From the single pole approximation(i.e., ω S = (cid:112) ω B + 2 Kω B , where ω S / B are the KS/boson orbital excitations and K is the matrixelement of the Pauli kernel), we see that the positivity of the kernel leads to a blue shift of thefirst excited state which is opposite to the sought behavior. For these reasons, we expect thenonadiabatic part of the kernel to play a very significant role.In the calculations presented below we start from the bosonic poles calculated by diagonalizingthe noninteracting boson Hamiltonian of Eq.(8), where the v P is determined by inversion usingEq.(55) and therefore it is exact within the numerical precision of the inversion procedure. Wethen compute the lowest-lying KS pole by solving the Casida equation [90] associated with theDyson equation Eq.(52) with an approximate Pauli kernel and compare to the exact KS polecomputed by diagonalization of the KS Hamiltonian.0 E n e r g y ( e V ) Na Boson poleKS poleAdiabatic Local term/1st column & rowFull correction/1st column & rowFull correction/Full matrix0 5 10 15 20Iteration02468 Ag FIG. 6. The convergence of the first KS pole with the non-adiabatic f P correction for an icosahedral Na cluster and a Ag nano-rod (structures in inset). The yellow solid line is the position of the bosonic pole, thegreen dotted line is the position of the KS pole calculated using only the adiabatic f P and the black dashedline is the position of the KS pole from diagonalizing the KS Hamiltonian. Red filled circles: uses only thelocal term and couples the first state with the others only. Blue half-filled circles: uses full correction andcouples the first state with the others only. Black empty circles: uses full correction and full coupling. We use the adiabatic TF as the adiabatic part of the Pauli kernel and Eq. (63) evaluated accord-ing to the prescriptions given in subsection IV C for the nonadiabatic part. The Dyson equationin Eq.(63) needs to be solved iteratively starting from the adiabatic result. The iterations evolveas follows: (1) at the n -th cycle we evaluate the KS response function via the Dyson equation, χ n S ( ω out n ) = χ B ( ω out n ) + χ B ( ω out n ) f P ( ω in n − ) χ n S ( ω out n ); (2) the kernel is updated with the value of thenew frequency ω in n = βω out n + (1 − β ) ω in n − , and the Dyson equation is solved again to find yetanother pole frequency. The β mixing parameter was determined adaptively from β = 1 to β = 0 . and a silver nanorod, Ag ,respectively. For the Na system, we carried out additional analyses: (1) we show the effect ofapplying the Pauli kernel on one row/column of the Casida matrix or on the full matrix (bluehalf-filled circles vs black empty circles) and (2) we show that applying only the space-local part1of the nonadiabatic Pauli kernel is not enough to substantially red shift the bosonic pole (red filledcircles).To ensure that the numerical artifacts in the inverted v P do not affect the final results, wealso performed the same calculations in Figure 6 with v P multiplied by a mask function m ( r ) =1 − / { n KS ( r ) /n cutoff ] } , where n KS ( r ) is the KS density used as the input density of theinversion, and n cutoff = 10 − a.u.. The mask function removed all artifacts from v P . The results ofthese calculations are almost identical to those in Figure 6.The above examples show that the proposed nonadiabatic correction to the Pauli kernel providesthe required red shift to the bosonic pole towards the KS pole, provided the kernel contains space-nonlocal terms. V. SUMMARY
In the first part of this work, we present orbital-free TD-DFT and derive properties (or con-ditions) that the Pauli energy, potential and kernel must satisfy. Furthermore, we provide proofsof theorems that also apply to conventional TD-DFT, such as the Runge-Gross and van Leuuwentheorems.In the second part of this work, we derive an approximation for the Pauli kernel based on theresponse properties of the uniform electron gas. The derived kernel is nonlocal in space as wellas nonlocal in time because it features an explicit frequency dependence. Pilot calculations showthat the nonadiabatic part of the kernel is much more important than the adiabatic part and thatspace-nonlocal terms in the kernel play a fundamental role. The proposed kernel is capable ofcorrectly red shifting the orbital free (bosonic) orbital excitation energies, bringing them close tothe KS orbital excitations. While the final result of the excitation energy of the lowest-lying KSexcited state for a sodium cluster and a silver nanorod are only in semiquantitative agreement withthe exact result, this work sheds light on a path to develop nonadiabatic Pauli kernels for modelingwith time-dependent orbital-free DFT systems out of equilibrium in a computationally cheap andsemiquantitatively accurate way.2
Appendix A: Derivation of Eq.(15)
We start with the equation of motion for the current density operator i ∂∂t j ( r , t )= (cid:68) φ B ( t ) (cid:12)(cid:12)(cid:12) [ˆ j ( r ) , ˆ H B ( t )] (cid:12)(cid:12)(cid:12) φ B ( t ) (cid:69) = (cid:68) φ B ( t ) (cid:12)(cid:12)(cid:12) ˆ j ( r ) ˆ H B ( t ) − ˆ H B ( t )ˆ j ( r ) (cid:12)(cid:12)(cid:12) φ B ( t ) (cid:69) = (cid:68) φ B ( t ) (cid:12)(cid:12)(cid:12) ˆ j ( r ) ˆ T − ˆ T ˆ j ( r ) + ˆ j ( r ) ˆ V B ( t ) − ˆ V B ( t )ˆ j ( r ) (cid:12)(cid:12)(cid:12) φ B ( t ) (cid:69) , (A1)where (cid:68) φ B ( t ) (cid:12)(cid:12)(cid:12) ˆ j ( r ) ˆ V B ( t ) − ˆ V B ( t )ˆ j ( r ) (cid:12)(cid:12)(cid:12) φ B ( t ) (cid:69) = 12 i N (cid:88) l (cid:90) d r . . . d r N φ ∗ B ( r , r , . . . , r N , t ) {∇ l [ δ ( r − r l ) v B ( r l , t )] + δ ( r − r l ) ∇ l v B ( r l , t ) − v B ( r l , t ) ∇ l δ ( r − r l ) − v B ( r l , t ) δ ( r − r l ) ∇ l } φ B ( r , r , . . . , r N , t )= N i { φ ∗ B ( r , t ) ∇ [ v B ( r , t ) φ B ( r , t )] − ∇ [ φ ∗ B ( r , t ) v B ( r , t ) φ B ( r , t )] + φ ∗ B ( r , t ) ∇ [ v B ( r , t ) φ B ( r , t )]+ ∇ [ φ ∗ B ( r , t ) v B ( r , t ) φ B ( r , t )] − φ ∗ B ( r , t ) v B ( r , t ) ∇ φ B ( r , t ) − φ ∗ B ( r , t ) v B ( r , t ) ∇ φ B ( r , t ) } = Ni φ ∗ B ( r , t ) φ B ( r , t ) ∇ v B ( r , t )= 1 i n ( r , t ) ∇ v B ( r , t ) , (A2)3and (cid:68) φ B ( t ) (cid:12)(cid:12)(cid:12) ˆ j ( r ) ˆ T − ˆ T ˆ j ( r ) (cid:12)(cid:12)(cid:12) φ B ( t ) (cid:69) = − i N (cid:88) l N (cid:88) k (cid:90) d r . . . d r N φ ∗ B ( r , r , . . . , r N , t ) {∇ l [ δ ( r − r l ) ∇ k ] + δ ( r − r l ) ∇ l ∇ k − ∇ k ∇ l δ ( r − r l ) − ∇ k [ δ ( r − r l ) ∇ l ] } φ B ( r , r , . . . , r N , t )= − N i { φ ∗ B ( r , t ) ∇ φ B ( r , t ) − ∇ [ φ ∗ B ( r , t ) ∇ φ B ( r , t )] + φ ∗ B ( r , t ) ∇ φ B ( r , t ) − φ ∗ B ( r , t ) ∇ φ B ( r , t ) + 3 ∇ [ φ ∗ B ( r , t ) ∇ φ B ( r , t )] − ∇ [ φ ∗ B ( r , t ) ∇ φ B ( r , t )] + ∇ [ φ ∗ B ( r , t ) φ B ( r , t )] − φ ∗ B ( r , t ) ∇ φ B ( r , t ) + 2 ∇ [ φ ∗ B ( r , t ) ∇ φ B ( r , t )] − ∇ [ φ ∗ B ( r , t ) ∇ φ B ( r , t )]+ ( N − { φ ∗ B ( r , t ) ∇ φ B ( r , t ) − ∇ [ φ ∗ B ( r , t ) φ B ( r , t )] + φ ∗ B ( r , t ) ∇ φ B ( r , t )+ ∇ [ φ ∗ B ( r , t ) φ B ( r , t )] − φ ∗ B ( r , t ) ∇ φ B ( r , t ) − φ ∗ B ( r , t ) ∇ φ B ( r , t ) } (cid:90) d r (cid:48) φ ∗ B ( r (cid:48) , t ) ∇ (cid:48) φ B ( r (cid:48) , t ) } = − N i { ∇ [ φ ∗ B ( r , t ) ∇ φ B ( r , t )] − ∇ [ φ ∗ B ( r , t ) ∇ φ B ( r , t )] + ∇ [ φ ∗ B ( r , t ) φ B ( r , t )] } = − N i {− ∇ φ ∗ B ( r , t )][ ∇ φ B ( r , t )] − ∇ φ ∗ B ( r , t )][ ∇ φ B ( r , t )] + ∇ [ φ ∗ B ( r , t ) φ B ( r , t )] } . (A3)Combining Eq. (A1) through Eq.(A3) we get ∂∂t j ( r , t ) = − n ( r , t ) ∇ v B ( r , t )) − F kin ( r , t ) , (A4)where F kin ( r , t ) = N { ∇ φ ∗ B ( r , t )][ ∇ φ B ( r , t )] + 4[ ∇ φ ∗ B ( r , t )][ ∇ φ B ( r , t )] − ∇ [ φ ∗ B ( r , t ) φ B ( r , t )] } . (A5) Appendix B: Derivation of Eq.(41)
Plug the Hamiltonian of the noninteracting boson system as ˆ O into Eq. (39), we obtaind E ( t )d t = (cid:42) ∂ ˆ H∂t (cid:43) = (cid:42) ∂ ˆ T∂t (cid:43) + (cid:42) ∂ ˆ V B ( t ) ∂t (cid:43) = (cid:90) d r ∂v B ( t ) ∂t n ( t ) . (B1)From Eq.(34) we obtain d E ( t )d t = d T B ( t )d t + (cid:90) d r v B ( r , t ) n ( r , t ) . (B2)4Combining Eq. (B1) and Eq. (B2) we getd T B ( t )d t = − (cid:90) d r ∂n ( r , t ) ∂t v B ( r , t ) . (B3) Appendix C: Derivation of non-adiabatic correction to the Pauli kernel
First we replace variables ω , q , η with their dimensionless counterparts ¯ ω = ω/qk F , ¯ η = q/ k F ,¯ γ = η/qk F . We can rewrite χ B and χ S χ B (¯ ω, ¯ η, ¯ γ ) = 6 π ¯ ηk F (cid:18) ω − ¯ η + i ¯ γ + 1 − ¯ ω − ¯ η − i ¯ γ (cid:19) − , (C1) χ S (¯ ω, ¯ η, ¯ γ ) = 2 π ¯ ηk F [Ψ (¯ ω − ¯ η + i ¯ γ ) − Ψ (¯ ω + ¯ η + i ¯ γ )] − . (C2)Using the Dyson equation Eq.(52) that relates χ B and χ S , we find the non-adiabatic and adiabaticPauli kernel f P and f P as following f P (¯ ω, ¯ η, ¯ γ ) = χ B (¯ ω, ¯ η, ¯ γ ) − χ S (¯ ω, ¯ η, ¯ γ ) , (C3) f P (¯ η ) = χ B (¯ ω = 0 , ¯ η, ¯ γ = 0) − χ S (¯ ω = 0 , ¯ η, ¯ γ = 0) . (C4)Expand f P and f P around ¯ η = 0 up to the second order of ¯ η , we obtain f P (¯ ω, ¯ η, ¯ γ ) = π k F (cid:26) (cid:34) − γ + 6 i ¯ γ ¯ ω + 3¯ ω + 2 i i + (¯ γ − i ¯ ω ) log( i ¯ γ +¯ ω − i ¯ γ +¯ ω ) (cid:35) + (cid:34) − − γ − i ( − ω )] [¯ γ − i (1 + ¯ ω )] [2 i + (¯ γ − i ¯ ω ) log( i ¯ γ +¯ ω − i ¯ γ +¯ ω )] (cid:35) ¯ η + O (¯ η ) (cid:27) , (C5) f P (¯ η ) = π k F (cid:20) −
83 ¯ η + O (¯ η ) (cid:21) . (C6)Take the limit ¯ γ → ω = 0 up to the second order of ¯ ω , we obtain f P (¯ ω, ¯ η ) = π k F (cid:20)(cid:18) −
83 ¯ η (cid:19) + iπ η )¯ ω + 112 (48 − π + 16¯ η − π ¯ η )¯ ω + O (¯ ω ) (cid:21) . (C7)The non-adiabatic correction to the Pauli kernel is f nad P (¯ ω, ¯ η ) = f P (¯ ω, ¯ η ) − f P (¯ η )= π k F (cid:20) iπ η )¯ ω + 112 (48 − π + 16¯ η − π ¯ η )¯ ω (cid:21) . (C8)Replace the variables back, we get f nad P ( ω, q ) = iπ (cid:18) k F q + 6 qk F (cid:19) ω + π (16 − π )4 k F q ω + π (16 − π )24 k F ω . (C9)5 ACKNOWLEDGMENTS