Downfolding approaches to electron-ion coupling: Constrained density-functional perturbation theory for molecules
DDownfolding approaches to electron-ion coupling: Constrained density-functionalperturbation theory for molecules
Erik G. C. P. van Loon, Jan Berges, ∗ and Tim O. Wehling Institut f¨ur Theoretische Physik, Universit¨at Bremen, Otto-Hahn-Allee 1, D-28359 Bremen, Germany andBremen Center for Computational Materials Science, Universit¨at Bremen, Am Fallturm 1a, D-28359 Bremen, Germany
Constrained electronic-structure theories enable the construction of effective low-energy modelsconsisting of partially dressed particles. Since the interpretation and physical content of these theo-ries is not straightforward, we investigate the minimal example of partially dressed phonons in smallmolecules to explore the properties of constrained density-functional perturbation theory (cDFPT).We show that the dipole selection rules determine whether the partially dressed phonons satisfyGoldstone’s theorem and we proof that electronic screening always lowers the phonon frequencies.We illustrate the theory with calculations for nitrogen and benzene.
I. INTRODUCTION
Electrons and nuclei together determine the micro-scopic properties of materials and molecules. The quan-titative understanding of this interplay is a formidabletask, since it is a many-body problem involving a largenumber of quantum particles.
Ab initio -derived low-energy models are a way to ad-dress this problem [1]. The electronic structure is dividedinto high- and low-energy states. The high-energy statesare integrated out, leaving an effective low-energy model.The “bare” particles that enter the low-energy model arein fact the partially dressed particles of the full system.These low-energy degrees of freedom can subsequentlybe analyzed in more detail. For example, a detailed andcomputationally expensive treatment of electronic [2] andelectron-ion [3, 4] correlations as well as complex dynam-ical and non-equilibrium phenomena is often only feasiblefor the low-energy model.
Ab initio -based low-energy models also form the basisfor including environmental effects like screening [5] andhybridization [6]. This two-step approach has the big ad-vantage that the changes in the environment only enterthe second step of the evaluation, which can be substan-tially cheaper to evaluate than the full calculation.One way for establishing low-energy models based on ab initio calculations are the so-called “constrained”methods. In these methods, the low-energy degrees offreedom are frozen, so that the effective interaction isscreened only by processes involving high-energy elec-trons. In this way, the constrained density-functionalperturbation theory (cDFPT) [4] creates a low-energymodel consisting of partially screened phonons, low-energy electrons, and an electron-phonon interaction.These three quantities can all be extracted from ab ini-tio calculations. Similar constrained theories exist forthe electron-electron interaction [7], in particular theconstrained random-phase approximation [8]. Together, ∗ Currently affiliated with the U Bremen Excellence Chair “Newmaterials on demand” at MAPEX Center for Materials and Pro-cesses, Universit¨at Bremen. these approaches have allowed for the investigation of thecombined effect of electron-electron and electron-phononinteractions, e.g., in the fullerides [4, 9–13]. The cDFPThas been applied to several materials with electron-phonon coupling [6, 14–16].Here, our aim is to improve the understanding of thegeneral structure and properties of downfolding theoriesfor electron-ion problems in general and of cDFPT inparticular. We show that Goldstone’s theorem does notgenerally apply to partially dressed phonons but thatsymmetry-based selection rules allow us to construct elec-tronic target spaces that satisfy Goldstone’s theorem.We also show that the electronic screening reduces thephonon frequencies and consider the basis transforma-tion between bare and dressed phonons, i.e., harmonicmode-mode coupling.To illustrate these findings, we have calculated theelectron-ion coupling in small molecules. The equationsof cDFPT are matrix relations in terms of both electronicand vibronic/phononic modes. In a crystalline solid, thismeans that all objects carry momentum labels in addi-tion to their mode label and accurate calculations requirea dense momentum grid. On the other hand, in moleculesthere is no momentum, the Hilbert space is finite, and allcalculations are substantially easier. This allows us toelucidate important aspects of cDFPT in unprecedenteddetail.The set-up of this paper is as follows: First, we con-struct a general framework for the calculation of partiallyand fully dressed phonon properties within the Born-Oppenheimer approximation. We then show how cDFPTfits in this general framework and proof several propertiesof cDFPT. Subsequently, we illustrate the theory withnumerical calculations for nitrogen and benzene.
II. METHOD
Let us start with some remarks on terminology: Wewill generally use the term phonons for the ionic displace-ment eigenmodes, even for molecules, where one mightalso call them vibrons. We will be calculating the clas-sical energy landscape corresponding to these displace- a r X i v : . [ c ond - m a t . s t r- e l ] F e b ments. We use the term ions to denote the nuclei andany core electrons that are fixed to the nuclei in theelectronic-structure calculations. We set (cid:126) = 1 and mea-sure all energies and frequencies in eV.We follow a variational approach to the electron-ioncoupling [17–19] and the screening of phonons in general,and we only specify density-functional theory (DFT) atthe end. We show that the relations between DFPT andcDFPT are particularly clear in this variational descrip-tion. A. Ab initio electronic structure
The starting point for our analysis is the Born-Oppenheimer approximation: ionic and electronic de-grees of freedom are formally separated. The electroniccoordinates are supposed to be much faster, so that wecan assume that they are always relaxed into an instan-taneous ground state corresponding to a specific ionicconfiguration.For a system consisting of N ions, there are 3 N ioniccoordinates R µ . The electronic degrees of freedom arewritten as ψ ∈ Ω where Ω is a sufficiently smooth man-ifold. Ω could consist, e.g., of wave functions or densitymatrices, depending on the electronic-structure methodthat is used. We assume that there is an energy func-tional E ( ψ ; R µ ) and that this functional is sufficientlysmooth to calculate all necessary derivatives. The en-ergy E ( R µ ) corresponding to a fixed set of ionic coordi-nates R µ is found as the minimum of the functional E with respect to ψ : E ( R µ ) = min ψ ∈ Ω E ( ψ ; R µ ). The elec-tronic configuration minimizing [20] the energy ψ ∗ ( R µ ),implicitly defined by E ( R µ ) = E ( ψ ∗ ; R µ ), depends on R µ .Therefore, a variation of R µ has two effects on the energy E ( R µ ): explicitly, and implicitly via ψ ∗ ( R µ ). The latterdescribes the electronic screening of ionic displacementsand is the main purpose of our investigations. B. Force
The first derivative of the energy is the force F . Thisis a 3 N -vector, containing the force on every ion as 3-vectors. The force corresponding to a specific ionic con-figuration R is − F µ ( R ) = dEdR µ (cid:12)(cid:12)(cid:12) R (1)= ∂ E ∂R µ + ∂ E ∂ψ ∂ψ ∗ ∂R µ (2)= ∂ E ∂R µ + 0 · ∂ψ ∗ ∂R µ . (3)The last line follows since ψ ∗ ( R ) is the minimum of E .This equation shows that the force can be obtained from ψ ∗ ( R ) without having to take into account changes in the electronic configuration. This is a manifestation ofthe 2 n + 1 theorem in electronic-structure theory [21–23].For establishing low-energy models for a given ionicconfiguration R , we will use constrained theories. Theyrestrict the electronic configuration space to some sub-space Ω (cid:48) ( R ) ⊂ Ω, with ψ ∗ ( R ) ∈ Ω (cid:48) . This defines anew energy E (cid:48) ( R µ ) = min ψ (cid:48) ∈ Ω (cid:48) E ( ψ (cid:48) ; R µ ) and a ψ (cid:48)∗ ( R µ )with E ( ψ (cid:48)∗ ( R µ ); R µ ) ≡ E (cid:48) ( R µ ). Since the constrainedvariational space is smaller, the following relations hold: E ( R µ ) ≤ E (cid:48) ( R µ ) , (4) E ( R ) = E (cid:48) ( R ) , (5) − F (cid:48) µ ( R ) = dE (cid:48) dR µ (cid:12)(cid:12)(cid:12) R (6)= ∂ E ∂R µ + 0 · ∂ψ (cid:48)∗ ∂R µ (7)= − F µ ( R ) . (8)The constrained theory gives the same forces at R , sincethe change in ψ does not enter the equation for the force. C. Phonons
Important information about the ionic degrees of free-dom is contained in the second derivative of the energy,a 3 N × N matrix. The eigenmodes of the dynamicalmatrix ˆ ω , with √ m µ m ν ˆ ω µν = d E/dR µ dR ν , where m µ are atomic masses, are called phonons and the associatedeigenvalues are the phonon frequencies/energies. Unlikethe forces, the dynamical matrix is different in the con-strained theory,ˆ ω µν = 1 √ m µ m ν d EdR µ dR ν , (9)ˆ ω (cid:48) µν = 1 √ m µ m ν d E (cid:48) dR µ dR ν , (10)ˆ ω ≤ ˆ ω (cid:48) . (11)The last inequality follows from Eqs. (4), (5), and (8).The inequality should be understood in the usual wayfor (symmetric) matrices: the difference ∆ˆ ω = ˆ ω (cid:48) − ˆ ω is a positive-definite matrix.Equation (11) shows that the phonons in the con-strained theory will always have a higher energy thanin the full theory. The constraints prevent the electronsfrom completely moving along with the ions, thus increas-ing the energy cost of ionic displacements. D. Goldstone’s theorem
Goldstone’s theorem states that every spontaneouslybroken continuous symmetry creates a massless ( ω = 0)bosonic excitation. In electronic-structure theory, theions break the three continuous translation symmetries,creating three acoustic phonon modes. In addition,molecules can have spontaneously broken rotation sym-metries and corresponding massless modes.Goldstone’s theorem for phonons easily follows fromour construction of the dynamical matrix. Assume thatthere is a continuous symmetry T λ parametrized by a realnumber λ . The total energy is invariant under this sym-metry, E ( T λ R µ ) = E ( R µ ), and it is therefore possible toconstruct a representation of T λ acting on the ψ ∈ Ω,with E ( T λ ψ ; T λ R µ ) = E ( ψ ; R µ ). Concretely, a transla-tion that acts on both the ionic coordinates R µ and theelectronic configuration ψ will leave the total energy un-changed.Since T λ is a continuous symmetry, we can take λ smalland write T λ R µ ≈ R µ + λ ∆ R µ . The translation definesa direction in displacement space and the energy is con-stant in this direction, E ( R µ ) = E ( R µ + λ ∆ R µ ). Thisdisplacement is thus an eigenmode of the dynamical ma-trix, with eigenvalue (energy) zero.This argument does not transfer to the constrainedtheory [4]. The constraints can break the continuoussymmetry explicitly: T λ ψ / ∈ Ω (cid:48) ( R ). In that case, theelectronic configuration cannot completely move alongwith the translation symmetry due to the constraints and E (cid:48) ( T λ R ) (cid:54) = E (cid:48) ( R ) = E ( R ) = E ( T λ R ). Due to thislast equality, we conclude E (cid:48) ( T λ R ) > E (cid:48) ( R ). How-ever, this inequality is not yet sufficient to conclude thatthe Goldstone modes will always acquire a finite energy in the constrained theory: it is still possible that con-strained and full theory agree to order δR and only dif-fer at order δR (or higher). In that case, the partiallydressed phonons will still satisfy Goldstone’s theorem.For cDFPT, we show in Sec. II G that dipole selectionrules can be used to construct a Ω (cid:48) that guarantees thatthe uniform translation modes stay massless. E. Bare and dressed frequencies
This brings us to the relation between the bare anddressed dynamical matrices. Both are second deriva-tives of the energy with respect to the displacement, withthe subtlety that the electronic configuration adjusts tothe displacement. This electronic response is where thedifference between the full and constrained theory (i.e.,DFPT and cDFPT) originates from. It is useful to con-sider this point in detail. We measure the atomic dis-placements δR with respect to some initial ionic posi-tions R with electronic configuration ψ = ψ ∗ ( R ). Todetermine the second derivative, it is sufficient to Taylorexpand the energy functional to second order in δR and δψ . Summation over repeated indices is implied and weuse Latin indices a, b for the electronic degrees of free-dom [24]. E ( ψ ; R + δR ) − E ( ψ ; R ) = ∂ E ∂R µ δR µ + 12 δψ a ∂ E ∂ψ a ∂ψ b δψ b + 12 δR µ ∂ E ∂R µ ∂R ν δR ν + δψ a ∂ E ∂ψ a ∂R µ δR µ . (12)There is no first-order contribution in δψ , as we saw whencalculating the force. For a given displacement δR , thereis a physical (minimal energy) solution ψ ∗ ( δR ) definedby ∂ E ( ψ ∗ ( δR ); R µ ) /∂ψ a = 0. At small δR µ , δψ ∗ b ( δR ) islinear, with0 = ∂ E ∂ψ a ∂ψ b δψ b + ∂ E ∂ψ a ∂R µ δR µ , (13) δψ ∗ b ( δR ) = − (cid:18) ∂ E ∂ψ a ∂ψ b (cid:19) − (cid:18) ∂ E ∂ψ a ∂R µ (cid:19) δR µ (14) ≡ −E − ψ a ψ b E ψ a R µ δR µ . (15)Here, we have introduced the subscript partial-derivativenotation E x = ∂ E /∂x , and the inverse is a matrix inver-sion. We reinsert this result into the energy functional,Eq. (12), E ( ψ ∗ ( δR ); R + δR ) − E ( ψ ; R ) = E R µ δR µ + 12 δR µ (cid:104) E R µ R ν − E R µ ψ a E − ψ a ψ b E ψ b R ν (cid:105) δR ν . (16) The term in brackets determines the interatomic forceconstants, √ m µ m ν ˆ ω µν = d EdR µ dR ν (17)= d E ( ψ ∗ ( R ); R ) dR µ dR ν (18)= E R µ R ν − E R µ ψ a E − ψ a ψ b E R ν ψ b . (19)The first term here is the bare energy cost of displacingions with fixed electronic configuration, the second termrepresents the screening by the electrons. Note that thesecond derivative E ψ a ψ b is the Hessian on the electronicmanifold and summation over the electronic degrees offreedom is implied. Since ψ ∗ is a local minimum of theenergy functional, the electronic Hessian E ψψ is positivedefinite and the second term in Eq. (19) can be inter-preted as an inner product (cid:104)E R µ ψ , E R ν ψ (cid:105) E − ψψ . This showsthat the µ = ν part of this term is always negative (dueto the − reduces the eigenvalues of ˆ ω . Screening reduces thephonon energies.The relation between dynamical matrices of the fulland the constrained theory is particularly clear inEq. (19). The second term contains an implicit inter-nal sum over the dimensions of the electronic manifold,and in the constrained theory this sum is restricted tothe submanifold. It can be useful to analyze this sumterm by term, i.e., to perform so-called fluctuation diag-nostics [16, 25], F. Phonons and density-functional theory
Many flavors of ab initio calculations exist, specified bytheir electronic coordinate space Ω and energy functional E . The relations and proofs given so far do not dependon the precise choice of E – as long as a single variationalfunctional is used consistently – although there will ofcourse be quantitative differences. In the remainder ofthis paper, we restrict ourselves to DFT, but one couldalso apply these considerations to, e.g., Hartree-Fock orvariational Monte Carlo.The determination of phonons and the electron-phononcoupling from DFT is done via density-functional pertur-bation theory (DFPT). Here, we will only state the mostrelevant formulas; more detailed derivations are foundelsewhere [4, 26].In the framework of DFT, the energy is a functional E ( ρ ; R ) of the electronic density ρ ( r ) and ionic coordi-nates R . The electronic properties are most easily ex-pressed using the Kohn-Sham basis {| m (cid:105)} . In this basis,the electronic density is represented by a density matrix ρ . In other words, using m, n to denote the electronic lev-els, the connection to the previous sections is ψ a ≡ ρ mn .Every degrees of freedom a in the electronic manifoldcorresponds to an electronic transition ( mn ) and whenthere are N KS Kohn-Sham states, the electronic manifoldΩ consists of N KS × N KS matrices and the number of elec-tronic degrees of freedom N el satisfies N el ≤ N KS × N KS ,where the inequality follows since not every matrix isalso a valid density matrix. See Appendix A for a min-imal example. The matrix structure is necessary sinceionic displacements generally lead not just to changesin Kohn-Sham energies but also to the hybridization ofKohn-Sham orbitals.The interatomic force constants can be expressed interms of the electronic ground-state density ρ ( r ; R ) forgiven R as d E ( ρ ; R ) dR µ dR ν = (cid:90) d r ∂V ext ( r ; R ) ∂R µ ∂ρ ( r ; R ) ∂R ν + (cid:90) d r ∂ V ext ( r ; R ) ∂R µ ∂R ν ρ ( r ; R ) + E R µ R ν , (20)where V ext is the bare external potential of an individ-ual electron amid the ensemble of ions and E R µ R ν is thesecond derivative of the ionic repulsion.The change of the electronic density follows from the Kohn-Sham equations and reads ∂ρ ( r ; R ) = 2 (cid:88) m,n ≤ N KS f ( (cid:15) m ) − f ( (cid:15) n ) (cid:15) m − (cid:15) n (cid:104) n | r (cid:105) (cid:104) r | m (cid:105) (cid:104) m | ∂ ˆ V eff ( R ) | n (cid:105) , (21)where V eff is the dressed self-consistent potential and | m (cid:105) , | n (cid:105) are Kohn-Sham eigenstates. Since ∂V eff depends on ∂ρ , this equation must be solved self-consistently. Withthat, we find the electronic contribution to the inter-atomic force constants (cid:90) d r ∂V ext ( r ; R ) ∂R µ ∂ρ ( r ; R ) ∂R ν = 2 (cid:88) m,n ≤ N KS f ( (cid:15) m ) − f ( (cid:15) n ) (cid:15) m − (cid:15) n (cid:104) n | ∂ ˆ V ext ( R ) ∂R µ | m (cid:105) (cid:104) m | ∂ ˆ V eff ( R ) ∂R ν | n (cid:105) . (22)This is nothing but the bare electronic susceptibility χ weighted with the bare and dressed (DFPT) deformation-potential matrix elements d µnm = (cid:104) n | ∂ ˆ V ext /∂R µ | m (cid:105) and˜ d νmn = (cid:104) m | ∂ ˆ V eff /∂R ν | n (cid:105) , respectively. The deformationpotential is related to the electron-phonon coupling asappears in a Hamiltonian via g = d/ √ ω (cid:48) m .The cDFPT is obtained by picking a target subspace N (cid:48) KS and restricting the summations in Eqs. (21) and (22)by excluding summands where both n and m lie in thetarget space. This yields the partially screened phononfrequencies. The relation between the dynamical matri-ces in DFPT ( ˆ ω ) and cDFPT (ˆ ω (cid:48) ) is [16] given by √ m µ m ν ˆ ω µν = √ m µ m ν ˆ ω (cid:48) µν + 2 (cid:88) m,n ∈ target d ∗ µmn f ( (cid:15) m ) − f ( (cid:15) n ) (cid:15) m − (cid:15) n ˜ d νmn . (23)The cDFPT formulation presented here seems to de-viate from the general result Eq. 19. In the latter, bothends of the self-energy diagram have the same electron-phonon vertex, which presently would read E R µ ρ mn .These are connected by a susceptibility (cid:15) − ρ ab ρ cd , which isa ( N el × N el ) × ( N el × N el ) matrix. Here, the full suscepti-bility is replaced by the bare susceptibility and the elec-tronic interactions are absorbed into one of the vertices,as in Fig. 1. The advantage of this approach is that thebare susceptibility is a diagonal ( N el × N el ) × ( N el × N el )matrix and it is sufficient to only consider the N el × N el diagonal elements. The figures in this manuscript onlyshow these N el × N el matrix elements. G. Saving Goldstone’s theorem via selection rules
The partially dressed cDFPT phonons are not guaran-teed to satisfy Goldstone’s theorem. However, here weshow that suitably chosen target spaces will guaranteemassless phonons corresponding to uniform translations.First, let us consider systems with an inversion symme-try, i.e., invariance under x (cid:55)→ − x . For uniform transla-tions, ∂V /∂R is odd under inversion. Thus, if all target-space orbitals are either even ( gerade ) or odd ( ungerade )under inversion, then the displacement-potential matrixelement d = (cid:104) m | ∂V /∂R | n (cid:105) is antisymmetric in total andtherefore zero. In that case, there is no contribution tothe phonon frequency from target-space electronic tran-sitions and the partially dressed phonon has the samefrequency as the fully screened mode, which is massless.In fact, an even stronger dipole selection rule appliesto uniform-translation phonons. If a uniform transla-tion λ is performed on the nuclear coordinates, then theKohn-Sham eigenfunctions are transformed by applyingexp( iP · λ ), since the total momentum operator P is thegenerator of uniform translations, | m λ (cid:105) = e iP λ | m (cid:105) , (24)and the Kohn-Sham eigenvalues are invariant under thistransformation by symmetry. The initial density matrixwas ρ = (cid:88) m f ( (cid:15) m ) | m (cid:105) (cid:104) m | , (25)and the density matrix after translation is ρ λ = (cid:88) m f ( (cid:15) m ) | m λ (cid:105) (cid:104) m λ | . (26)We expand up to linear order in λ , ρ λ = (cid:88) m f ( (cid:15) m )(1 + iP λ ) | m (cid:105) (cid:104) m | (1 − iP λ ) ,ρ λ − ρ = (cid:88) m f ( (cid:15) m ) (cid:0) iP λ | m (cid:105) (cid:104) m | − | m (cid:105) (cid:104) m | iP λ (cid:1) , ( ρ λ − ρ ) ab = iλ (cid:0) f ( (cid:15) b ) − f ( (cid:15) a ) (cid:1) (cid:104) a | P | b (cid:105) . (27)In the general analysis of Goldstone’s theorem, westated that the partially dressed phonon satisfies Gold-stone’s theorem if ρ λ = T λ ρ ∈ Ω (cid:48) for small λ , where Ω (cid:48) is the cDFPT submanifold consisting of density matriceswith the target space subblock frozen to the initial value.In other words, ρ λ is part of the cDFPT submanifold aslong as ( ρ λ − ρ ) ab = 0 for all target states | a (cid:105) , | b (cid:105) . Theprevious derivation shows that this is the case for uni-form translation modes if (cid:104) a | P | b (cid:105) is zero for any pair oftarget states | a (cid:105) , | b (cid:105) , i.e., the Goldstone modes are pre-served if there are no long-wavelength dipole-allowed [27]transitions possible within the target space. III. POSITIVITY OF FLUCTUATIONDIAGNOSTICS
As we have seen in Eq. (11), dressed phonons havea lower energy than bare phonons, since electrons canmove along to screen the ionic charges. We have also > < ˜ dµ k k > < d > < ˆ χ ˆ U >< ˜ d k k k k = +FIG. 1. Relation between dressed and bare electron-phononcoupling. stated that ∆ˆ ω can be interpreted mathematically asan inner product of electron-phonon vertices, with themetric given by the electronic susceptibility. Here, we willmake this proof more explicit for cDFPT. This analysisshows that every contribution ∆ˆ ω µµ,mn in the fluctuationdiagnostics is positive. This is helpful, since it meansthat no cancellations occur and the relative contributionof specific fluctuations can be quantified easily.The proof is based on the relation [4, 28, 29] between˜ d and d , also illustrated in Fig. 1,˜ d µk k = d µk k (cid:20) II − ˆ χ ˆ U (cid:21) k k ,k k , (28)where k i labels the electronic states, (cid:2) ˆ χ (cid:3) k k ,k k = δ k k δ k k f ( (cid:15) k ) − f ( (cid:15) k ) (cid:15) k − (cid:15) k (29)is the Lindhard bubble in the electronic eigenbasis, and U is the electronic interaction kernel [30]. We once againsee that all elements have either two or four electronicstate labels, which motivates us to consider these labelspairwise, as transitions. If there are N el electronic states,then d µ is a vector in C N and both ˆ χ and ˆ U are N × N matrices (or rank-4 tensors) [31]. The inverse in Eq. (28)should be understood as a matrix inversion in this spaceand I is the identity matrix.The phonon self-energy ∆ˆ ω µν = (ˆ ω − ˆ ω ) µν is a matrix in phonon-branch space. The fluctuation di-agnostics looks at the contribution coming from a pair ofelectronic states k , k , which we denote by ∆ˆ ω µν,k k .Looking at a symmetrized combination of the contri-butions from k and k , we find √ m µ m ν ∆ˆ ω µν,k k (30)= − ˜ d µk k χ k k ,k k d ∗ νk k − d µk k χ k k ,k k ˜ d ∗ νk k (31)= − (cid:88) k k d µk k (cid:20) ˆ χ I − ˆ U ˆ χ (cid:21) k k ,k k d ∗ νk k − . . . (32)= − (cid:88) k k ,k k d µk k (cid:20) ˆ χ I − ˆ U ˆ χ ˆ P k k (cid:21) k k ,k k d ∗ νk k − . . . (33) ≡ (cid:104) d µ , ˆ P k k d ν (cid:105) χ + (cid:104) ˆ P k k d µ , d ν (cid:105) χ , (34)where we have introduced the projection onto k k :[ ˆ P k k ] k k ,k k = δ k k ,k k δ k k ,k k . In the last line,we have introduced the bilinear form corresponding tothe operator ˆ χ , with ˆ χ = − ˆ χ I − ˆ U ˆ χ . (35)In Appendix B, we proof that ˆ χ is a positive-definitesymmetric real matrix, so that the corresponding bilin-ear form is an inner product. ˆ χ can be understood asthe susceptibility and the fixed sign of ˆ χ is then relatedto thermodynamic stability, as discussed below Eq. (19).From this, our desired results follow directly, namely that∆ˆ ω µµ,k k ≥ ω µµ = (cid:80) k k ∆ˆ ω µµ,k k = (cid:104) d qµ , d qν (cid:105) χ / √ m µ m ν ≥ k and k label the electronic eigenstates.Here, these are the molecular orbitals. In lattice systems,momentum is a good quantum number and the electronicstates are labeled by momentum k and a band index. Themomentum analysis is simplified by setting k = k + q and by then observing that the entire equation is diagonalin q . IV. COMPUTATIONAL DETAILS
The ab initio calculations in this paper are realizedusing
Quantum ESPRESSO [32, 33]. We apply thegeneralized gradient approximation (GGA) by Perdew,Burke, and Ernzerhof (PBE) [34, 35] and optimizednorm-conserving Vanderbilt pseudopotentials [36] fromthe
PseudoDojo pseudopotential table [37]. Core elec-trons are incorporated into the pseudopotential. We usea large Gaussian occupation smearing of 0.1 Ry. Forcesare minimized to below 10 µ Ry/Bohr.
V. NITROGEN MOLECULE
We consider an N molecule, consisting of two iden-tical N ions aligned along the z -axis. Explicitly, theatomic coordinates are R µ = (0 , , − a/ , , , a/ a = 1 . × ω µν is a diagonal matrix. For thecalculation of the bare phonons, we have fixed all elec-tronic energy levels of Fig. 2. The resulting cDFPT andDFPT spectra are shown in Fig. 2. Due to rotation sym-metry in the x - y -plane, the eigenmodes with x and y displacements come in degenerate pairs.The five DFPT modes with vanishing energy are Gold-stone modes: two spontaneously broken rotation symme-tries (around x and y ) and three spontaneously broken electronicenergies2520151050510 ( e V ) ( e V ) x + x , y + y z + z x x , y y z z FIG. 2. Electronic and phononic spectra of an N molecule.Note that the 1s electrons are incorporated into the pseudopo-tential in our calculations, so they are not shown here. translation symmetries ( x , y , z ). In contrast, cDFPTconsiders a system where the electronic density is fixed.These constraints have already explicitly broken the sym-metries, see Sec. II D, so there is no spontaneous symme-try breaking and there are no bare Goldstone modes.This is visible in the cDFPT spectrum, all five barephonon modes have a finite energy.As predicted in Sec. II C, the frequency of the dressedphonons in Fig. 2 is reduced compared to the bare modes,although this effect is small and hardly visible for the z − z mode. A. Fluctuation diagnostics
In fluctuation diagnostics, one looks at the contribu-tion of individual electronic states to total change inphonon energy, thereby extracting the relevant screen-ing processes. Since there is no phonon mode mixing inN by symmetry, we can analyze the six phonon-diagonalelements ∆ˆ ω µµ one by one.The first ingredient is the bare susceptibility χ , shownin Fig. 3, quantifying the electronic transitions possible inthe system. Due to the Pauli principle, only transitionsacross the Fermi level are allowed and the susceptibilityis largest for pairs of states close to the Fermi level.In addition to the purely electronic χ , there is alsothe displacement potential. Figure. 4 shows d (note thatthis is shorthand for the product of a bare and dresseddeformation potential). We observe that d is real, pos-itive, and a symmetric matrix in electronic space. Be-cause of the rotation symmetry, the displacement poten-tial is identical for the x and y phonon modes (note thatwe have summed over degenerate electronic states). Themagnitude of d for the bond-stretching mode is approx-imately one order of magnitude larger than for the othermodes. The displacement potential for this mode is qual- x translation y translation z translation rotation around y rotation around x bond stretching x + x y + y z + z x − x y − y z − z (+1,0,0,+1,0,0) (0,+1,0,0,+1,0) (0,0,+1,0,0,+1) (+1,0,0, − − − . - . e V - . e V - . e V - . e V . e V . e V . e V -22.1 eV-7.3 eV-5.4 eV-4.0 eV4.2 eV5.9 eV7.2 eV FIG. 3. Bare susceptibility χ of N . The susceptibility isnon-zero only for transitions across the Fermi level. itatively different, since it is largely diagonal, whereas theother modes only have off-diagonal elements. Physically,a diagonal displacement potential means that ionic dis-placements change the energy of the Kohn-Sham orbitals,whereas off-diagonal elements describe displacement in-duced hybridization between orbitals. Most matrix ele-ments of the electron-phonon coupling are zero, reflectingselection rules.The total phonon renormalization √ m µ m ν ∆ˆ ω µµ is ob-tained by multiplying χ and d and summing over theintermediate electronic states. Fluctuation diagnosticsconsiders these summands one by one. The results so farshow that χ is restricted to transitions across the Fermilevel. The bond-stretching mode has no correspondingelectron-phonon coupling, so its energy is barely renor-malized. On the other hand, the five other modes aresubstantially renormalized. In fact, this renormalizationis responsible for breaking Goldstone’s theorem.From Fig. 4, it is clear that only very few combinationsof electronic states contribute to the phonon renormal- ization. For the three uniform translations, it is worth-while to analyze in more detail which electronic transi-tions are responsible for the phonon renormalization andthe absence of Goldstone’s theorem in the bare phonons.In all cases, it is the electronic state at 4.2 eV (LUMO)that is responsible for the screening, combined with ei-ther the − . − . x + x is odd under the reflection x (cid:55)→ − x , thetwo-fold degenerate 4.2 eV states are a combination ofp x and p y orbitals and thus have a p x component that isodd under this reflection, whereas − . z orbitals and is even under the reflection.This makes (cid:104) p x | ( x + x ) | sp z (cid:105) even and the transition isallowed. Note that (cid:104) p x | ( x + x ) | sp z (cid:105) is also even underthe reflections y (cid:55)→ − y and z (cid:55)→ − z . The same argumentholds for y + y and p y . Finally, for z + z , the stateat − . x and p y orbitals and the combination (cid:104) p x | ( z + z ) | p x (cid:105) is even under all three inversions and thus allowed. Thestructure of the upper panels in Fig. 4 reflects the dipoleselection rule for uniform translations: most matrix ele-ments are zero by symmetry.This also provides a recipe for choosing the electronictarget space in such a way that specific partially screeneduniform translation modes indeed have zero energy. If thestate at − . (cid:48) , then the partially screened phononmodes would already include this coupling. No furthercoupling to x + x is dipole-allowed, so the partiallyscreened x + x phonon needs to have zero energy. Sim-ilarly, integrating out the − . z + z phonon. This con-struction of Goldstone-preserving target spaces based ondipole selection rules can also be applied to crystallinematerials.Finally, we observe that both − χ and d are alwayspositive, so that every summand in the fluctuation diag-nostics gives a positive contribution. This is a numericalconfirmation of the earlier proof of positivity. VI. BENZENE
We now move on to a more complicated molecule, ben-zene (C H ). With twelve atoms, this molecule has 36phonon modes in total. All atoms lie in a single plane d , x + x d , y + y d , z + z . e V . e V . e V . e V . e V . e V . e V d , x x . e V . e V . e V . e V . e V . e V . e V d , y y . e V . e V . e V . e V . e V . e V . e V d , z z ( e V / Å ) FIG. 4. Displacement potential d in N . Every panel stands for the coupling to a single phonon mode, there is no mode-modecoupling in N . The colored matrix shows which electronic states couple to the phonon. Degenerate electronic states have beensummed over. ( z = 0). By symmetry, in-plane and out-of-plane dis-placements decouple. We restrict our attention to out-of-plane displacements that are furthermore symmetric un-der 180 degree rotation around the z -axis. This leaves asubspace consisting of six phonon modes. Figure 5 showsthe DFPT dynamical matrix ˆ ω in the basis of cDFPTeigenmodes. The presence of off-diagonal elements showsthat, unlike in N , there is mode-mode coupling. In par-ticular, we find coupling between modes that differ inthe relative direction of the C and H atoms [39]. In thecDFPT eigenbasis, the motion of the C and H ions isalmost completely decoupled. It is the electronic chem-ical bonding that is responsible for the coherent motionof both types of atoms and this bonding is frozen out incDFPT. An example of this is the z -translation DFPTGoldstone mode. This mode is a linear combination ofthe cDFPT modes at 0.0125 eV (translation of C atoms)and 0.0439 eV (translation of H atoms).Comparing the DFPT energies inside the matrix withthe cDFPT eigenenergies (to the left) shows that themagnitude of the screening is substantial. Note that neg-ative off-diagonal elements do appear in the matrix; pos-itivity is only guaranteed for the eigenvalues of ˆ ω ,ˆ ω , and ˆ ω − ˆ ω . A. Partially dressed
So far, we have considered the effect of electronicscreening on phonons by either allowing or forbiddingscreening from all electronic states. The resulting modesare called dressed and bare phonons, respectively. For es-tablishing low-energy models, one is frequently interestedin an intermediate object, the partially dressed phonons.In that case, screening by most electrons is allowed andonly transitions within a small “target” subspace close tothe Fermi level are excluded.For benzene, this low-energy subspace is spanned bythe six p z orbitals – analogous to the tight-bindingdescription of graphene – and there is a considerableCoulomb interaction between these electrons. Minimal-ist models of benzene consisting of six electronic statestherefore regularly feature as a testbed for investigatingelectron-electron interactions [40–46]. Here, we insteadfocus on the electron-ion interactions of this subspace.Figure 6 shows the phonon frequencies of both thefully screened DFPT phonons and the partially screenedcDFPT phonons. Their difference shows the effect ofscreening by the p z orbitals on the phonons. An im-portant observation is that out-of-plane phonons are not (eV ) FIG. 5. Selected phonons in benzene. The dressed (DFPT)dynamical matrix in the basis of bare (cDFPT) eigenmodes.The dynamical matrix is block diagonal and only a 6 × x , yz DFPT0.0000.0050.0100.0150.0200.0250.0300.0350.0400.045 ( e V ) cDFPT FIG. 6. DFPT and cDFPT phonons in benzene. x, y ( z )labels in-plane (out-of-plane) modes. The cDFPT phononsexclude electronic screening from within the six p z states.The out-of-plane modes are not affected by this constraintdue to the mirror symmetry selection rule. screened by the p z electrons, their frequencies are identi-cal in DFPT and cDFPT. This is caused by a symmetryselection rule, all p z orbitals are odd under mirror sym-metry and this implies d = 0.By comparing the partially and fully dressed dynam-ical matrices, we find that the most substantial renor-malization occurs for the phonon mode shown in Fig. 7,an in-plane mode where neighboring C atoms move inopposite directions. The figure also shows the fluctua-tion diagnostics of this mode. In the middle, the matrix d only has contributions coming from excitations be-tween the HOMO and LUMO levels, namely π ↔ π ∗ and π ↔ π ∗ . These are the only combinations of molec-ular orbitals that are allowed to be coupled to this modeby mirror symmetry. Since these excitations are also en-ergetically favorable for the susceptibility χ (right), theylead to a large overall renormalization of this phononmode. VII. CONCLUSION
We have constructed a general framework for ab ini-tio downfolding electron-ion problems within the Born-Oppenheimer approximation and have shown how thecDFPT fits into this framework. We have shown ana-lytically that electronic screening lowers the phonon en-ergies. Even the fluctuation diagnostics – the contribu-tions from specific electronic states – are sign definite.The constrained theory can explicitly break symmetries,thereby reducing the number of Goldstone modes. Dipoleselection rules can be used to construct a low-energy elec-tronic subspace where even the partially dressed Gold-stone modes are guaranteed to have zero energy. Wenote that the cDFPT phonon dispersion in Fig. 2(a) ofRef. 16 indeed satisfies Goldstone’s theorem.We have illustrated the theorem with cDFPT formolecules (nitrogen and benzene), since these are amongthe simplest and clearest examples of electron-ion cou-pling. In particular, they provide a clear view on theorbital structure of the theory.
ACKNOWLEDGMENTS
The authors would like to thank J¨org Kr¨oger, SusanK¨oppen, and Filippo Balzaretti for useful discussions. Fi-nancial support by the Deutsche Forschungsgemeinschaft(DFG) through GRK 2247 and computational resourcesof the North-German Supercomputing Alliance (HLRN)are gratefully acknowledged. EvL is supported by theCentral Research Development Fund of the University ofBremen.0 d FIG. 7. Left: An in-plane phonon (DFPT eigenmode) with substantial coupling to the p z electrons. Middle: d for this mode.Right: Electronic susceptibility χ . d and χ are both visualized in terms of the eigenstates of the electronic p z space. Thecontribution to the phonon self-energy is obtained as the pointwise product of these two squares. Appendix A: Manifolds and derivatives
In DFPT, the electronic configuration is described bythe density operator ˆ ρ . The Kohn-Sham electronic en-ergy functional is E (ˆ ρ ) = Trˆ ρH , where H is the Kohn-Sham Hamiltonian. Let us initially assume that H isindependent of ˆ ρ . One could anticipate that the Hes-sian matrix, the second derivative of E with respect to ˆ ρ ,would be zero. However, at this point it is important torealize that the derivatives are taken on the manifold ofdensity matrices and the Hessian on a curved manifoldcontains additional terms.To see how this works, we consider a minimal sys-tem of two electronic levels with single-particle energies (cid:15) a < (cid:15) b and corresponding states | a (cid:105) , | b (cid:105) . The Hamil-tonian is H = (cid:15) a | a (cid:105) (cid:104) a | + (cid:15) b | b (cid:105) (cid:104) b | . The ground state isgiven by the ground-state density operator ˆ ρ . If (cid:15) a and (cid:15) b are on the same side of the Fermi level, the systemis completely filled or empty. The most interesting sit-uation occurs when (cid:15) a < E f < (cid:15) b . In that case, theground-state ˆ ρ = | a (cid:105) (cid:104) a | is the projection operator onto a and the total number of electrons is 1. The mani-fold of allowed density matrices is given by all densitymatrices with total density 1. It is easy to see thatˆ ρ θ = (cos θ | a (cid:105) + sin θ | b (cid:105) )(cos θ (cid:104) a | + sin θ (cid:104) b | ) lies in thismanifold for any θ . The corresponding energy is E (ˆ ρ θ ) = Trˆ ρ θ H (A1)= (cid:15) a cos θ + (cid:15) b sin θ, (A2)which has first derivative zero at θ = 0, and the secondderivative at θ = 0 is d E dθ (cid:12)(cid:12)(cid:12) θ =0 = 2( (cid:15) b − (cid:15) a ) = 2( χ ab,ab ) − (A3)= d E d ˆ ρ ab + d E d ˆ ρ ba . (A4)The inverse bare susceptibility at T = 0 in Eq. (A3)enters the phonon self-energy. Generalizing the result to higher-dimensional spaces shows that off-diagonal terms ∂ ab ∂ cd vanish, so χ ab,cd is a diagonal matrix, as stated inthe main text.So far, we have assumed that the Kohn-Sham Hamil-tonian H is independent of ˆ ρ . In reality, DFT is a self-consistent theory and the Kohn-Sham Hamiltonian con-tains the Hartree and exchange-correlation potentials. Itis generally possible to write d E d ˆ ρ αβ d ˆ ρ γδ (cid:12)(cid:12)(cid:12) ˆ ρ = ( χ ) − αβ,γδ − ˆ U αβ,γδ , (A5)which here merely acts as the definition of ˆ U . ˆ U accountsfor the electron-electron interactions that are not presentin the auxiliary Kohn-Sham system. In cDFPT calcula-tions, ˆ U is formally incorporated into ˜ d (see Fig. 1) andnever calculated explicitly. Still, this formal relation isnecessary for some of the proofs in the main text. Appendix B: Sign of phonon self-energy
From the definition, ˆ χ is diagonal (as a N × N matrix) and thus also symmetric. The diagonal ele-ments are given by the Lindhard expression and arenegative or zero, since f is a decreasing function.Thus, − ˆ χ is positive semi-definite. The matrix ˆ U isreal and symmetric. From this, it follows that ˆ χ = − ˆ χ / ( I − ˆ U ˆ χ ) is also real and symmetric, which canbe checked order by order as ( ˆ χ ) ˆ U ˆ χ . . . ˆ χ ˆ U ˆ χ ) T =( ˆ χ ) T ˆ U T ( ˆ χ ) T . . . ( ˆ χ ) T ˆ U T ( ˆ χ ) T = ˆ χ ˆ U ˆ χ . . . ˆ χ ˆ U ˆ χ .As long as the geometric series is convergent, which cor-responds to thermodynamic stability, it does not changethe sign of eigenvalues and ˆ χ will be positive semi-definiteand symmetric.A further detail that is necessary for our proof: P is anorthogonal projection with respect to the original innerproduct but not with respect to the inner product definedby ˆ χ , so it is not immediately trivial that (cid:104) d, P d (cid:105) χ ≥ χ is written without hat for convenience. Since1 P is a projection operator, we can write d = α + β with α ≡ P d , P α = P d = P d = α and β ≡ d − α , P β = P d − P α = α − α = 0.We wish to show that d T χP d = α T χα + β T χα ? ≥ . (B1)If β T χα >
0, this is proven immediately, since α T χα > x = α + λβ with a real number λ so that0 ≤ x T χx = α T χα + 2 λβ T χα + λ β T χβ = α T χα + β T χα. (B2)The last equality is used to solve for λ and a real solu-tion λ can be found as long as the discriminant D/ β T χα ) − ( β T χα )( β T χβ ) is positive. But we were study-ing the case β T χα <
0, so indeed
D > λ can be chosenappropriately, and the proof of the lemma is finished. [1] M. Imada and T. Miyake, Electronic structure calcu-lation by first principles for strongly correlated elec-tron systems, J. Phys. Soc. Jpn. , 112001 (2010),arXiv:1009.3851.[2] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko,O. Parcollet, and C. A. Marianetti, Electronic struc-ture calculations with dynamical mean-field theory, Rev.Mod. Phys. , 865 (2006), arXiv:cond-mat/0511085.[3] G. Giovannetti, M. Casula, P. Werner, F. Mauri, andM. Capone, Downfolding electron-phonon Hamiltoniansfrom ab initio calculations: Application to K picene,Phys. Rev. B , 115435 (2014), arXiv:1406.4108.[4] Y. Nomura and R. Arita, Ab initio downfolding forelectron-phonon-coupled systems: Constrained density-functional perturbation theory, Phys. Rev. B , 245108(2015), arXiv:1509.01138.[5] M. R¨osner, C. Steinke, M. Lorke, C. Gies, F. Jahnke, andT. O. Wehling, Two-dimensional heterojunctions fromnonlocal manipulations of the interactions, Nano Lett. , 2322 (2016).[6] J. Hall, N. Ehlen, J. Berges, E. van Loon, C. van Efferen,C. Murray, M. R¨osner, J. Li, B. V. Senkovskiy, M. Hell,M. Rolf, T. Heider, M. C. Asensio, J. Avila, L. Plucinski,T. Wehling, A. Gr¨uneis, and T. Michely, Environmentalcontrol of charge density wave order in monolayer 2H-TaS , ACS Nano , 10210 (2019).[7] V. I. Anisimov, J. Zaanen, and O. K. Andersen, Bandtheory and Mott insulators: Hubbard U instead of stoner I , Phys. Rev. B , 943 (1991).[8] F. Aryasetiawan, M. Imada, A. Georges, G. Kotliar,S. Biermann, and A. I. Lichtenstein, Frequency-dependent local interactions and low-energy effectivemodels from electronic structure calculations, Phys. Rev.B , 195104 (2004), arXiv:cond-mat/0401620.[9] O. Gunnarsson, Superconductivity in fullerides, Rev.Mod. Phys. , 575 (1997), arXiv:cond-mat/9611150.[10] M. Capone, M. Fabrizio, C. Castellani, and E. Tosatti,Strongly correlated superconductivity, Science , 2364(2002), arXiv:cond-mat/0207058.[11] Y. Nomura, S. Sakai, M. Capone, and R. Arita, Unifiedunderstanding of superconductivity and Mott transitionin alkali-doped fullerides from first principles, Sci. Adv. , e1500568 (2015), arXiv:1505.05849.[12] Y. Nomura, S. Sakai, M. Capone, and R. Arita,Exotics-wave superconductivity in alkali-doped ful-lerides, J. Phys. Condens. Matter , 153001 (2016),arXiv:1512.05755.[13] Y. Nomura, Ab Initio Studies on Superconductivity in Alkali-Doped Fullerides (Springer Singapore, 2016).[14] R. Arita, T. Koretsune, S. Sakai, R. Akashi, Y. Nomura,and W. Sano, Nonempirical calculation of superconduct-ing transition temperatures in light-element supercon-ductors, Adv. Mater. , 1602421 (2017).[15] D. Novko, Broken adiabaticity induced by Lifshitz tran-sition in MoS and WS single layers, Commun. Phys. ,30 (2020), arXiv:1907.04766.[16] J. Berges, E. G. C. P. van Loon, A. Schobert, M. R¨osner,and T. O. Wehling, Ab initio phonon self-energies andfluctuation diagnostics of phonon anomalies: Lattice in-stabilities from Dirac pseudospin physics in transitionmetal dichalcogenides, Phys. Rev. B , 155107 (2020),arXiv:1911.02450.[17] X. Gonze, D. C. Allan, and M. P. Teter, Dielectric tensor,effective charges, and phonons in α -quartz by variationaldensity-functional perturbation theory, Phys. Rev. Lett. , 3603 (1992).[18] A. Putrino, D. Sebastiani, and M. Parrinello, General-ized variational density functional perturbation theory,J. Chem. Phys. , 7102 (2000).[19] K. Refson, P. R. Tulip, and S. J. Clark, Variationaldensity-functional perturbation theory for dielectrics andlattice dynamics, Phys. Rev. B , 155114 (2006).[20] We assume that the minimum is unique.[21] X. Gonze and J.-P. Vigneron, Density-functional ap-proach to nonlinear-response coefficients of solids, Phys.Rev. B , 13120 (1989).[22] X. Gonze, Perturbation expansion of variational princi-ples at arbitrary order, Phys. Rev. A , 1086 (1995).[23] X. Gonze, Adiabatic density-functional perturbation the-ory, Phys. Rev. A , 1096 (1995).[24] Ω is a manifold and thus locally equivalent to a vectorspace, so the derivative δψ can be decomposed into com-ponents.[25] O. Gunnarsson, T. Sch¨afer, J. P. F. LeBlanc, E. Gull,J. Merino, G. Sangiovanni, G. Rohringer, and A. Toschi,Fluctuation diagnostics of the electron self-energy: Ori-gin of the pseudogap physics, Phys. Rev. Lett. ,236402 (2015), arXiv:1411.6947.[26] S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Gi-annozzi, Phonons and related crystal properties fromdensity-functional perturbation theory, Rev. Mod. Phys. , 515 (2001), arXiv:cond-mat/0012092.[27] The name dipole selection rule originates in the relationbetween the matrix elements of the momentum operatorand the position operator. If ˆ H = ˆ p / m + V (ˆ r ), then[ ˆ H, ˆ r ] = − i ˆ p/m and ( E a − E b ) (cid:104) a | ˆ r | b (cid:105) = (cid:104) a | [ ˆ H, ˆ r ] | b (cid:105) = − i/m (cid:104) a | ˆ p | b (cid:105) .[28] F. Giustino, Electron-phonon interactions from firstprinciples, Rev. Mod. Phys. , 015003 (2017),arXiv:1603.06965.[29] L. Hedin and S. Lundqvist, Effects of electron-electronand electron-phonon interactions on the one-electronstates of solids, in Solid State Physics: Advances in Re-search and Applications , edited by F. Seitz, D. Turnbull,and H. Ehrenreich (Academic Press, New York, London,1969).[30] This electronic interaction is present since we need toconsider the self-consistent response of the electronic sys-tem, See Ref. 4 for more details. In the current proof, weonly require that this denominator does not change thesign of the eigenvalues. In other words, the electronic sys-tem should be thermodynamically stable [47].[31] M. Kaltak,
Merging GW with DMFT , Ph.D. thesis, Uni-versit¨at Wien (2015).[32] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car,C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococ-cioni, I. Dabo, A. D. Corso, S. d. Gironcoli, S. Fabris,G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis,A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari,F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello,L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero,A. P. Seitsonen, A. Smogunov, P. Umari, and R. M.Wentzcovitch,
Quantum ESPRESSO : A modular andopen-source software project for quantum simulations ofmaterials, J. Phys. Condens. Matter , 395502 (2009),arXiv:0906.2569.[33] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau,M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni,D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo,A. D. Corso, S. d. Gironcoli, P. Delugas, R. A. DiStasio,A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer,U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawa-mura, H.-Y. Ko, A. Kokalj, E. K¨u¸c¨ukbenli, M. Lazzeri,M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. Otero-de-la-Roza, L. Paulatto, S. Ponc´e,D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P.Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser,P. Umari, N. Vast, X. Wu, and S. Baroni, Advancedcapabilities for materials modelling with QuantumESPRESSO , J. Phys. Condens. Matter , 465901(2017), arXiv:1709.10010.[34] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalizedgradient approximation made simple, Phys. Rev. Lett. , 3865 (1996).[35] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalizedgradient approximation made simple [Phys. Rev. Lett. 77, 3865 (1996)], Phys. Rev. Lett. , 1396 (1997).[36] D. R. Hamann, Optimized norm-conserving Vander-bilt pseudopotentials, Phys. Rev. B , 085117 (2013),arXiv:1306.4707.[37] M. J. van Setten, M. Giantomassi, E. Bousquet, M. J.Verstraete, D. R. Hamann, X. Gonze, and G. M.Rignanese, The PseudoDojo: Training and gradinga 85 element optimized norm-conserving pseudopoten-tial table, Comput. Phys. Commun. , 39 (2018),arXiv:1710.10138.[38] Here, the phonon energies stand only for the secondderivative of the classical potential-energy surface of themolecule in the Born-Oppenheimer approximation.[39] In total, there are twelve atoms and thus also 12 out-of-plane modes. The C rotational symmetry divides theseup into six pairs of modes, with mode-mode couplingallowed only within these pairs. Three of the six pairsare shown in Fig. 5.[40] R. Pariser and R. G. Parr, A semi-empirical theory ofthe electronic spectra and electronic structure of complexunsaturated molecules. II, J. Chem. Phys. , 767 (1953).[41] J. A. Pople, The electronic spectra of aromatic moleculesII: A theoretical treatment of excited states of alternanthydrocarbon molecules based on self-consistent molecularorbitals, Proc. Phys. Soc. A , 81 (1955).[42] A. Valli, G. Sangiovanni, O. Gunnarsson, A. Toschi,and K. Held, Dynamical vertex approximation fornanoscopic systems, Phys. Rev. Lett. , 246402 (2010),arXiv:1003.2630.[43] M. Sch¨uler, M. R¨osner, T. O. Wehling, A. I. Lichtenstein,and M. I. Katsnelson, Optimal Hubbard models for mate-rials with nonlocal Coulomb interactions: Graphene, sil-icene, and benzene, Phys. Rev. Lett. , 036601 (2013),arXiv:1302.1437.[44] H. J. Changlani, H. Zheng, and L. K. Wagner, Density-matrix based determination of low-energy model Hamil-tonians from ab initio wavefunctions, J. Chem. Phys. , 102814 (2015), arXiv:1504.03704.[45] P. Pudleiner, P. Thunstr¨om, A. Valli, A. Kauch, G. Li,and K. Held, Parquet approximation for molecules:Spectrum and optical conductivity of the Pariser-Parr-Pople model, Phys. Rev. B , 125111 (2019),arXiv:1812.04962.[46] Y. in ’t Veld, M. Sch¨uler, T. O. Wehling, M. I. Katsnel-son, and E. G. C. P. v. Loon, Bandwidth renormalizationdue to the intersite Coulomb interaction, J. Phys. Con-dens. Matter , 465603 (2019), arXiv:1901.11257.[47] E. C. Stoner, Collective electron ferromagnetism II. En-ergy and specific heat, Proc. R. Soc. Lond. A169