Disordered Crystals from First Principles II: Transport Coefficients
DDisordered Crystals from First Principles II: TransportCoefficients
Thomas D. K¨uhne , Julian Heske and Emil Prodan Chair of Theoretical Chemistry, University of Paderborn, Paderborn, Germany Department of Physics, Yeshiva University, New York, New York, USA
Abstract
This is the second part of a project on the foundations of first-principle calculations of the electron transport incrystals at finite temperatures, aiming at a predictive first-principles platform that combines ab-initio moleculardynamics (AIMD) and a finite-temperature Kubo-formula with dissipation for thermally disordered crystallinephases. The latter are encoded in an ergodic dynamical system (Ω , G , d P ), where Ω is the configuration spaceof the atomic degrees of freedom, G is the space group acting on Ω and d P is the ergodic Gibbs measure relativeto the G -action. We first demonstrate how to pass from the continuum Kohn-Sham theory to a discrete atomic-orbitals based formalism without breaking the covariance of the physical observables w.r.t. (Ω , G , d P ). Then weshow how to implement the Kubo-formula, investigate its self-averaging property and derive an optimal finite-volume approximation for it. We also describe a numerical innovation that made possible AIMD simulationswith longer orbits and elaborate on the details of our simulations. Lastly, we present numerical results on thetransport coefficients of crystal silicon at different temperatures. Contents
Preprint submitted to Elsevier August 3, 2020 a r X i v : . [ phy s i c s . c o m p - ph ] J u l Acknowledgements 32
1. Introduction
As the scale of the fabrication processes of electronic components is continuously re-duced, the quantum mechanical aspects of the charge transport become more importantand ab-initio quantum simulations will be required for an accurate and predictive charac-terization. Since most electronic components operate at room and higher temperatures,these ab-initio simulations have to take into account the thermal motion of the atoms.Since the dynamics of the electrons is orders of magnitude faster than that of the ioniccores, the quantum dynamics of the electrons takes place in a highly disordered environ-ment. This can result in qualitatively different dynamical behaviors, notably the absenceof quantum diffusion or Anderson localization [1], that can not be captured by empiricalmodels or idealistic zero temperature simulations. This and other effects will be investi-gated from first-principles in this work. To the best of our knowledge, the work reportedhere is the first attempt to simulating quantum charge transport at finite-temperaturefrom first-principles.Most of us think of crystals as condensed phases of matter, where the atoms areperiodically arranged in space. However, the crystalline phase persists all the way to themelting point so clearly that oversimplifying picture is highly misleading. In fact, defin-ing the crystalline phase is a deep and highly non-trivial problem in condensed matterphysics. On the formalism side, the works of Bellissard on the homogeneous phases ofmatter represent a milestone [2, 3]. They taught us that, in crystals, the space groupsymmetry G manifested at zero temperature is replaced at finite temperatures by an er-godic G -action w.r.t. the Gibbs measure on the space Ω of thermally disordered atomicconfigurations. Furthermore, the invariance w.r.t. G of the electronic Hamiltonians man-ifested at zero temperature is replaced at finite temperatures by the covariance w.r.t. the G -action. Ergodicity and convariance w.r.t. the G -action explain why the measurementsof the macroscopic physical observables, including the transport coefficients, do not fluc-tuate from one configuration to another and why the symmetry w.r.t. the full spacegroup is restored at the macroscopic level. For example, the latter is manifested in the G -symmetric X-ray diffraction patterns observed all the way to the melting point [4].Another manifestation is the stability of the topological phases of matter stabilized bypoint symmetries in conditions where thermal disorder breaks these symmetries [5].In our previous work [6], we took the task of quantifying the ergodic dynamical sys-tem (Ω , G , τ, d P ) that defines a crystalline phase. Using crystal silicon (Si) as a workingexample, we devised an algorithm that extracts this data from the output of conventional ab-initio molecular dynamics (AIMD) simulations [7, 8]. In particular, we were able toquantify and parametrize the Gibbs measure for crystalline Si at various temperatures.In this work, our focus is mostly on the electronic degrees of freedom, which are simu-lated with hybrid Gaussian-plane wave based density functional theory (DFT) electronicstructure codes [9]. In the first part of our work, we demonstrate how to generate effec-tive lattice models that encode the entire output of the electronic structure codes andwhere the covariance w.r.t. (Ω , G , τ, d P ) is explicitly manifested. Particular attentionwill be dedicated to the tight-binding expressions of the Kohn-Sham (KS) Hamiltonian,position and charge current operators. 2or the charge transport, we adopt the non-commutative Kubo-formula derived bySchulz-Baldes and Bellissard [10–12]. One extremely important aspect of their formalismis that it includes dissipation. More precisely, given a dissipation mechanism encoded ina scattering operator, the formalism produces a dissipation super-operator that is organi-cally incorporated in the Kubo-formula (see section 4.2). Various dissipation mechanismsand their corresponding super-operators have been analyzed in [13] and they certainlycan be evaluated from first principles. One should be aware that dissipation has an im-portant role in shaping the I-V characteristics of both metals and semiconductors andthis is why a Kubo-formula that incorporates dissipation is so valuable.The Kubo-formula derived by Schulz-Baldes and Bellissard has been numerically im-plemented in the past for disordered tight-binding model Hamiltonians [14–19] and othertypes of aperiodic Hamiltonians [20–22]. One of the main findings of these works isthe rapid convergence of the results with system size. For example, in systems withknown quantized transport coefficients, such as 2-dimensional Hall systems, the non-commutative Kubo-formula reproduced the quantization with two digits of precisioneven on small 10 ×
10 highly disordered tight-binding lattices. This is a convincing factthat this approach is highly suited for the applications we seek in this work, given thatthe super-cells that can be handled by first-principles simulations are inherently small.The simulations we report here for crystalline Si at different temperatures are prelim-inary and certainly not converged w.r.t. either the system’s size, or the atomic orbitalbasis, but they are certainly converged w.r.t. the thermal disorder sampling. Also, thedissipation super-operator is treated in the relaxation time approximation where it be-comes proportional to the identity map. The simulations produced expected outputs forthe available electronic structures and enabled us to test several important qualitativeaspects of the charge transport. One aspect is the formation of a dynamical band gapwhere the quantum diffusion is absent and this dynamical gap was found to be muchlarger than the spectral gap. The former defines the reference for the activated behaviorof the conductivity, while the latter for the charge carriers. Since these are two differentreference energy levels, the Anderson localization phenomenon can lead to substantialquantitative effects that were overlooked so far. We also found that the conductivitytensor is extremely sensitive to the dissipation relaxation time. Given this sensitivity,we believe that the prevalent dissipation mechanism in crystalline Si at room tempera-ture can be identified with high precision by future simulations which incorporate firstprinciples dissipation super-operators.Based on previous tight-binding model simulations, we initially estimated that at least1000 disordered atomic configurations will be necessary and, as such, we performed largetime scale AIMD simulations to acquire that amount of data. However, our calculationsrevealed that the average over the atomic configurations of the conductivity tensor canbe achieved with a relatively small number of configurations, which can be as low as 50.In fact, with reasonable level of dissipation, the thermal fluctuations are almost entirelysuppressed for the largest crystal we simulated, which is a direct manifestation of theself-averaging property of the Kubo-formula. This finding assures us that, in the futuresimulations, we can reduce the time scale of AIMD simulations, hence, enabling us tofurther increase the crystal size and to better optimize the orbital basis.3 . Thermal Disorder from First Principles
In order to fix our notations and provide the context for the present calculations, webriefly recall our main results reported in [6]. Therein, we describe the ergodic dynamicalsystem (Ω , G , τ, d P ), which completely characterizes the crystalline phase of Si at finite-temperature, where Ω is the atomic configuration space, G is the space group, d P is theGibbs measure and τ is an ergodic action of G on Ω. The crystal structure of Si is summarized in Fig. 2.1. Its space group is G = F d ¯3m[23, 24], whose structure is summarized by the following exact sequence of groups1 → B → G → P → , (2.1)capturing the extension of the point group P ⊂ O (3) by the group of discrete translations B . The latter can be pictured as the Bravais lattice of the crystal (hence our notation B ), i.e. the discrete sub-group of R defined by the centers of the primitive cells B = (cid:8) n a + n a + n a , n = ( n , n , n ) ∈ Z (cid:9) , (2.2)with the generators a i supplied in Fig. 2.1. The point group P of crystalline Si is O h ,the full symmetry group of the cube.Let us recall that a space group is called symmorphic if the exact sequence (2.1) issplit. Silicon’s cubic-diamond lattice is an example of a non-symmorphic space group.Nevertheless, every element g of G can be presented in the form g = ( p | a ), with p ∈ P and a ∈ R . Note that for a symmorphic space group, a can be always drawn from B , butthis is not the case here. Such space group elements act on the points of the Euclideanspace as ( p | a ) x = p x + a , x ∈ R . (2.3)They also act on any subset L of the Euclidean space, such as a lattice, via g L = { g x , x ∈ L } . (2.4)The multiplication of the elements takes the form( p | a )( p (cid:48) | a (cid:48) ) = ( pp (cid:48) | pa (cid:48) + a ) (2.5)and the inverse of an element is( p | a ) − = ( p − | − p − a ) . (2.6)The ideal or zero temperature Si lattice will be denoted by L . This lattice is leftinvariant by the space group G . In fact, the asymmetric unit cell of the diamond cubicstructure contains a single atom [23, 24], which means that the entire lattice can bereconstructed from one single point by acting with the full space group: L = (cid:8) g · x , g ∈ G (cid:9) . While x can be any point of the Euclidean space, we will fix x at the origin.4 a a a Figure 2.1:
Si crystalizes in a diamond cubic lattice (Fd-3m), whose conventional unit cell isshown in this diagram. This cubic unit cell is symmetric to the full point group and containseight Si atoms. The diamond cubic lattice results from the inter-penetration of two face-centeredcubic (fcc) lattices. The fcc lattice can be generated by translating a primitive cell that containsjust one atom. Hence, silicon’s diamond cubic lattice can be generated by translating the sameprimitive cell, but with one additional Si atom inside it. This primitive cell is shown in red,together with its two atoms (blue disks) and the generating primitive vectors a = b (ˆ z + ˆ y ), a = b (ˆ x + ˆ z ) and a = b (ˆ y + ˆ x ), with b = 5 .
431 ˚A. The magnitude of the primitive vectors is a = b/ √ When the temperature is finite, the Si atoms undergo a thermal motion and theinstantaneous snapshots of the Si lattice can be labeled by a configuration space Ω. In[6], Ω was found to be well represented by a product of identical ballsΩ = (cid:89) g ∈ G B g , B g = B , (2.7)where a point ω = { ω g } g ∈ G of Ω encodes the displacements of the atoms from their equi-librium positions. Thermal motion defines an ergodic dynamical system ω ( t ) ( t = time)over Ω and an instantaneous snapshot of crystalline Si supplies a thermally disorderedlattice L ω = (cid:8) g x + ω g , g ∈ G (cid:9) ⊂ R . (2.8)For these disordered lattices, the invariance of L under the space group is replaced bythe covariance relation g L ω = L τ g ω , ∀ ω ∈ Ω , g ∈ G . (2.9)5he action τ of the space group on Ω, appearing above, can be computed as follows. If g = ( p | a ) ∈ G , then g L ω = { g ( g (cid:48) x + ω g (cid:48) ) , g (cid:48) ∈ G } (2.10)= { p ( g (cid:48) x ) + p ω g (cid:48) + a , g (cid:48) ∈ G } . After regrouping, p ( g (cid:48) x ) + p ω g (cid:48) + a = gg (cid:48) x + p ω g (cid:48) , (2.11)and, after the change of variable g (cid:48) → g − g (cid:48) , we have g L ω = { g (cid:48) x + p ω g − g (cid:48) , g (cid:48) ∈ G } = L τ g ω . (2.12)We now can identify the action as τ g ω = τ g { ω g (cid:48) } g (cid:48) ∈ G = { ω (cid:48) g (cid:48) } g (cid:48) ∈ G , ω (cid:48) g (cid:48) = p ω g − g (cid:48) . (2.13)One can verify that τ g τ g = τ g g , as it should be for a group action.The Gibbs measure over the configuration space Ω can be computed from the atomicorbits in an AIMD simulation, more precisely, from the histograms encoding the numberof times an orbit intersects the elementary volumes of Ω. In [6], the Gibbs measure ofthe Si crystal was found to be extremely well characterized by a multivariate normaldistribution of zero meand P ( ω ) = ρ ( ω )d ω, ρ ( ω ) = √ Det(2 π ˆΣ) e − ω T ˆΣ − ω , (2.14)where ω is seen here as a 1-column matrix and the variance matrix ˆΣ was quantified in[6] as a function of temperature. The Gibbs measure is invariant and ergodic w.r.t. the τ -action. In fact, the Gibbs measure found in [6] is ergodic relative to the subgroup B of translations, which is in fact a generic property of homogeneous systems at thermo-dynamic equilibrium [25, Chap. 6]. The crystalline phase of Si at finite-temperature isentirely defined by the ergodic dynamical system (Ω , G , τ, d P ).The observations of the last paragraph will play an important role for the self-averaging properties of the transport coefficients. Let us stress again that, due to thewell separated scales in the dynamics of the atomic and electronic degrees of freedom,the quantum state of the electrons evolves in a static atomic potential. This is a ther-mally disordered potential and, as we shall see, the physical observables, such as theHamiltonians or charge currents, become indexed by points of the configuration spaceΩ. In this new context, the notion of a symmetric observable is replaced by that of acovariant observable. The macroscopic measurements of these observables, however, areindependent of the thermally disordered configuration. This remarkable property is aconsequence of the covariance and of the ergodic character of the Gibbs measure. The electron-phonon and electron-electron scattering processes are rare and sudden dynamicalevents, which are included via Poisson processes as explained later. . Tight-Binding Form of the Physical Observables Our goal for this section is to formulate discrete representations of the Hamiltoni-ans and other physical observables in the context of Gaussian-based implementations ofthe KS program. Special attention will be given to the transformation of the physicalobservables under the space group G of the crystal. The formally exact KS theory for condensed matter systems, at its most fundamentallevel, is formulated over the Hilbert space L ( R ) of square integrable wave functions[26]. For the Si crystal in a thermally disordered configuration ω , the KS-Hamiltoniantakes the form H ω KS = − (cid:126) m ∇ r − (cid:88) x ∈ L ω Ze | r − x | + V xc [ n ω ]( r ) , (3.1)where V xc is a local potential encoding the exchange and correlation (XC) effects. Thelatter has a functional dependence on the electron density n ω ( r ), which is to be deter-mined self-consistently. As the notation suggests, the electron density has a dependenceon the atomic configuration ω . In fact, this becomes even more apparent if we reformulate(3.1) as a fixed point problem n ω ( r ) = (cid:104) r | Φ FD (cid:0) H ω KS ; T, µ (cid:1) | r (cid:105) , (3.2)where Φ FD is the Fermi-Dirac distribution at temperature T and chemical potential µ [27]. In the following, we assume that this equation has a unique solution for almost allthermally disordered configuration (see [28] and [29]).The starting point of our study is the covariant property of the KS-Hamiltonian underthe space group transformations. To understand the origin of this property, we need to goall the way to the Euclidean group E of transformations and recall that the XC potentialenjoys the following property V xc [ n ◦ e ]( e − r ) = V xc [ n ]( r ) , e ∈ E , (3.3)for any density function n and point r ∈ R , which can be inferred from the universalityand uniqueness assumptions on V xc [30]. Certainly, this can be verified directly for thelocal density approximation (LDA) to V xc [31]. The action of E on R lifts to a unitaryaction on the Hilbert space L ( R ) via (cid:0) T e ψ (cid:1) ( r ) = ψ ( e − r ) , e ∈ E , ψ ∈ L ( R ) . (3.4)Now, recall that G is just a subgroup of the Euclidean group, hence (3.4) describes theaction of G as well. Then, under such unitary actions, the KS-Hamiltonian behaves as T g H ω KS T † g = − (cid:126) m ∇ g − r − (cid:88) x ∈ L ω Ze | g − r − x | + V xc [ n ω ]( g − r ) (3.5)= − (cid:126) m ∇ r − (cid:88) x ∈ L g ω Ze | r − x | + V xc [ n ω ◦ g − ]( r ) . (3.6)7e learn from here that, if n ω is the solution of (3.2) for configuration ω , then n ω ◦ g − is the solution of (3.2) for configuration g ω . In other words, the self-consistent solutionsof the KS equations enjoy the covariant property n g ω = n ω ◦ g − . (3.7)In turn, this assures us that the converged KS-Hamiltonian satisfies the covariant relation T g H ω KS T † g = H τ g ω KS , ∀ g ∈ G . (3.8)It will be extremely important to preserve this characteristics in our tight-binding ap-proximation. As we already mentioned, (3.8) together with the ergodicity of the spacegroup action ensure the self-averaging of the transport coefficients. In Gaussian-based approaches, the atom located at position x ∈ L ω carries a finite-dimensional local Hilbert space H x = Span (cid:8) φ n ( r − x ) , n = 1 , . . . N (cid:9) , (3.9)where φ n : R → C are optimized atomic orbitals (see section 5.2 for details). It isimportant to realize that the same set of functions φ n are used for all x ∈ L ω . The totalHilbert space for the Gaussian-based computations is the linear subspace H ω = Span (cid:8) H x , x ∈ L ω (cid:9) ⊂ L ( R ) . (3.10)As the notation suggests, this subspace depends on the configuration ω ∈ Ω of the atoms.As we shall see, it is isomorphic to the tight-binding Hilbert space C N ⊗ (cid:96) ( L ω ) = Span (cid:8) ξ ⊗ | x (cid:105) , ξ ∈ C N , x ∈ L ω (cid:9) (3.11)of square summable linear combinations of ξ ⊗ | x (cid:105) basis vectors. The scalar product forthis space is defined by the orthonormality condition (cid:104) x | x (cid:48) (cid:105) = δ x , x (cid:48) , ∀ x , x (cid:48) ∈ L ω . (3.12)All our physical observables will be mapped over this tight-binding Hilbert space and allthe calculations will be ultimately performed on C N ⊗ (cid:96) ( L ω ).Our goal for this section is to explain in details how to transfer the observablesbetween the Hilbert spaces. We start with the consideration of the overlap coefficients S ij xx (cid:48) ( ω ) = (cid:90) R d r φ ∗ i ( r − x ) φ j ( r − x (cid:48) ) , (3.13)which can be found among the outputs of standard AIMD simulations. Using thesecoefficients, we form the self-adjoint, positive and invertible operator S ω : C N ⊗ (cid:96) ( L ω ) → C N ⊗ (cid:96) ( L ω ) , (3.14) S ω = (cid:88) x , x (cid:48) ∈ L (cid:98) S x , x (cid:48) ( ω ) ⊗ | x (cid:105)(cid:104) x (cid:48) | , (cid:98) S x , x (cid:48) ( ω ) is the overlap matrix with the entries S ij x , x (cid:48) ( ω ) defined in (3.13). Then,the isomorphism between H ω and C N ⊗ (cid:96) ( L ω ) is supplied by the unique linear map U ω that acts on the generators as φ n ( r − x ) (cid:55)→ (cid:112) S ω ξ n ⊗ | x (cid:105) , n = 1 , . . . N. (3.15)Above, ξ n ∈ C N is a column vector, whose entries are one at position n and zero for allothers, whereas √ S ω is the square root operator defined via the functional calculus. Letus verify that the map indeed preserves the scalar product. We have (cid:0)(cid:112) S ω ξ i ⊗ | x (cid:105) , (cid:112) S ω ξ j ⊗ | x (cid:48) (cid:105) (cid:1) = (cid:0) ξ i ⊗ | x (cid:105) , S ω ξ j ⊗ | x (cid:48) (cid:105) (cid:1) (3.16)= ξ Ti (cid:98) S x , x (cid:48) ( ω ) ξ j = S ij x , x (cid:48) ( ω ) . As a consequence, (cid:0)(cid:112) S ω ξ i ⊗ | x (cid:105) , (cid:112) S ω ξ j ⊗ | x (cid:48) (cid:105) (cid:1) = (cid:90) R d r φ ∗ i ( r − x ) φ j ( r − x (cid:48) ) , (3.17)for all i, j = 1 , . . . , N , as desired.To preserve the covariance of the physical observables w.r.t. the space group G andthe disordered configurations, it is important to choose the atomic orbitals as such thatthey span a linear space, which is closed under representations of the O (3) group. Thus,all φ n are assumed to transform under rotations as( φ , . . . , φ N ) ◦ r − = ( φ , . . . , φ N ) (cid:98) D ( r ) , r ∈ O (3) ⊂ E , (3.18)where { (cid:98) D ( r ) , r ∈ O (3) } is a family of N × N matrices supplying a N -dimensional unitaryrepresentation of the rotation group that is not necessarily irreducible (see Section 5.2for details). Then, if g = ( p | a ) ∈ G , the overlap matrix satisfies the relation (cid:98) D ( p ) † (cid:98) S g x , g x (cid:48) ( τ g ω ) (cid:98) D ( p ) = (cid:98) S x , x (cid:48) ( ω ) , (3.19)for any x and x (cid:48) in L ω . This is an important relation for which we provide the derivationbelow. Indeed, both g x and g x (cid:48) belong to g L ω = L τ g ω and S ij g x , g x (cid:48) ( τ g ω ) = (cid:90) R d r φ ∗ i ( r − g x ) φ j ( r − g x (cid:48) ) (3.20)= (cid:90) R d ( g r ) φ ∗ i ( g r − g x ) φ j ( g r − g x (cid:48) ) . Since g r − g x = p ( r − x ) and d ( g r ) = d r , we can continue as S ij g x , g x (cid:48) ( τ g ω ) = (cid:90) R d r φ ∗ i (cid:0) p ( r − x (cid:1) φ j (cid:0) p ( r − x (cid:48) (cid:1) (3.21)= (cid:90) R d r D ( p − ) ∗ ki φ ∗ k ( r − x ) φ s ( r − x (cid:48) ) D ( p − ) sj = D ( p ) ik (cid:90) R d r φ ∗ k ( r − x ) φ s ( r − x (cid:48) ) D ( p ) ∗ js , ω ∈ Ω, we define a Hilbert space isomorphism between C N ⊗ (cid:96) ( L ω ) and C N ⊗ (cid:96) ( L τ g ω ) as T g (cid:0) ξ ⊗ | x (cid:105) (cid:1) = (cid:98) D ( p ) ξ ⊗ | g x (cid:105) , g = ( p | a ) ∈ G , x ∈ L ω . (3.22)Note that we use here the same notation as in (3.4) because these two maps can beeasily differentiated from the context. Now, by examining the rule of multiplication in(2.5) for the space group, it is immediate to see that T respect this binary operation: T g T g (cid:48) = T gg (cid:48) . Furthermore, it follows directly from (3.19) and the definition (3.14) that T g S ω T † g = S τ g ω , g ∈ G . (3.23)Indeed, if g = ( p | a ), then T g S ω T † g = (cid:88) x , x (cid:48) ∈ L ω (cid:98) D ( p ) (cid:98) S x , x (cid:48) ( ω ) (cid:98) D ( p ) † ⊗ | g x (cid:105)(cid:104) g x (cid:48) | (3.24)= (cid:88) x , x (cid:48) ∈ L τ g ω (cid:98) D ( p − ) † (cid:98) S g − x , g − x (cid:48) ( ω ) (cid:98) D ( p − ) ⊗ | x (cid:105)(cid:104) x (cid:48) | . Furthermore, from (3.19), (cid:98) D ( p − ) † (cid:98) S g − x , g − x (cid:48) ( ω ) (cid:98) D ( p − ) = (cid:98) S x , x (cid:48) ( τ g ω ) , (3.25)hence (3.23) follows.We conclude this section with the observation that the map U ω , defined in (3.15),satisfies the covariant relation U τ g ω T g = T g U ω . Indeed, for x ∈ L ω , φ j ( g − r − x ) = φ j (cid:0) p − ( r − g x ) (cid:1) = N (cid:88) k =1 φ k ( r − g x ) D ( p ) kj , (3.26)while T g (cid:112) S ω ξ j ⊗ | x (cid:105) = (cid:113) S τ g ω (cid:98) D ( p ) ξ j ⊗ | g x (cid:105) = N (cid:88) k =1 (cid:113) S τ g ω ξ k ⊗ | g x (cid:105) D ( p ) kj . (3.27)Then, by applying rule (3.15) on each terms of the two sums, one can convince oneselfthat we have the following correspondence φ j ( g − r − x ) (cid:55)→ T g (cid:112) S ω ξ j ⊗ | x (cid:105) , ∀ x ∈ L , (3.28)under the U map. Let A be an operator defined over L ( R ). Our goal here is to investigate how to definea canonical approximation as an operator A ω over the effective Hilbert space C N ⊗ (cid:96) ( L ω ).10he natural requirement is the matching of all the available matrix elements under the U ω map (3.15), i.e. (cid:104) φ n ( · − x ) | A | φ m ( · − x (cid:48) ) (cid:105) = (cid:104) n, x | (cid:112) S ω A ω (cid:112) S ω | m, x (cid:48) (cid:105) , (3.29)for all n, m = 1 , N and x , x (cid:48) ∈ L ω . For convenience, above and throughout, we use thenotation | n, x (cid:105) for ξ n ⊗ | x (cid:105) . Henceforth, let (cid:98) A x , x (cid:48) ( ω ) be the matrix with the entries A n,m x , x (cid:48) ( ω ) = (cid:90) R d r φ ∗ n ( r − x )( Aφ m )( r − x (cid:48) ) , (3.30)which is just the explicit form of the coefficients appearing in the left side of (3.29). Weform first the operator (cid:101) A ω = (cid:88) x , x (cid:48) ∈ L ω (cid:98) A x , x (cid:48) ( ω ) ⊗ | x (cid:105)(cid:104) x (cid:48) | , (3.31)over C N ⊗ (cid:96) ( L ω ). Then, the solution to (3.29) is supplied by A ω = S − ω (cid:101) A ω S − ω , (3.32)as it readily follows from a direct calculation. We call (3.32) the canonical tight-bindingoperator associated to the operator A that is defined in the continuum KS theory. Notethat under this correspondence, the identity operator is sent to the identity operator.Now assume that the continuum observable depends on ω in a covariant fashion. Insuch a case, we can repeat the calculations leading to (3.19) to prove (cid:98) D ( p ) † (cid:98) A g x , g x (cid:48) ( τ g ω ) (cid:98) D ( p ) = (cid:98) A x , x (cid:48) ( ω ) . (3.33)This automatically implies that (cid:101) A ω is a covariant operator under the space group trans-formations and, since A ω in (3.32) is a product of covariant operators, A ω is also acovariant operator: T g A ω T † g = A τ g ω . (3.34)Below, we apply this standard procedure to several observables of interest.As we learned in section 3.1, the continuum KS-Hamiltonian is a covariant observable.Furthermore, among the standard outputs of AIMD simulations are the matrix elements W ij x , x (cid:48) ( ω ) = (cid:90) R d r φ ∗ i ( r − x )( H ω KS φ j )( r − x (cid:48) ) . (3.35)This is precisely the data one needs to define the tight-binding Hamiltonian. Followingthe above procedure, we define first the operator (cid:101) H ω = (cid:88) x , x (cid:48) ∈ L ω (cid:99) W x , x (cid:48) ( ω ) ⊗ | x (cid:105)(cid:104) x (cid:48) | , (3.36)which then supplies the tight-binding expression of the KS-Hamiltonian H ω = S − ω (cid:101) H ω S − ω , T g H ω T † g = H τ g ω , g ∈ G . (3.37)11e now focus on the position operator X . At the continuum level of the theory, thematrix elements of the position operator are R ij x , x (cid:48) ( ω ) = (cid:90) R d r φ ∗ i ( r − x ) r φ j ( r − x (cid:48) ) . (3.38)Note that these matrix elements depend too on the disordered configuration. They,however, satisfy a different covariant relation (cid:98) D ( g ) † (cid:98) R g x , g x (cid:48) ( τ g ω ) (cid:98) D ( g ) = p (cid:98) R x , x (cid:48) ( ω ) + a (cid:98) S x , x (cid:48) ( ω ) , (3.39)for all g = ( p | a ) ∈ G . The above relation follows from an exercise similar to thatbelow (3.19). Before going any further, let us explain the notation. Note that (cid:98) R x , x (cid:48) ( ω )is actually a 3-component column vector with matrices as entries. Then, p in frontof it, which is an ordinary 3 × a in the second term is viewed as an ordinary 3-component vector from R such that a (cid:98) S x , x (cid:48) ( ω ) becomes a 3-component vector with matrix entries. Now, as before,we define an operator on C N ⊗ (cid:96) ( L ω ) (cid:101) R ω = (cid:88) x , x (cid:48) ∈ L (cid:98) R x , x (cid:48) ( ω ) ⊗ | x (cid:105)(cid:104) x (cid:48) | , (3.40)which satisfies the covariance relation T g (cid:101) R ω T † g = p − (cid:101) R τ g ω + ( p − a ) S τ g ω , (3.41)as it follows directly from (3.39). Then X ω = S − ω (cid:101) R ω S − ω (3.42)maps the position operator from L ( R ) to C N ⊗ (cid:96) ( L ω ). Furthermore, the mappedposition operator satisfies the covariance relation T g X ω T † g = p − X τ g ω + ( p − a ) I, g = ( p | a ) ∈ G , (3.43)where I is the identity operator. The above relation follows directly from (3.41).Note that, although X ω is not entirely a covariant operator, the commutator [ X ω , A ω ]is covariant whenever A ω is, i.e. T g [ X ω , A ω ] T † g = p − [ X τ g ω , A τ g ω ] , g = ( p | a ) ∈ G . (3.44)This will become relevant when we will analyze the charge current operator. Over the Hilbert space L ( R ), the trace per volume of a bounded operator A withcontinuous kernel (cid:104) r | A | r (cid:48) (cid:105) is defined asTr V { A } = lim V → R V (cid:90) V d r (cid:104) r | A | r (cid:105) , (3.45)12here, for consistency with the space group, we require that the limit be taken over finitevolumes V that are invariant under the point group action. Our goal here is to supplyits canonical translation over the effective Hilbert space C N ⊗ (cid:96) ( L ω ). For this, let A ω be the tight-binding operator associated to A . We claim thatTr V { A ω } = 1 V lim V →∞ | L ω ∩ V | (cid:88) x ∈ L ∩ V N (cid:88) n =1 (cid:104) n, x | A ω | n, x (cid:105) (3.46)supplies the canonical expression. Above, V is the volume per Si atom, which is justhalf of that of the primitive cell, and | · | denotes the cardinal of a set.Indeed, let us note that (cid:82) V d r (cid:104) r | A | r (cid:105) coincides with the trace of A , when A is re-stricted over L ( V ). This trace can be alternatively computed as (cid:80) i (cid:104) ψ i | A | ψ i (cid:105) , with { ψ i } being an arbitrary orthonormal basis of L ( V ). But, up to errors that are irrelevant inthe thermodynamic limit and when N is large, the finite-volume trace can be computedusing the partial basis { U ∗ ω | n, x (cid:105)} n =1 ,N x ∈ L ω ∩ V . As a consequence, if N a is the total numberof atoms in V , thenTr V { A } = 1 V lim N a →∞ N a N (cid:88) n =1 (cid:88) x ∈ L ∩ V (cid:104) n, x | U ω AU ∗ ω | n, x (cid:105) , (3.47)which coincides with (3.46).The trace per volume, which is defined in (3.46), is a genuine trace over the al-gebra of operators we encounter in this work. For example, it displays the standardproperty Tr V { A ω B ω } = Tr V { B ω A ω } . An extremely important property of Tr V is theself-averaging when evaluated on covariant operators, i.e. those operators satisfying therelations T g A ω T † g = A τ g ω . (3.48)Indeed, using the invariance of the trace under conjugations, we haveTr V { A ω } = 1 | H | (cid:88) g ∈ H ⊂ G Tr V { T g A ω T † g } = 1 | H | (cid:88) g ∈ H ⊂ G Tr V { A τ g ω } , (3.49)where H is a finite subset of G invariant to the point group. Since τ acts ergodicallyover Ω, in the limit H → G , Birkhoff’s ergodic theorem assures us that the last termcoincides with the ensemble average [32]. Hence,Tr V { A ω } = (cid:90) Ω d P ( ω ) Tr V { A ω } . (3.50)Let us point out that intensive thermodynamic variable as measured in laboratories, suchas the transport coefficients, are all computed as traces per volumes of covariant observ-ables. The aforementioned self-averaging property assures us that these macroscopicvariables do not fluctuate from one disordered configuration to another, as long as thecorresponding physical observables are covariant. This is the main reason why we payspecial attention to the covariant properties of the physical observables in our theory.13 . Transport Coefficients With the mappings from the previous section, we can formulate the theory of quantumcharge transport directly on the Hilbert space C N ⊗ (cid:96) ( L ω ). The goal of this section isto supply the key elements of this theory and to formulate the Kubo-formula for theconductivity tensor. The purpose of this section is to review the theory of charge transport in the presenceof dissipation, as developed by Schulz-Baldes and Bellissard [10–12].Let us recall that the physical observable corresponding to the 3-component vectorof the electron charge current density is J ω = − eı (cid:126) [ X ω , H ω ] , (4.1)where e = 1 . × − C is the charge of the electron. Based on the last remark in section3.3, J ω is a covariant operator, i.e. T g J ω T g = p − J τ g ω . (4.2)Under the action of an externally applied electric field E , the measured current-densityis j E = lim T →∞ T (cid:90) T d t Tr V (cid:8) J ω ρ ω ( t ) (cid:9) , (4.3)where ρ ω ( t ) is the time-evolved density matrix. The time evolution is w.r.t. the time-dependent Hamiltonian H ( t ) = H ω + e E · X ω + V ω ( t ) , (4.4)which incorporates the externally applied electric field E , as well as dissipation via thescattering potential V ω ( t ) = (cid:88) j ∈ Z δ ( t − t j ) W ω , (4.5)with W ω assumed to be covariant. The collision times η = { t j } j ∈ Z are generated via aPoisson process with fixed collision time-scale τ c . Such processes are known to be self-averaging, hence the time and the space averages in (4.3) do not depend on the particularrealization of the collision times, nor on the disordered configuration. In other words, j E defined in (4.3) is a genuine macroscopic thermodynamic coefficient.For the reason state above, one can use in (4.3) an effective quantum time evolution,which is averaged over the Poisson processes η . A computation of this average can befound in [14]. It takes the form U eff ( t ) AU eff ( t ) ∗ = e − t (cid:126) (Γ ω +L E ,ω ) [ A ] , (4.6)where Γ is the collision super-operator, acting on the physical observables asΓ ω [ A ] = (cid:126) τ c ( A − e ı (cid:126) W ω A e − ı (cid:126) W ω ) , (4.7)14nd L E ,ω is the super-operatorL E ,ω [ A ] = ı [ H ω , A ] − e E · ı [ X ω , A ] . (4.8)The electrons are assumed initially at the thermal equilibrium, hence the initial den-sity matrix takes the form ρ ω ( t = 0) = Φ FD ( H ω ; T, µ ) , (4.9)where Φ FD ( (cid:15) ; T, µ ) is the Fermi-Dirac distribution at temperature T and µ the chemicalpotential. The density matrix is evolved via the time propagator (4.6), hence ρ ω ( t ) = U ( t ) Φ FD ( H ω ; T, µ ) U ( t ) ∗ . (4.10)Since, the two parameters T and µ are kept fixed, we will omit writing them explicitly. With the inputs supplied in the previous section, (4.3) can be evaluated explicitly: J ω E = e (cid:126) Tr V (cid:110)(cid:2) X ω , H ω (cid:3)(cid:0) Γ ω + L E ,ω (cid:1) − (cid:2) E · X ω , Φ FD ( H ω ) (cid:3)(cid:111) . (4.11)In the linear regime, this leads to a Kubo-formula with dissipation for the conductivitytensor σ αβ ( T, µ ; ω ) = − πG Tr V (cid:110)(cid:2) X αω , H ω (cid:3) (Γ ω + L ω ) − (cid:2) X βω , Φ FD ( H ω ) (cid:3)(cid:111) , (4.12)where α and β indicate space directions, L ω is the limit of L E ,ω as E → G = e h =7 . × − S is the conductance quantum. Note that the super-operator (Γ ω + L ω ) − acts on the observable appearing at its right. We now discuss the self-averaging properties of the conductivity tensor. Using thecovariant properties of the operators appearing in the Kubo formula, one finds σ αβ ( T, µ ; τ g ω ) = − πG p αα (cid:48) Tr V (cid:110) T g (cid:2) X α (cid:48) ω , H ω (cid:3) (4.13)(Γ ω + L ω ) − (cid:2) X β (cid:48) ω , Φ( H ω ) (cid:3) T † g (cid:111) p β (cid:48) β , for any g = ( p | a ) ∈ G , where repeated indices are summed over their range. Using theinvariance of the trace under conjugation, we find the simple rule of transformationˆ σ ( T, µ ; τ g ω ) = p ˆ σ ( T, µ ; ω ) p − . (4.14)As such, the conductivity tensor is invariant under the translations (1 | t ) ∈ B . Ourobservation in section 2.2 that this subgroup of G acts ergodically on Ω, become extremelyimportant because it assures us that ˆ σ ( T, µ ) is self-averaging and does not fluctuate from15ne disordered configuration to another. Indeed, given that ˆ σ ( T, µ ; ω ) = σ ( T, µ ; τ (1 | t ) ω )for any (1 | t ) ∈ B ), we can writeˆ σ ( T, µ ; ω ) = lim V → R | V ∩ B | (cid:88) t ∈ V ˆ σ ( T, µ ; τ (1 | t ) ω ) (4.15)= (cid:90) Ω d P ( ω (cid:48) ) ˆ σ ( T, µ ; ω (cid:48) ) , where the last equality follows from Birkhoff’s theorem [32]. Now, the only way toreconcile the above conclusion and (4.14), is to admit the invariance of the conductivitytensor under the point group action p − ˆ σ ( T, µ ; ω ) p = ˆ σ ( T, µ ; ω ) , ∀ p ∈ P . (4.16)The remarkable conclusion is that the invariance w.r.t. the full space group G of thenon-averaged conductivity tensor is exact even though this symmetry is broken locallyby the thermal motion of atoms. We mention that in our numerical calculations, weevaluate the isotropic part of the conductivity tensor σ ( T, µ ; ω ) =
13 3 (cid:88) α =1 σ αα ( T, µ ; ω ) , (4.17)which is manifestly invariant under the action of the entire space group.Let us stress that the above self-averaging property manifests itself only in the strictthermodynamic limit. For finite samples, there will be fluctuations w.r.t. the thermallydisordered configurations. This is because the group of symmetry transformations getsreduced when dealing with finite samples and, as a consequence, τ g ω does not explore the whole Ω when g is given all allowed values. For a finite Si crystal of cubic shape,which is built by repeating the unit cell, the rank of the group of symmetries is equalto the number N a of atoms in the crystal. Given the invariance of σ ( µ, T ; ω ) w.r.t.these transformations, when we evaluate σ ( µ, T ; ω ) for one disordered configuration, wein fact evaluate the conductivity for all τ g ω configurations. In other words, with just onecalculation, we sample N a points of Ω. Hence, if we repeat the calculation of σ ( µ, T ; ω )for a number N c of different configurations, we effectively sampled Ω at N a × N c points.Because of this amplification effect that stems from the invariance of σ ( µ, T ; ω ) relative tothe space symmetries, we expect that a good disorder average can be achieved even witha small number of disordered configurations. This is indeed observed in our simulations. There are two fundamental difficulties when attempting to evaluate (4.12) on a com-puter. The first one stems from the incompatibility between the covariant relation (3.43)for the position operator and the periodic boundary conditions. The second difficultycomes from inverting the super-operator Γ ω +L ω . Both these issues have been resolved in Up to subsets of measure zero X , A ] m,n x , x (cid:48) ( ω ) = (cid:90) R d r φ ∗ m ( r − x )([ X , A ] φ n )( r − x (cid:48) ) (4.18)= (cid:90) R d r (cid:90) R d r (cid:48) ( r − r (cid:48) ) φ ∗ m ( r − x ) A ( r , r (cid:48) ) φ n ( r (cid:48) − x (cid:48) )=( x − x (cid:48) ) (cid:90) R d r (cid:90) R d r (cid:48) φ ∗ m ( r − x ) A ( r , r (cid:48) ) φ n ( r (cid:48) − x (cid:48) )+ (cid:90) R d r (cid:90) R d r (cid:48) ( r − r (cid:48) ) φ ∗ m ( r ) A ( r , r (cid:48) ) φ n ( r (cid:48) ) . We can summarize the above calculation as[ X , A ] m,n x , x (cid:48) ( ω ) = ( x − x (cid:48) ) A m,n x , x (cid:48) ( ω ) + [ X , A ] m,n , ( ω ) . (4.19)While the right hand side makes perfect sense for an infinite samples, when the simulationproceeds over a finite crystal with periodic boundary conditions, there is an obviousproblem with the first term. In [19, 33], it was found that the optimal adaptation to theperiodic boundary conditions is through the following substitution: x − x (cid:48) → x − x (cid:48) − (cid:20)(cid:20) x − x (cid:48) L/ (cid:21)(cid:21) L, (4.20)where L is the size of the periodic super-cell of the simulation and [[ · ]] denotes theinteger part of a real number. The second term in (4.19) is a local term and there isno need for a modification when finite crystals with periodic boundary conditions areconsidered. With the proper matrix elements at hand, the finite-volume tight-bindingoperators corresponding to the commutators with the position operator are derived viathe procedure detailed in section 3.3 without any modifications. To alert the readerabout the substitution (4.20), we write the modified commutators of these tight-bindingoperators as (cid:98) X ω , A ω (cid:99) .We now focus on the super-operator Γ ω + L ω . We will only consider here the so calledrelaxation time approximation where the dissipation super-operator is proportional withidentity: Γ ω = Γ id, with Γ a positive number. Now we recall that L ω acts on operators A ω over C N ⊗ (cid:96) ( L ω ) via L ω [ A ω ] = ı [ H ω , A ω ]. Observe that, if (cid:0) (cid:15) ωa , ψ ωa (cid:1) a =1 ,...,N | L ω | (4.21)is an eigen-system for H ω , thenL ω (cid:2) | ψ ωa (cid:105)(cid:104) ψ ωb | (cid:3) = ı ( (cid:15) ωa − (cid:15) ωb ) | ψ ωa (cid:105)(cid:104) ψ ωb | . (4.22)In other words, (cid:16) (cid:15) ωa − (cid:15) ωb , | ψ ωa (cid:105)(cid:104) ψ ωb | (cid:17) a,b =1 ,...,N | L ω | (4.23)17s an eigen-system for L ω . This observation together with the fact that any operator canbe decomposed as A ω = (cid:88) a,b (cid:104) ψ ωa | A ω | ψ ωb (cid:105) | ψ ωa (cid:105)(cid:104) ψ ωb | (4.24)provide a straightforward way to invert the super-operator:(Γ id + L ω ) − [ A ω ] = (cid:88) a,b (cid:104) ψ ωa | A ω | ψ ωb (cid:105) Γ + ı ( (cid:15) ωa − (cid:15) ωb ) | ψ ωa (cid:105)(cid:104) ψ ωb | . (4.25)Finally, we can give a direct translation of the Kubo-formula (4.12) at finite-volume: σ αβ = (cid:68) − πG Vol (cid:88) a,b (cid:10) ψ ωa (cid:12)(cid:12) (cid:98) X αω , H ω (cid:99) (cid:12)(cid:12) ψ ωb (cid:11)(cid:10) ψ ωb (cid:12)(cid:12) (cid:98) X βω , Φ FD ( H ω ) (cid:99) (cid:12)(cid:12) ψ ωa (cid:11) Γ + ı ( (cid:15) ωa − (cid:15) ωb ) (cid:69) ω . (4.26)This expression is useful when the matrix elements of the Fermi operator are available.Since this quantity is not among the standard outputs of AIMD simulations, we processthis expression one step further as in [18]: σ αβ = (cid:68) − πG Vol (cid:88) a,b Φ FD ( (cid:15) ωa ) − Φ FD ( (cid:15) ωb ) (cid:15) ωa − (cid:15) ωb (4.27) × (cid:10) ψ ωa (cid:12)(cid:12) (cid:98) X αω , H ω (cid:99) (cid:12)(cid:12) ψ ωb (cid:11)(cid:10) ψ ωb (cid:12)(cid:12) (cid:98) X βω , H ω (cid:99) (cid:12)(cid:12) ψ ωa (cid:11) Γ + ı ( (cid:15) ωa − (cid:15) ωb ) (cid:69) ω . This is the expression we coded as a post-processing subroutine to the AIMD simulations.The inputs for this expression are the matrix elements of the Kohn-Sham Hamiltonians(3.35) and the overlap coefficients (3.13), as well as the xyz-coordinates of the atoms.
5. Numerical Implementation
In this section, we first present a novel electronic structure method that is onlyscaling quadratically with system size, thus facilitating second-generation Car-ParrinelloAIMD simulations of even longer length and time scales than previously thought fea-sible [7, 8]. More importantly, this approach permits to efficiently compute the exactfinite-temperature density matrix ρ ω of a given Kohn-Sham Hamiltonian H ω “on-the-fly”during the AIMD. Thereafter, the computational details of our simulatios are describedin detail. Following Alavi and coworkers [35, 36], we begin with the (Helmholtz) free energyfunctional F = Ξ + µN e + V dc , (5.1)where N e = 2 N is the number of electrons and Ξ the grand-canonical potential (GCP)for noninteracting fermions. The latter reads asΞ = − β ln det (cid:16) e β ( µS ω − H ω ) (cid:17) = − β Tr ln (cid:16) e β ( µS ω − H ω ) (cid:17) , (5.2)18ith given by β − = kT ( k = Boltzmann constant). Yet, in the low-temperature limitlim β →∞ Ξ = 2 N (cid:88) a =1 (cid:15) ωa − µN e , (5.3)the so-called band-structure energy, which is given by the sum of the lowest N doublyoccupied eigenvalues (cid:15) ωa of H ω , can be recovered andlim β →∞ F = 2 N (cid:88) a =1 (cid:15) ωa + V dc (5.4)holds. Therein, V dc accounts for double counting terms, as well as for the nuclearCoulomb interaction.In the present case of fully self-consistent KS-DFT calculations V dc [ n ω ( r )] = − (cid:90) d r (cid:90) d r (cid:48) n ω ( r ) n ω ( r (cid:48) ) | r − r (cid:48) |− (cid:90) d r n ω ( r ) δ Ξ XC [ n ω ( r )] δn ω ( r ) + Ξ XC [ n ω ( r )] + E II , (5.5)where the first term on the right hand side is the double counting correction of the Hartreeenergy, while Ξ XC [ n ω ( r )] is the finite-temperature XC grand-canonical functional and E II the nuclear Coulomb interaction. Except for the latter term, Eq. (5.5) accounts forthe difference between Ξ and the GCP for the interacting spin- Fermi gas, i.e.Ξ int [ n ω ( r )] = − β ln det (cid:16) e β ( µS ω − H ω KS ) (cid:17) − (cid:90) d r (cid:90) d r (cid:48) n ω ( r ) n ω ( r (cid:48) ) | r − r (cid:48) |− (cid:90) d r n ω ( r ) δ Ξ XC [ n ω ( r )] δn ω ( r ) + Ξ XC [ n ω ( r )] . (5.6)As before, in the low-temperature limit Ξ int [ n ω ( r )] + µN e equals to the band-structureenergy, whereas Ξ XC [ n ω ( r )] corresponds to the familiar XC energy, so that in this limit F = Ξ+ µN e + V dc = Ξ int [ n ω ( r )]+ µN e + E II is equivalent to the Harris-Foulkes functional[37, 38]. Equally than the latter, F is explicitly defined for any n ω ( r ) and obeys exactlythe same stationary point as the finite-temperature functional of Mermin [39].Whereas it is well known how to calculate V dc with linear-scaling computational effort,the computation of all occupied orbitals by diagonalization requires O ( N ) operations.Due to the fact that the band-structure term can be equivalently expressed in terms of ρ ω , the total energy can be written as E KS [ n ω ( r )] = 2 N (cid:88) a =1 (cid:15) ωi + V dc = Tr[ ρ ω H ω KS ] + V dc [ n ω ( r )] . (5.7)As a consequence, the cubic-scaling diagonalization of H ω KS can be bypassed by directlycalculating ρ ω rather than all ε i ’s. 19n order to make further progress, let us now factorize the operator of Eq. (5.2) into P terms. Given that P is even, which we shall assume in the following, Krajewski andParrinello derived the following identity1 + e β ( µS ω − H ω KS ) = P (cid:89) l =1 (cid:16) − e iπ P (2 l − e β P ( µS ω − H ω KS ) (cid:17) = P/ (cid:89) l =1 M ∗ l M l , (5.8)where the matrices M l , with l = 1 , . . . , P , are defined as M l := 1 − e iπ P (2 l − e β P ( µS ω − H ω KS ) , (5.9)while ∗ denotes complex conjugation [40]. Analog to numerical path-integral calculations[41], it is possible to exploit the fact that if P is large enough, so that the effectivetemperature β/P is small, the exponential operator e β P ( µS ω − H ω KS ) can be approximatedby a Trotter decomposition or simply by a high-temperature expansion, i.e. M l = 1 − e iπ P (2 l − (cid:16) β P ( µS ω − H ω KS ) (cid:17) + O (cid:18) P (cid:19) . (5.10)However, as we will see, here no such approximation is required, which is in contrast tothe original approach [40]. In any case, the GCP can be rewritten asΞ = − β ln P/ (cid:89) l =1 det ( M ∗ l M l ) = β P/ (cid:88) l =1 ln (det ( M ∗ l M l )) − . (5.11)As is customary in lattice gauge theory [42, p. 17], where the minus sign problemis avoided by sampling a positive definite distribution, the inverse square root of thedeterminant can be expressed as an integral over a complex field φ l , which has the samedimension M as the full Hilbert space, i.e.det ( M ∗ l M l ) − / = 1(2 π ) M (cid:90) dφ l e − φ ∗ l M ∗ l M l φ l . (5.12)Inserting Eq. (5.12) into Eq. (5.11) we end up with the following field-theoretic expressionfor the GCP: Ξ = β P/ (cid:88) l =1 ln (cid:34) π ) M (cid:90) dφ l e − φ ∗ l M ∗ l M l φ l (cid:35) = β P/ (cid:88) l =1 ln (cid:90) dφ l e − φ ∗ l M ∗ l M l φ l + const., (5.13)where φ l are appropriate vectors.All physical relevant observables can then be determined as functional derivatives ofthe GCP w.r.t. an appropriately chosen external parameter. For example, N e = − ∂ Ξ /∂µ and lim β →∞ Ξ + µN e = 2 (cid:80) N a =1 (cid:15) ωi , so that E = lim β →∞ F = 2 N (cid:88) a =1 (cid:15) ωa + V dc = ∂ ( β Ξ) ∂β − µ ∂ Ξ ∂µ + V dc . (5.14)20ince the functional derivative of the constant in Eq. (5.13) is identical to zero, all physicalinteresting quantities can be computed via ∂ Ξ ∂λ = − β P/ (cid:88) l =1 (cid:82) dφ l d (cid:80) i,j =1 ( φ l ) ∗ i (cid:16) ∂ ( M ∗ l M l ) ∂λ (cid:17) ij ( φ l ) j e − φ ∗ l M ∗ l M l φ l (cid:82) dφ l e − φ ∗ l M ∗ l M l φ l (5.15a)= − β P/ (cid:88) l =1 d (cid:88) i,j =1 (cid:18) ∂ ( M ∗ l M l ) ∂λ (cid:19) ij ( M ∗ l M l ) − ij (5.15b)= − β P/ (cid:88) l =1 Tr (cid:20) ( M ∗ l M l ) − ∂ ( M ∗ l M l ) ∂λ (cid:21) = − β P (cid:88) l =1 Tr (cid:20) M − l ∂M l ∂λ (cid:21) . (5.15c)Thereby, the left-hand side of Eq. (5.15c) holds because of Montvay and M¨unster [42,p. 18], whereas the right-hand side is due to the fact that beside being positive definite M ∗ l M l is also symmetric.Comparing Eq. (5.7) with Eq. (5.3), it is easy to see that the GCP and hence allphysical significant observables can be written as the trace of a product consisting ofthe Fermi matrix ρ ω . Specifically, Ξ = Tr[ ρ ω H ω KS ] − µN e , but because at the same time N e = Tr[ ρ ω S ω ] holds, the former can be simplified toΞ = Tr[ ρ ω ( H ω KS − µS ω )] , (5.16)where S ω = − ∂H ω KS /∂µ and ρ ω = ∂ Ξ /∂H ω KS . As a consequence, the GCP and all itsfunctional derivatives can be reduced to evaluate ρ ω based on Eq. (5.15c) with λ = H ij .Using the identity ∂M l ∂H ij = − P { ( M l − β + β ( M l − } , (5.17)for this particular case, Eq. (5.15c) eventually equals to ρ ω = ∂ Ξ ∂H ω KS = P P/ (cid:88) l =1 (cid:16) − (cid:0) M ∗ l M l (cid:1) − (cid:17) = P P (cid:88) l =1 (cid:0) − M − l (cid:1) . (5.18)In other words, the origin of the method is the notion that the density matrix, thesquare of the wavefunction at low temperature and the Maxwell-Boltzmann distributionat high temperature, can be decomposed into a sum of M − l matrices, each at highereffective temperature β/P and hence always sparser than ρ ω . Yet, contrary to the originalapproach [40], neither a Trotter decomposition nor a high-temperature expansion forEq. (5.9) has been used, so far everything is exact for any P .In particular, unlike Eq. (5.11), the determination of Ω = ∂ ( β Ω) /∂β does no longerinvolve the calculation of the inverse square root of a determinant, but just the inverseof M l , which is not only very sparse, since it obeys the same sparsity pattern as H ω KS ,but is furthermore also always better conditioned as the latter. Hence, all M − l matricesare substantially sparser than ρ ω and thus can be efficiently computed [43, 44]. In fact,for quasi one-dimensional systems, M l is tridiagonal that permits for an exact linear-scaling calculation of its inverse using a recursive scheme [45]. For all other dimensions21
00 k 600 k900 k 1200 k 1500 k
Figure 5.1:
The orbits of the atoms under the thermal motion at different temperatures. Thesimulation is for a crystal containing 216 Si atoms and each orbit is sampled at 1001 points.In these renderings, the crystal is viewed from atop of xy-plane. The units of the graphs areAngstroms. D , M l can be sought of being block-tridiagonal, where the dimensionality of each blockis d = N − (1 /D ) , eventually leading to a computational effort, which scales like N d = N − /d . Since this is only marginally better than the initial O ( N ) scaling for a generalmatrix inversion (or diagonalization), we compute M − l by solving the N e sets of linearequations M l Φ lj = ψ j , where { ψ j } is a complete set of basis functions [46]. Using apreconditioned biconjugate gradient method [47], the inverse can be exactly computedas M − l = (cid:80) N e j =1 φ lj ψ lj within O ( N ) operations. Furthermore, the formal analogy of thedecomposition to the Trotter factorization immediately suggests the possibility to applysome of the here presented ideas with benefit to numerical path-integral calculations [41].The same applies for a related area where these methods are extensively used, namelythe lattice gauge theory to quantum chromodynamics [48], whose action is rather similarto the one of Eq. (5.12). We now return to our specific simulations. Our models of crystalline silicon consistedof 216 and 1000 Si atoms in a cubic simulation box with periodic boundary conditions.For each system size, five simulations have been conducted, at T = 300 K, 600 K, 900 K,1200 K and 1500 K, respectively. All of our calculations were performed in the canonicalNVT ensemble using the second-generation Car-Parrinello AIMD method of K¨uhne and22
00 k 600 k900 k 1200 k 1500 k
Figure 5.2:
Same as Fig. 5.2 but for a simulation with 1000 Si atoms. coworkers [7, 8]. Throughout, the experimental density of crystalline silicon was assumed,which, at ambient conditions, is semiconducting and four-fold coordinated. In each run,we carefully equilibrated the system for 250 ps before accumulating statistics duringadditional 1.25 ns, resulting in a total AIMD simulation time of 15 ns.All simulations were performed at the DFT level using the mixed Gaussian and planewave code CP2K/
Quickstep [9]. In this approach, the KS orbitals are expanded in con-tracted Gaussians functions, while the electronic charge density is represented by planewaves [34]. A density cutoff of 100 Ry was employed for the latter, whereas for the formera dimer-optimized minimal basis set was used [49] of s - and p -type. As such, N = 4 in(3.9) and the linear space spanned by these wave functions is indeed a representationspace for the O (3) group. The unknown exact XC potential is substituted by the LDA[31], whereas the interactions between the valence electrons and the ionic cores are de-scribed by separable norm-conserving Goedecker-Teter-Hutter pseudopotentials [50, 51].For the sake of simplicity, the first Brillouin zone of the super cell is sampled at theΓ-point only.
6. Numerical results
In this section, we first present and analyze the output of our AIMD simulations andthen we report the output of the charge-transport post-processing subroutine detailed insection 4.4. 23 (a)(b) D spec D mob C o n f i g u r a t i o n s C o n f i g u r a t i o n s Figure 6.1:
Level statistics for a (a) 3 × × × × T = 900 K. The background displays theenergy spectra for various thermally disordered configurations. The red curve represents the variance of thelevel spacings ensembles collected at different energies. The yellow line indicates the variance of the Gaussianorthogonal ensemble. Both the red and yellow curves have their y-axis on the right side. Also shown are thespectral and the mobility gaps, as inferred from the data. In Figs. 5.1 and 5.2 we report the orbits of the atoms at different temperatures, assimulated with 216 and 1000 Si atoms, respectively. As one can see, the orbits wonderaround the equilibrium positions and the data reveals that crystal Si is quite disorderedeven at the room temperature. Let us point out that many electronic devices operate at600 K or higher under heavy loads. At these temperatures, the thermal disorder is quitepronounced. As it is well known, in such conditions, some of the wave functions canand will become affected by the phenomenon called Anderson localization [52]. Whena wave function becomes Anderson localized, its contribution to the Kubo-formula iszero. Almost as a rule [53], these localized states occur close to the edges of the energyspectrum and, for 3-dimensional crystals, it is predicted there exist mobility edges in theenergy spectrum, one in the conduction and one in valence bands, beyond which the wavefunctions remain extended. These mobility edges define the so called mobility gap andthe expectation of the charge current operator is zero when one only populates electronstates with eigen-energy within this gap. It becomes clear that the activated behaviorof the conductivity is determined by the mobility gap and not by the spectral gap. Assuch, it is extremely important to detect the mobility edges for our crystals. For this, weemploy a technique called the level statistics analysis, which has been successfully usedin the past for this very purpose [19, 54].We exemplify the process for the temperature T = 900 K, where the disorder is quitepronounced and the effects described above are more visible. Before we start, we need toexamine the spectral characteristics of the Hamiltonians. For this, we have diagonalized24 Figure 6.2: (a) Fluctuation (black dots) and average value (red line) of the chemical potential value cor-responding to the charge neutrality. (b) Same as (a) for the intrinsic charge carrier density. The data wasextracted from the spectra shown in Fig. 6.1(b) for a crystal containing 1000 Si atoms at temperature 900 K. the tight-binding KS-Hamiltonian for 1001 selected thermally disordered configurations.The result is a sequence of 1001 discrete sets of eigenvalues { (cid:15) ωa } , which we renderedon a horizontal line for each configuration and then we stuck these lines vertically. Theresulting collection of spectra then appears as the fuzzy dots seen in the background ofFig. 6.1. We recall that for covariant systems of Hamiltonians in the thermodynamic,the spectrum is in fact non-fluctuating in the sense that, if we pick any energy intervaland ask what is the probability (w.r.t. to ω ) for at least one eigenvalue to fall withinthis interval, the answer will either 0 or 100 percent. This is a consequence of the factthat T g H ω T † g have the same spectrum for all g ∈ G . In the same time, T g H ω T † g = H τ g ω and the orbit { τ g ω, g ∈ G } samples Ω entirely. The spectrum of a covariant familyof Hamiltonians is defined as the intersection of all closed subsets of the real axis thatcontain all eigenvalues with 100% probability. Rendering the spectra as in Fig. 6.1 helpsone identify this non-fluctuating spectral set.We now examine the spectra more closely. The first issue we want to address isthe appearance of the spectral gaps inside the valence and conduction bands. These areartificial features due to relative small size of the system. For a periodic system simulatedwith periodic boundary conditions ( i.e. at Γ-point) on a finite super-cell containing manyunit cells, this gaps will be explained by the coarse sampling of the Brillouin zone of theunit cell. As one can see in Fig. 6.1, many of the gaps disappear when the size of theSi crystal is increased from 3 × × × × . – m Figure 6.3:
Fluctuation (black dots) and average value (red line) of the impurity density as a function ofchemical potential. The inset shows a restricted range of the same data. The numerical values were extractedfrom the spectra shown in Fig. 6.1(b) for a crystal containing 1000 Si atoms at temperature 900 K. and it indicates that the local orbital basis is too coarse. We have verified that, indeed,increasing the local orbital basis converges the spectral gap to the standard KS-DFTvalue of 0.7 eV. We recall that the experimental value is 1.1700 eV at 4.2K [55] and thatthe experimental band gap displays a temperature dependence which has been assessedquite precisely [56, 57]. Our simulations, however, are performed with the same super-cell regardless of the temperature, hence we cannot relate them to that experimentalfact. Our conclusions based on the spectral data reported in Fig. 6.1 is that the presentsimulations are not yet precise enough for quantitative predictions. As such, we willfocus in the following only at qualitative aspects.We now turn our focus on the level statics analysis, which was performed in thefollowing way. We picked an arbitrary energy (cid:15) and, for each of the 1000 thermallydisordered configurations considered in Fig. 6.1, we identified the unique eigenvalues (cid:15) ωa and (cid:15) ωa +1 that satisfy the constraint (cid:15) ωa < (cid:15) < (cid:15) ωa +1 . Then we computed the levelspacings ∆ (cid:15) = (cid:15) ωa + j +1 − (cid:15) a + j ω , letting j take 11 consecutive values between − (cid:15) . These level spacings were subsequently normalized bytheir average.As done in [54], one can examine the histograms of these ensembles and determinewhat kind of distributions they manifest. Since the KS-Hamiltonians are real, we expectthe outcome to be either a Poisson distribution P ( s ) = e − s or a Gaussian orthogonalensemble (GOE), P GOE = π se − π s . These distributions are expected when the local-ization length of the wave functions with energy close to (cid:15) is smaller/larger than thesize of the super-cell, respectively [58]. If the super-cell is large enough, one can derivefrom these distributions the localized or de-localized character of the wave functions.Overlapped over the spectra in Fig. 6.1 is the variance (cid:104) s (cid:105) − (cid:104) s (cid:105) of the level spacingensembles collected at different energies, as well as the variance value of 0.273 computedfrom P GOE . Since the variance of the Poisson ensemble is 1, we can easily identify fromFig. 6.1 the character of the wave functions, in particular, the mobility gap. As one cansee, it extends well beyond the spectral gap.26 = 0.01 kT G =0.1 kT G = kT Figure 6.4:
Conductivity as a function of chemical potential for 216 Si atoms, T = 300 K and differentvalues of the dissipation coefficient. Shown in blue are the un-processed output for 1000 thermally disorderedconfigurations. The red curves represent the average values. The charge-neutrality point is defined by the precise value of the chemical potential µ where the charge neutrality of the crystal is achieved. Since in our calculations eachionic core carries 4 e charge, the charge neutrality condition reads:4 N a Vol = (cid:42) (cid:88) a
11 + exp( (cid:15) ωa − µ kT ) (cid:43) ω . (6.1)We want to point out that, for covariant systems, µ and the quantity inside the averagebrackets in (6.1) are self-averaging in the thermodynamic limit. However, for our finite-size crystals, these quantities will display fluctuations from one thermally disorderedconfiguration to another and the size of these fluctuations is a good indicator of howclose is the simulation to the thermodynamic limit. A rendering of the fluctuationsas well as the average value of the chemical potential µ at the neutrality point and900 K temperature are reported in Fig. 6.2(a). The data reveal an extremely low level offluctuations, characterized by a standard deviation of 0 .
063 % around the average value µ = 0 . = 0.01 kT G =0.1 kT G = kT Figure 6.5:
Same as Fig. 6.4 for T = 900 K. determined by the depletion of the valence states due to the thermal excitations: n h ( T, µ ) = (cid:42) (cid:88) (cid:15) ωa ≤ µ (cid:34) −
11 + exp( (cid:15) ωa − µkT ) (cid:35)(cid:43) ω . (6.2)The concentration of the mobile electrons is determined by the population of the con-duction states due to the thermal excitations: n e ( T, µ ) = (cid:42) (cid:88) (cid:15) ωa ≥ µ
11 + exp( (cid:15) ωa − µkT ) (cid:43) ω . (6.3)Note that at the neutrality point, we have the equality: n h ( T, µ ) = n e ( T, µ ) . (6.4)The common value of the two densities is called the intrinsic density ( n i ) of chargecarriers and it is one of the most important characteristic of Si semiconductor. Itsexperimental value at 300 K has been determined with great precision [59–61] to be n i =9 . ÷ . × cm − . Experimental data on the dependence of n i with the temperaturehas been summarized in [62, Fig. 14], from where we extracted the experimental value of1 . × cm − at T = 900 K. For our simulations, the fluctuations and the average valueof n i are reported in Fig. 6.2(b). As one can see, despite of extremely low fluctuations in µ , there are substantial fluctuations in the n i data, which reflect the extreme sensitivity28 = 0.01 kT G =0.1 kT G = kT Figure 6.6:
Direct conductivity as a function of the chemical potential for 1000 Si atoms, T = 300 K andΓ = 0 . kT (top), Γ = 0 . kT (middle) and Γ = kT (bottom). Shown in blue are the un-processed outputfor 1000 thermally disordered configurations. The red curves represent the average values. of n i on the energy spectrum. Quantitatively, the standard deviation in Fig. 6.2(b) is6 % and the average value is n i = 1 . × . This value is much lower than theexperimental value mentioned above, the main reason being the over-estimation of theband gap by our simulations.In our study, we will consider not only neutral but also Si crystals that are away fromneutrality point by letting the chemical potential µ be a variable. Experimentally, thevariation of the chemical potential can be achieved via gate potentials for thin films orvia impurity doping for bulk samples. Either way, such variations lead to changes in theelectronic structure of the crystal, which should be recomputed every time the dopinglevel is changed. Since in our simulations we use the same electronic structure, specificallythe one computed at the neutrality point, the results we present here are relevant onlyfor lightly doped or weakly gated samples where the changes in the electronic structureare expected to insignificant. To make contact with the experiment, one has to rely onthe impurity density value rather than on the chemical potential, because the former isthe parameter that can be controlled in laboratory. The impurity density is evaluatedfrom: n ( µ, T ) = | n e ( µ, T ) − n h ( µ, T ) | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:42) (cid:88) a
11 + exp( (cid:15) ωa − µkT ) (cid:43) ω − N a Vol (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (6.5)For completeness, we show in Fig. 6.3 the relation between the impurity density andchemical potential µ , as derived from the spectra shown in Fig. 6.1 and (6.5). Let us recall29 = 0.01 kT G =0.1 kT G = kT Figure 6.7:
Same as Fig. 6.6 for T = 900 K. In addition, the bottom panel contains a zoom-in. that a light to moderate doping corresponds to the experimental values n exp < cm − ,which in terms of the chemical potential means, approximately, that | µ − µ | < .
01 Ha.Let us end this section by specifying that the extrinsic hole density from Eqs. (6.2)will be used in the next section to generate the hole mobility via the relation σ = en h µ h .The hole mobility µ h will be mapped as a function of the acceptor concentration (6.5). In this section we present and analyze the numerical results on the direct conductivity(4.17). Figs. 6.4 and 6.5 report the data for a 216 Si atoms crystal at temperatures T = 300 K and T = 900 K, respectively. Similarly, Figs. 6.6 and 6.7 report the datafor a 1000 Si atoms crystal at temperatures T = 300 K and T = 900 K, respectively.The chemical potential has been varied throughout the entire energy spectral range andthe dissipation parameter Γ was sampled at three different values in these simulations.Overall, the results show a good correlation between the conductivity plots and thespectrum of the Kohn-Sham Hamiltonians. There is a significant difference between theoutputs for 216 and 1000 atom crystals, indicating that the results are not convergedyet w.r.t. the system’s size. There is also a significant dependence on the dissipationparameter Γ . It is interesting to notice that for its largest values, the fluctuations ofthe direct conductivity are drastically suppressed. As explained in [19], the convergenceto thermodynamic limit is faster for larger Γ , and this explains the suppression of thefluctuations observed in these figures. 30 = kT G = 0.1kT G = 0.01kT Figure 7.1:
Dependence with the acceptor concentration of (a) direct conductivity and (b) hole mobility,as computed for the 1000-atom Si crystal at T = 900 K. The fine lines represent unprocessed data comingfrom individual disordered configurations, hence they are a measure of the fluctuations. The thicker redlines represent the averages over 50 configurations. The simulations have been carried for three values of thedissipation parameter, which are specified in each panel. The conductivity results are in very good agreement with the mobility gap prediction,which for T = 900 K can be found in Fig. 6.1. Indeed, the conductivity is obviously notinfluenced at all by the first band of spectrum, which was determined in Fig. 6.1 to beAnderson localized. This is quite obvious in the inset of Fig. 6.7, which shows a zoominto the region around the spectral gap. Let us point out again that, on the other hand,the intrinsic and extrinsic charge carriers are highly influenced by the presence of thisAnderson localized band.In Fig. 7.1(a), we focus on the behavior of direct conductivity inside and around theinsulating gap, especially on the hole side. We chose to investigate only at the crystal with1000 Si atoms and T = 900 K because it is the most converged system. At T = 300 K,the conductivity curves display a pronounced dependence on the spectral details whichare not yet converged, hence the analysis will no be reliable. As it is customary, thetransport coefficient has been plotted as a function of the acceptor concentration ratherthan chemical potential. The behavior of σ seen in Fig. 7.1(a) is as expected. Deep insidethe insulating gap, the direct conductivity saturate at a value proportional to Γ and, asthe chemical potential moves towards the valence band, an activated behavior takes over.In Fig. 7.1(b), we report the hole mobility as a function of the acceptor concentration.The functional shape is in good agreement with the measured one (see [55, Fig. 21.8]).Let us point out that, when the chemical potential is inside a spectral gap, the mobilityis proportional with the dissipation Γ , hence with inverse of the relaxation time. Incontrast, when the chemical potential is inside a spectral band, as is the case of a metal,the mobility is proportional with relaxation time, hence inverse proportionally with Γ .As such, one should not be surprised by the behavior with Γ seen in Fig. 7.1(b).31 . Conclusions We have derived disordered tight-binding models based on AIMD outputs and for-mulated a Kubo-formalism that preserves the self-averaging property of the transportcoefficients. The Kubo-formalism was coded as a post-processing subroutine to a stan-dard AIMD code and preliminary results on the transport coefficients of crystals Si wereobtained at various temperatures. According to our study, the thermal disorder can havemeasurable effects even at room temperature.
8. Acknowledgements
Emil Prodan acknowledges financial support from USA National Science Founda-tion through grant DMR-1823800. Part of this project has received funding from theEuropean Research Council (ERC) under the European Union’s Horizon 2020 researchand innovation programme (grant agreement No 716142). The authors would like tothank the Paderborn Center for Parallel Computing (PC ) for the generous allocation ofcomputing time on OCuLUS and the FPGA-based supercomputer NOCTUA. References [1] E. Abrahams, P.W. Anderson, D. Licciardello, T. Ramakrishnan,
Scaling theory of localization:Absence of quantum diffusion in two dimensions , Phys. Rev. Lett. 42, 673-676 (1979).[2] J. Bellissard,
Noncommutative Geometry of Aperiodic Solids , in H. Ocampo, E. Pariguan, S. Pay-cha, editors,
Geometric and Topological Methods for Quantum Field Theory (World Sci. Publ.,River Edge, NJ, 2003).[3] J. Bellissard,
Delone sets and material science: a program , Progress in Mathematics , 403426(2015).[4] D. W. Fitting, W. P. Dub´e, T. A. Siewert,
Monitoring the solidification of single-crystal castingsusing high-energy X-ray diffraction , Journal of Metals 51, No. 7 (1999).[5] J. Song, E. Prodan,
Quantization of topological invariants under symmetry-breaking disorder , Phys.Rev. B , 195119 (2015).[6] T. D. K¨uhne, E. Prodan, Disordered crystals from first principles I: Quantifying the configurationspace , Annals of Physics , 120-149 (2018).[7] T. D. K¨uhne, M. Krack, F. R. Mohamed, M. Parrinello,
Efficient and Accurate Car-Parrinello-likeApproach to Born-Oppenheimer Molecular Dynamics , Phys. Rev. Lett. , 066401 (2007).[8] T. D. K¨uhne, Second generation Car-Parrinello molecular dynamics , WIREs Comput. Mol. Sci. ,391-406 (2014).[9] T. D. K¨uhne et al., CP2K: An Electronic Structure and Molecular Dynamics Software Package -Quickstep: Efficient and Accurate Electronic Structure Calculations , J. Chem. Phys. , 194103(2020).[10] J. Bellissard, A. van Elst, H. Schulz-Baldes,
The non-commutative geometry of the quantum Halleffect , J. Math. Phys. , 5373-5451 (1994).[11] H. Schulz-Baldes, J. Bellissard, A kinetic theory for quantum transport in aperiodic media , J. Stat.Phys. , 991-1026 (1998).[12] H. Schulz-Baldes, J. Bellissard, Anomalous transport: A mathematical framework , Rev. Math.Phys. , 1-46 (1998).[13] G. Androulakis, J. Bellissard, C. Sadel, Dissipative dynamics in semiconductors at low tem- pera-ture , J. Stat. Phys. , 448-486 (2012).[14] E. Prodan,
Quantum transport in disordered systems under magnetic fields: A study based onoperator algebras , Applied Mathematics Research eXpress Vol. 2013, 176255 (2013).[15] Y. Xue, E. Prodan,
The noncommutative Kubo-formula: Applications to Transport in DisorderedTopological Insulators with and without Magnetic Fields , Phys. Rev. B , 155445 (2012).[16] Y. Xue, E. Prodan, Quantum criticality at the Chern-to-normal insulator transition , Phys. Rev. B , 115141 (2013).
17] J. Song, E. Prodan,
Characterization of the quantized Hall insulator phase in the quantum criticalregime , Euro. Phys. Lett. , 37001 (2014).[18] E. Prodan, J. Bellissard,
Mapping the current-current correlation function near a quantum criticalpoint , Annals of Physics , 1-15 (2016).[19] E. Prodan,
A Computational Non-Commutative Geometry Program for Disordered TopologicalInsulators , (Springer, Berlin, 2017).[20] E. Cances, P. Cazeaux, M. Luskin,
Generalized Kubo formulas for the transport properties ofincommensurate 2D atomic heterostructures , J. Math. Phys. , 063502 (2017).[21] D. Massatt, S. Carr, M. Luskin, C. Ortner, Incommensurate heterostructures in momentum space ,Multiscale Model Simul. , 429-451 (2018).[22] D. Massatt, S. Carr, M. Luskin, Ecient computation of Kubo conductivity for incommensurate 2Dheterostructures , Eur. Phys. J. B , 60 (2020).[23] J. Singh, Physics of semiconductors and their heterostructures , (McGraw-Hill, New York, 1993).[24] H. Wondratschek, U. M¨uller, editors,
International tables for crystallography Vol. A1: Symmetryrelations between space groups , (Kluwer Acad. Publ., Dordrecht, 2004)[25] D. Ruelle,
Statistical mechanics: rigorous results , (W. A. Benjamin, London, 1969).[26] W. Kohn, L. J. Sham,
Self-consistent equations including exchange and correlation effects , Phys.Rev. , A1133 (1965).[27] E. Prodan, P. Nordlander,
On the Kohn-Sham equations with periodic background potentials , J.Stat. Phys. , 967-992 (2003).[28] E. Canc´es, S. Lahbabi, and M. Lewin,
Mean-eld electronic structure models for disordered materials ,in Proceeding of the International Congress on Mathematical Physics, Aalborg (Denmark), (2012).[29] S. Lahbabi,
Mathematical study of quantum and classical models for random materials in the atomicscale , Mathematical Physics [math-ph]. Universit de Cergy Pontoise, (2013).[30] W. Kohn, L. J. Sham,
Self-consistent equations including exchange and correlation effects , Phys.Rev. , A1133 (1965).[31] D. M. Ceperley, B. J. Alder,
Ground State of the Electron Gas by a Stochastic Method , Phys. Rev.Lett. , 566-569 (1980).[32] G. D. Birkhoff, Proof of the ergodic theorem , Proc. Natl. Acad. Sci. USA , 656-660 (1931).[33] C. Bourne, E. Prodan, Non-Commutative Chern Numbers for Generic Aperiodic Discrete Systems ,J. Phys. A: Math. & Theor. , 235202 (2018).[34] G. Lippert, J. Hutter, M. Parrinello, A Hybrid Gaussian Plane Wave Density Functional Scheme ,Mol. Phys. , 477-487 (1997).[35] A. Alavi, D. Frenkel, Grand-canonical simulations of solvated ideal fermions. Evidence for phaseseparation , J. Chem. Phys. , 9249 (1992).[36] A. Alavi, J. Kohanoff, M. Parrinello, D. Frenkel, Ab Initio Molecular Dynamics with ExcitedElectrons , Phys. Rev. Lett. , 2599 (1994).[37] J. Harris, Simplified method for calculating the energy of weakly interacting fragments , Phys. Rev.B , 1770 (1985).[38] W. M. C. Foulkes, R. Haydock, Tight-binding models and density-functional theory , Phys. Rev. B , 12520 (1989).[39] N. D. Mermin, Thermal Properties of the Inhomogeneous Electron Gas , Phys. Rev. 137, A1441(1965).[40] F. R. Krajewski, M. Parrinello,
Stochastic linear scaling for metals and nonmetals , Phys. Rev. B , 233105 (2005).[41] C. John, T. Spura, S. Habershon, T. D. K¨uhne, Quantum ring-polymer contraction method: In-cluding nuclear quantum effects at no additional computational cost in comparison to ab initiomolecular dynamics , Phys. Rev. E , 043305 (2016).[42] I. Montvay, G. M¨unster, Quantum Fields on a Lattice , (Cambridge Monographs on MathematicalPhysics, Cambridge, 1994).[43] M. Ceriotti, T. D. K¨uhne, M. Parrinello,
An efficient and accurate decomposition of the Fermioperator , J. Chem. Phys. , 024707 (2008)[44] M. Ceriotti, T. D. K¨uhne, M. Parrinello,
A hybrid approach to Fermi operator expansion , AIPConf. Proc. , 658-661 (2009).[45] E. M. Godfrin,
A method to compute the inverse of an n-block tridiagonal quasi-Hermitian matrix ,J. Phys.: Condens. Matter , 7843 (1991).[46] D. Richters, T. D. K¨uhne, Self-consistent field theory based molecular dynamics with linear system-size scaling , J. Chem. Phys. , 134109 (2014).[47] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery,
Numerical Recipes: The Art of cientific Computing , (Cambridge University Press, Cambridge, 2007).[48] J. B. Kogut, The lattice gauge theory approach to quantum chromodynamics , Rev. Mod. Phys. ,775 (1983).[49] J. VandeVondele, J. Hutter, Gaussian basis sets for accurate calculations on molecular systems ingas and condensed phases , J. Chem. Phys. , 114105 (2007).[50] S. Goedecker, M. Teter, J. Hutter,
Separable Dual-Space Gaussian Pseudopotentials , Phys. Rev. B , 1703-1710 (1996).[51] M. Krack, Pseudopotentials for H to Kr Optimized for Gradient-Corrected Exchange-CorrelationFunctionals , Theor. Chem. Acc. , 145-152 (2005).[52] E. Abrahams, P. W. Anderson, D. Licciardello, T. Ramakrishnan,
Scaling theory of localization:absence of quantum diffusion in two dimensions , Phys. Rev. Lett. , 673 (1979).[53] M. Aizenman, Localization at weak disorder: some elementary bounds , Rev. Math. Phys. , 1163(1994).[54] E. Prodan, Disordered topological insulators: A non-commutative geometry perspective , J. Phys.A: Math & Theor , 113001 (2011).[55] S. Kasap, P. Capper, eds., Springer handbook of electronic and photonic materials , (Springer,Berlin, 2006).[56] Y. P. Varshni,
Temperature dependence of the energy gap in semiconductors , Physica , 149(1967).[57] K. P. O’Donnell, X. Chen, Temperature dependence of semiconductor band gaps , Appl. Phys. Lett. , 2924 (1991).[58] K. B. Efetov, Supersymmetry in disorder and chaos, (Cambridge University Press, Cambridge, UK,1997).[59] A. B. Sproul, M. A. Green, Improved value for the silicon intrinsic carrier concentration from 275to 375 K , J. Appl. Phys. , 846 (1991).[60] K. Misiakos, D. Tsamakis, Accurate measurements of the silicon intrinsic carrier density from 78to 340 K , J. Appl. Phys. , 3293 (1993).[61] P. P. Altermatt, A. Schenk, F. Geelhaar, G. Heiser, Reassessment of the intrinsic carrier density incrystalline silicon in view of band-gap narrowing , J. Appl. Phys. , 1598 (2003).[62] C. D. Thurmond, The standard thermodynamic functions for the formation of electrons and holesin Ge, Si, GaAs , and GaP, Journal of The Electrochemical Society , 1133 (1975)., 1133 (1975).