Thermodynamic electric quadrupole moments of nematic phases from first-principles calculation
TThermodynamic electric quadrupole moments of nematic phases from first-principles calculation
Taisei Kitamura, ∗ Jun Ishizuka, † Akito Daido, ‡ and Youichi Yanase
1, 2, § Department of Physics, Graduate School of Science, Kyoto University, Kyoto 606-8502, Japan Institute for Molecular Science, Okazaki 444-8585, Japan (Dated: February 22, 2021)The electronic nematic phase emerging with spontaneous rotation symmetry breaking is a central issue ofmodern condensed matter physics. In particular, various nematic phases in iron-based superconductors andhigh- T c cuprate superconductors are extensively studied recently. Electric quadrupole moments (EQMs) areone of the order parameters characterizing these nematic phases in a unified way, and elucidating EQMs is akey to understanding these nematic phases. However, the quantum-mechanical formulation of the EQMs incrystals is a nontrivial issue because the position operators are non-periodic and unbound. Recently, the EQMshave been formulated by local thermodynamics, and such thermodynamic EQMs may be used to characterizethe fourfold rotation symmetry breaking in materials. In this paper, we calculate the thermodynamic EQMsin iron-based superconductors LaFeAsO and FeSe as well as a cuprate superconductor La CuO by a first-principles calculation. We show that owing to the orbital degeneracy the EQMs in iron-based superconductorsare mainly determined by the geometric properties of wave functions. This result is in sharp contrast to thecuprate superconductor, in which the EQMs are dominated by distortion of the Fermi surface. I. INTRODUCTION
In recent years, the nematic phases which spontaneouslybreak fourfold rotation ( C ) symmetry are attracting a lotof interest in condensed matter physics. For example, iron-based superconductors such as LaFeAsO , FeSe , andBaFe As undergo nematic order with an electronic origin,although it is accompanied by the tetragonal-orthorhombicstructural phase transition. The relations of nematic order,superconductivity, and magnetism have been of central in-terest in the research of iron-based superconductors in thepast decade . The nematic order in cuprate superconduc-tors is also a topic of interest, motivated by recent experimen-tal indications . Although vast studies have been devoted,comprehensive clarification of the nematic order and its rela-tion to the pseudogap phase and superconductivity remains anongoing issue. For an origin of the nematic order in cuprates,the charge density wave (CDW) order and bond order aswell as other exotic order such as loop current and pair densitywave have been proposed. Because of a possible inter-play with the high- T c superconductivity in the two categories,namely, iron-based and cuprate superconductors, the nematicphase in strongly correlated electron systems is a central issueof modern condensed matter physics. In particular, the sim-ilarities and differences between iron-based superconductorsand cuprate superconductors are issues to be solved.Although various nematic order parameters have been pro-posed, a ubiquitous feature is the C -symmetry breaking.The electric quadrupole moments (EQMs) are one of thefundamental quantities which naturally characterize the C -symmetry breaking . Therefore studying EQMs is a key tounderstanding the nematic phases, and may also give an in-sight into the relations to superconductivity.EQMs are originally introduced in classical electromag-netism . However, contrary to their apparently simple form,a naive extension to periodic crystals is problematic due to thedifficulties in the treatment of the position operator. This isalso the case for magnetic multipole moments. In contrast tothe modern theory of electric polarization where the approach based on Wannier functions is successful , EQMs of theWannier functions are gauge dependent unless appropriatesymmetries are preserved . On the other hand, EQMs havebeen well formulated as topological invariants of higher ordertopological crystalline insulators . However, obtained re-sults are valid only in the presence of crystalline symmetriessuch as the C symmetry. Therefore, these approaches arenot useful for evaluating EQMs that emerge with spontaneous C -symmetry breaking.Although the above quantum mechanical approach is stillan ongoing issue, recent thermodynamic approaches to theelectric/magnetic multipole moments successfully ob-tained gauge-invariant and unit-cell independent formulas.In particular, the EQMs obtained by the thermodynamic ap-proach, which we call thermodynamic EQMs, are well definedeven without crystalline symmetries such as the C symme-try, in contrast to the EQMs formulated as topological invari-ants. Therefore, using the thermodynamic EQMs, now wecan study the nematic phases of iron-based superconductorsand cuprate superconductors in a unified manner.In this paper, after showing the failures of the formulationsby electromagnetism and Wannier function methods, we cal-culate the thermodynamic EQMs of two iron-based super-conductors, LaFeAsO, FeSe, and a cuprate superconductorLa CuO using first-principles calculations. For the nematicorder parameters, the orbital order and momentum-dependentorbital polarization are assumed for LaFeAsO and FeSe, re-spectively, in accordance with theoretical proposals .For La CuO , we examine the d x − y -wave bond order aswell as the orbital order of O2 p x and O2 p y orbitals for acomparison. The results reveal a difference in the EQMs be-tween the iron-based superconductors and cuprate supercon-ductors. The Fermi-sea term in the thermodynamic EQM isdominant in iron-based superconductors, owing to the uniqueband structures and associated geometric properties. In con-trast, the thermodynamic EQM in the cuprate superconductoris dominated by the Fermi-surface term arising from the dis-tortion of the band structure. a r X i v : . [ c ond - m a t . s t r- e l ] F e b II. ELECTROMAGNETIC EQMS IN CRYSTALS
We begin with the definition of EQMs in classical electro-magnetism. With the multipole expansion, the scalar potentialis described as , φ ( r ) = 14 πε (cid:90) d r (cid:48) ∞ (cid:88) l =0 ρ ( r (cid:48) ) r (cid:48) l r l +1 P l (ˆ r · ˆ r (cid:48) ) . (1)Here, ε is the dielectric constant, P l ( x ) are the Legendre poly-nomials, ˆ r is the unit vector of r , i.e. ˆ r = r / | r | , and ρ ( r ) isthe charge density. For simplicity, we adopt the units withelectric charge e = 1 and lattice volume V cell = 1 . We focuson the component of l = 2 , which is represented as φ (2) ( r ) = 14 π(cid:15) (cid:88) ij r (4 r i r j − r δ ij ) Q ij , (2)with the electromagnetic EQMs Q ij = (cid:90) d r r i r j ρ ( r ) . (3)The electromagnetic EQMs are determined by the quadrupoledistribution of the charge density. The other multipole mo-ments are defined as well in a similar manner.Here, we try to evaluate the electromagnetic EQMs of elec-trons on the periodic crystal lattice. Difficulties due to thenon-periodic and unbound position operators will becomemanifest. We consider the spinless case for simplicity. Anidea to avoid the unboundedness of the position operator,which can be infinite in the thermodynamic limit, is to con-sider the EQMs in a unit cell. The EQMs are redefined as Q cell ij = (cid:90) cell dr d r i r j (cid:104) ψ † ( r ) ψ ( r ) (cid:105) = 1 V (cid:88) nm (cid:88) kk (cid:48) (cid:90) cell d r r i r j (cid:104) c † n ( k ) c m ( k (cid:48) ) (cid:105)× e − i ( k − k (cid:48) ) · r u ∗ n ( k , r ) u m ( k (cid:48) , r ) , (4)with ψ † ( r ) being the creation operator of electrons at r . In the second line, we expanded ψ ( r ) by the peri-odic part of the Bloch wave function u n ( k , r ) : ψ ( r ) =1 √ V (cid:80) n (cid:80) k c n ( k ) e i k · r u n ( k , r ) . Here, m, n are band in-dices, k is the wave vector, V is the volume of the system, and d is the dimension of the system. The momentum sum is car-ried out in the first Brillouin zone. In this notation, u n ( k , r ) are normalized as (cid:82) dr d u ∗ m ( k , r ) u n ( k , r ) = δ nm . For freeelectron systems with band dispersion (cid:15) n ( k ) , the relation (cid:104) c † n ( k ) c m ( k (cid:48) ) (cid:105) = f ( (cid:15) n ( k )) δ nm δ k , k (cid:48) , (5)leads to Q cell ij = (cid:88) n (cid:90) BZ dk d (2 π ) d (cid:90) cell dr d r i r j f ( (cid:15) n ( k )) u ∗ n ( k , r ) u n ( k , r ) . (6)Here, f ( (cid:15) ) = ( e (cid:15)/T + 1) − represents the Fermi distributionfunction for the temperature T . Fe A Fe B Fe A Fe B x x y y o o(a) (b) FIG. 1. Unit cell of LaFeAsO and choices of real space coordinates.For the case (a), the electromagnetic EQM, Q x − y = Q x − Q y ,vanishes, while Q x − y is finite for the case (b). Thus, the EQMs ofthe unit cell are not well defined. The EQMs of the unit cell, Eq. (6), depends on the originof real space coordinates. For concreteness, we consider thetwo-sublattice systems as for LaFeAsO we study in Sec. IV.In a tight-binding model, the orbitals of electrons are localizedon atoms and we have a simple relation u ∗ n ( k , r ) u n ( k , r ) = (cid:104) u n ( k ) | r (cid:105) (cid:104) r | u n ( k ) (cid:105) = (cid:88) R (cid:88) sub = A,B | (cid:104) r sub | u n ( k ) (cid:105) | δ ( r − R − r sub ) , (7)where r A , r B are the sublattice positions within a unit celland R are the lattice points. Only the home unit cell R = 0 contributes to the integral (6).Examples of the coordinates are shown in Fig. 1. A EQMcharacterizing the nematic order is the inbalance between the x and y directions: Q cell x − y = Q cell x − Q cell y . This vanishesfor the coordinates in Fig. 1(a), since Fe and Fe are on thelines of x − y = 0 , while Q cell x − y is finite for the case ofFig. 1(b). Indeed, we obtain Q cell x − y (cid:39) − . at T = 0 . for LaFeAsO in the absence of the nematic order, when wechoose coordinates in Fig. 1(b). Because of these undesir-able properties, (1) Q cell ij depend on the coordinates, and (2) Q cell ij can be finite even in the absence of the nematic order,the electromagnetic EQMs of the unit cell are not suitable forquantifying the nematic order.We also discuss the EQMs of Wannier functions defined bythe moment of position operators: Q ( n ) ij = (cid:104) W n | ˆ r i ˆ r j | W n (cid:105) = − (cid:88) m ( (cid:54) = n ) (cid:90) BZ dk d (2 π ) d (cid:104) u n ( k ) | ∂ k i u m ( k ) (cid:105) (cid:104) u m ( k ) | ∂ k j u n ( k ) (cid:105)− (cid:90) BZ dk d (2 π ) d (cid:104) u n ( k ) | ∂ k i u n ( k ) (cid:105) (cid:104) u n ( k ) | ∂ k j u n ( k ) (cid:105) , (8)with the Wannier function of the n -th band, | W R n (cid:105) = (cid:90) dk d (2 π ) d e − i k · ( R − ˆ r ) | u n ( k ) (cid:105) . (9)Clearly, the second term of Eq. (8) is gauge dependent. There-fore, the EQMs of Wannier functions are unsuitable for eval-uating the EQMs of nematic phases.Instead of the above quantum-mechanical approaches, weadopt the thermodynamic approach as we explain in the nextsection. Note that the first term of Eq. (8) appears in the ther-modynamic EQMs . It is a part of the geometric term andgiven by the quantum metric , namely, the real part of the quantum geometric tensor. III. THERMODYNAMIC EQMS
In this section, we introduce the thermodynamic EQMs.The EQMs are recently formulated by the variation of the freeenergy density : dF ( r ) = ρ ( r ) dφ ( r ) + p i d [ ∂ i φ ( r )] + Q ij d [ ∂ i ∂ j φ ( r )] + O ( d [ ∇ φ ( r )] , [ dφ ( r )] ) . (10)While the charge density ρ ( r ) defined by the differential withrespect to the scalar potential φ ( r ) is regarded as thermody-namic electric monopole moment, the electric dipole p i andquadrupole moments Q ij are given by the differential withrespect to the spatially nonuniform scalar potential. The ther- modynamic EQMs are defined as the change of free energydensity by the nonuniform electric field, and therefore, it isnaturally related to the quadrupole charge distribution.The expressions for the thermodynamic EQMs Q ij are ob-tained as Q ij = (cid:88) n (cid:90) BZ d d k (2 π ) d (cid:34) g ijn ( k ) f ( (cid:15) n ( k )) − X ijn ( k ) (cid:90) ∞ (cid:15) n ( k ) d(cid:15)f ( (cid:15) ) − m − n ( k ) ij f (cid:48) ( (cid:15) n ( k )) (cid:35) . (11)Here, g ijn ( k ) is the quantum metric : g ijn ( k ) = (cid:88) m (cid:54) = n A inm ( k ) A jmn ( k ) + c . c ., (12)which is a counter-part of the Berry curvature. The quantummetric is the real part of the quantum geometric tensor ,while the imaginary part is the Berry curvature. Equation (11)reveals that the momentum integral of the quantum metricgives a part of the thermodynamic EQMs, while the integralof the Berry curvature is an intrinsic part of the anomalousHall conductivity . In the second term, X ijn ( k ) is given by X ijn ( k ) = − (cid:88) m (cid:54) = n A inm ( k ) A jmn ( k ) + c . c .(cid:15) n ( k ) − (cid:15) m ( k ) , (13)while m − n ( k ) ij = ∂ k i ∂ k j (cid:15) n ( k ) in the third term is the in-verse effective mass tensor. We adopted the following nota-tions ( ˆ H − µ ) | ψ n ( k ) (cid:105) = (cid:15) n ( k ) | ψ n ( k ) (cid:105) , (14) | u n ( k ) (cid:105) = e − i kr | ψ n ( k ) (cid:105) , (15) A inm ( k ) = − i (cid:104) u n ( k ) | ∂ k i u m ( k ) (cid:105) . (16) ˆ H and µ are the noninteracting Hamiltonian and the chemicalpotential, respectively. Note that Eq. (11) is valid only whenall bands are isolated. General expressions valid in the pres-ence of band touchings are provided in Ref. 67. Although weadopt single-particle Hamiltonian, the expressions can also beused for the many-body states after the mean-field approxima-tion. The EQMs formulated based on the thermodynamics aregauge-invariant and independent of unit-cell choices. Thus,difficulties of the electromagnetic EQMs in periodic crystalshave been solved. Therefore, we evaluate the thermodynamicEQMs in the representative nematic phases of cuprates andiron-based superconductors and discuss their origin.In Eq. (11), the first two terms are the contribution fromthe Fermi sea. These terms reflect the geometric properties ofthe Bloch electrons (12) and (13), and thus, they are calledthe Fermi-sea terms in the following. On the other hand,the third term is a property of Fermi surfaces and vanishesin insulators at zero temperature. This term is given by theanisotropy of the inverse effective mass tensor on Fermi sur-faces and called the Fermi-surface term. In the following sec-tions, we discuss the difference between iron-based supercon-ductors and cuprate superconductors from the perspective ofthe thermodynamic EQMs, based on the decomposition intothe Fermi-sea and Fermi-surface contributions. For more de-tails of Eq. (11), please see Ref. 67. IV. THERMODYANMIC EQMS FROMFIRST-PRINCIPLES CALCULATION
In this section, we calculate and discuss the thermodynamicEQMs using the first-principles calculation. Here, we focuson three representative high- T c superconductors, LaFeAsO,FeSe, and La CuO . The first-principles electronic structurecalculations are performed with using the WIEN
2k code , andthe tight-binding models based on the maximally localizedWannier functions are constructed by the WANNIER .The tight-binding Hamiltonian is represented as H LDA = (cid:88) k (cid:88) l,m,i,j,σ t lmij ( k ) c † liσ ( k ) c mjσ ( k ) , (17)where c † liσ ( k ) ( c liσ ( k )) is the creation (annihilation) operatorof the electrons with wave vector k , orbital l , sublattice i , andspin σ . The matrix elements t lmij ( k ) are given by the Fouriertransform of the hopping integrals, which are obtained fromthe WIEN
2k code. Here, we neglect the spin-orbit couplingfor simplicity, and evaluate the EQMs per spin. To calculatethe EQMs in the nematic phases we take into account phe-nomenological molecular fields for the nematic order param-eter. The total Hamiltonian is H = H LDA + ∆( T )Γ , (18)where ∆( T ) is the order parameter at the temperature T , and Γ is the molecular field of the bond order, orbital order, andso on. We specify Γ for each compound later. To calculate thetemperature dependence of the EQMs, we assume ∆( T ) = (cid:26) T > T s )∆ (cid:112) − T /T s ( T ≤ T s ) , (19)where T s is the phase transition temperature of the nematicorder accompanied by C -symmetry breaking. We set ∆ = T s = 0 . , roughly in accordance with the nematic transi-tion temperatures. Because all materials studied in this pa-per have quasi-two-dimensional electronic structures, we ig-nore the hopping integral along the z -direction. Thus, two-dimensional multi-orbital and multi-sublattice tight-bindingmodels are analyzed. Two EQMs, namely, Q x − y and Q xy can be finite in two-dimensional systems. In the following, weconsider the nematic phases accompanied by a finite Q x − y .Note that Q xy = 0 because of the mirror symmetry. A. LaFeAsO
Here we evaluate the EQMs in LaFeAsO. We construct 10-orbital tight-binding models for Fe3 d electrons. The presenceof two iron atoms in a unit cell doubles the number of orbitals,as × . It has been suggested that the origin of thenematic order in LaFeAsO is the orbital order of d xz and d yz electrons of irons . Thus, the molecular field isassumed as Γ = (cid:88) k (cid:88) σ,i =1 , (cid:104) c † d xz iσ ( k ) c d xz iσ ( k ) − c † d yz iσ ( k ) c d yz iσ ( k ) (cid:105) , (20)where d xz and d yz denote atomic orbitals.We examine two crystal structures. One is the tetragonalcrystal with the space group P /nmm and the other is the or-thorhombic crystal with the space group Cmma . The formeris realized in the high-temperature phase of LaFeAsO, while - (cid:1) (cid:1) - (cid:1) (cid:1) k x k y - (cid:1) (cid:1) - (cid:1) (cid:1) k x k y - (cid:1) (cid:1) - (cid:1) (cid:1) k x k y - (cid:1) (cid:1) - (cid:1) (cid:1) k x k y - (cid:1) (cid:1) - (cid:1) (cid:1) k x k y - (cid:1) (cid:1) - (cid:1) (cid:1) k x k y - (cid:1) (cid:1) - (cid:1) (cid:1) k x k y - (cid:1) (cid:1) - (cid:1) (cid:1) k x k y (a) (b)(c) (d) k Y k X FIG. 2. Fermi surfaces of tetragonal LaFeAsO [(a)-(b)] and or-thorhombic LaFeAsO [(c)-(d)]. (a) and (c) the normal state at T = 0 . , where ∆( T ) = 0 . (b) and (d) the orbital ordered stateat T = 0 . . the latter emerges below the structural transition temperature T s ≈ K. Lattice parameters are derived from Ref. 2. Al-though the structural transition is associated with the nematicorder, we independently study the two phenomena to clarifythe origin of the EQMs. The Fermi surfaces in the orbital-ordered state (normal state) on the tetragonal crystal latticeare shown in Fig. 2(b) (Fig. 2(a)), indicating sizable distortionof the Fermi surfaces due to the orbital order. The same plotsfor the orthorhombic crystal lattice are shown in Figs. 2(c)and 2(d). We see that the effects of the orthorhombic crystaldistortion on the Fermi surfaces are not prominent.
Fe1Fe2 𝑑 !" orbital 𝑑 orbital Y X
FIG. 3. Unit cell of LaFeAsO and FeSe. The d xz and d yz orbitalsare illustrated. The x -axis is different from the principal X axis ofthe crystal by ◦ . We take the unit of energy so that the largest hopping inte-gral of LaFeAsO is t d xz ( yz ) d xz ( yz ) = 1 , which is the near-est neighbour hopping integral between parallel d xz ( yz ) or-bitals. The other hopping integrals and temperature are scaledas t/t d xz ( yz ) d xz ( yz ) . The unit cell and the shape of the twoorbitals are illustrated in Fig. 3. Although we may expect alarger hopping integral of red orbitals, for the 10-orbital mod-els of d electrons on Fe ions, the hopping between the blueorbitals is larger than that between the red orbitals, since thehopping through p orbitals of As ions is dominant. Thus, wechoose this hopping integral as the unit of energy. We takethe same unit also for FeSe, because the hopping parameteris not significantly different between the iron-based supercon-ductors. Note that the axes for orbitals in Fig. 3 are rotatedby ◦ from the principal X and Y axes in accordance withthe conventional notation. We adopt the ( x, y ) axes and corre-sponding wave vector ( k x , k y ) for calculating the thermody-namic EQMs, although the Fermi surfaces are drawn by usingthe ( k X , k Y ) axes. In this notation a finite EQM Q x − y isinduced by the orbital order of d xz and d yz orbitals. -0.002-0.001 0 0.001 0.002 0.003 0.004 0.005 0.006 0 0.02 0.04 0.06 0.08 0.1 Q x - y T fermi-sea termfermi-surface term Q x -y FIG. 4. Temperature dependence of the thermodynamic EQM Q x − y in the tetragonal LaFeAsO. We show the contributions fromthe Fermi-sea term and the Fermi-surface term by circles with blueline and squares with red line, respectively. The total EQM is shownby triangles with yellow line. -0.004-0.002 0 0.002 0.004 0.006 0.008 0.01 0.012 0 0.02 0.04 0.06 0.08 0.1 Q x - y T fermi-sea termfermi-surface term Q x -y FIG. 5. Thermodynamic EQM Q x − y in the orthorhombicLaFeAsO. The temperature dependences of the thermodynamic EQMsin the tetragonal and orthorhombic LaFeAsO are shown inFig. 4 and Fig. 5, respectively. We show the contributionsfrom the Fermi-sea term and the Fermi-surface term by bluelines and red lines, respectively. Total thermodynamic EQMsare shown by yellow lines. The quantities are indicated by thesame color and symbols in all later results.In both tetragonal and orthorhombic structures, the Fermi-sea term is dominant for the EQMs. Thus, the EQM in the ne-
FIG. 6. Fermi-sea term of the tetragonal LaFeAsO from each k points. We set T = 0 . . Large contributions from the momen-tum near the Fermi surfaces around the Γ = (0 , and M = ( π, π ) points are observed. matic phase of LaFeAsO mainly has a geometric origin. Com-paring the results for the tetragonal and orthorhombic crystals,we notice the additional contribution to the thermodynamicEQMs from the orthorhombic crystal deformation. We seesizable EQMs at T = 0 . in Fig. 5, where ∆( T ) = 0 .Here we show that the dominant geometric origin of theEQMs is a unique property of LaFeAsO by comparing the10-orbital model from first-principles with a toy model. Twoof us previously calculated the thermodynamic EQMs usinga toy model constructed for only the d xz and d yz orbitals .The Fermi-sea term with a geometric origin is much largerin the 10-orbital model than the toy model. This is becausethe Fermi-sea term is enhanced by the band degeneracy. Tosee this we show the momentum-resolved Fermi-sea term inFig. 6. Dominant contributions come from the k points nearthe Fermi surfaces. Thus, in LaFeAsO, the geometrically non-trivial property of wave functions due to the band degeneracyaround the Fermi surfaces gives rise to the sizable Fermi-seaterm (see also Appendix A). The EQMs from the momen-tum near the Γ and M points are plotted in Fig. 7. Fig-ures 7(a) and 7(b) show the EQM arising from the momentum | k x | , | k y | < π/ and that from | k x − π | , | k y − π | < π/ ,respectively. As shown in Figs. 6 and 7, the main contri-bution comes from the Fermi surfaces around the M point.On the other hand, the degenerate band structure is too sim-plified in the toy model, and in particular, the contributionfrom the momentum near the M point is almost overlooked.While the electron and hole Fermi surfaces are additive for theFermi-sea term, the Fermi-surface term is partially cancelled.Given that the multiple band degeneracy near the Fermi sur-faces is a unique property of LaFeAsO, the geometric originof the EQM is also regarded as a characteristic property ofLaFeAsO. To support this argument, we calculate the chemi-cal potential dependence of the EQMs in the 10-orbital modeland find that the Fermi-surface term is comparable or largerthan the Fermi-sea term in most cases except for the realisticparameter of LaFeAsO (see Appendix B). -0.0015-0.001-0.0005 0 0.0005 0.001 0.0015 0.002 0.0025 0 0.02 0.04 0.06 0.08 0.1 T 𝑄 ! ! " ! 𝑄 ! ! " ! (a) -0.001-0.0005 0 0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 0.0045 0 0.02 0.04 0.06 0.08 0.1 𝑄 ! ! " ! T 𝑄 ! ! " ! (b) FIG. 7. Contributions to the thermodynamic EQM Q x − y inLaFeAsO from momentum spaces (a) − π/ < k x , k y < π/ and(b) π/ < k x , k y < π/ . Not only the Fermi-surface term butalso the Fermi-sea term mainly come from the momentum near theFermi surfaces (see also Fig. 6). Note that the geometrical origin of the EQMs is not aunique consequence of the orbital order. Indeed, the or-thorhombic lattice distortion causes not only the orbital po-larization but also the bond anisotropy, and the magnitude ofthe bond order is larger than that of the orbital order. We seea sizable Fermi-sea term due to the lattice distortion (Fig. 5 at T = 0 . ), which is much larger than the Fermi-surface term. B. FeSe
Conducting a first-principles calculation for FeSe, we con-struct the 10-orbital tight-binding model similar to LaFeAsO.Lattice parameters given in Ref. 5 are adopted, and the spacegroup is P /nmm . It is known that in FeSe tiny Fermi-surfaces obtained by angle-resolved photoemission spec-troscopy (ARPES) measurements are not reproduced bythe first-principle calculation: larger Fermi surfaces and an ex-tra Fermi surface of the d xy orbital appear. To reproduce theexperimentally observed Fermi surfaces of FeSe, we take intoaccount additional hopping parameters in addition to thosegiven by the WIEN
2k code, in a similar manner to Refs. (see Appendix. C for details). The additional hopping pa-rameters may stem from the self-energy correction . Differ-ent from LaFeAsO, sign reversal of the orbital polarization inthe momentum space between the Γ and M points has beenobserved by ARPES for FeSe and studied theoreti-cally . To reproduce this property of the nematic order, wetake into account the molecular fields of the bond order in ad-dition to the orbital order (see Appendix. D for details). Thetotal molecular field is given by Γ = Γ orb + Γ bond , (21) Γ orb = (cid:88) k (cid:88) σ,i =1 , (cid:104) c † d xz iσ ( k ) c d xz iσ ( k ) − c † d yz iσ ( k ) c d yz iσ ( k ) (cid:105) , (22) Γ bond = (cid:88) k (cid:18) cos k y − k x − cos k y + k x (cid:19) × (cid:88) σ,l = xz,yz (cid:104) c † d l σ ( k ) c d l σ ( k ) + c † d l σ ( k ) c d l σ ( k ) (cid:105) . (23)The Fermi surfaces of the 10-orbital model with addi-tional hopping parameters and nematic order parameter ∆( T ) - (cid:1) (cid:1) - (cid:1) (cid:1) k x k y - (cid:1) (cid:1) - (cid:1) (cid:1) k x k y - (cid:1) (cid:1) - (cid:1) (cid:1) k x k y - (cid:1) (cid:1) - (cid:1) (cid:1) k x k y (a) (d)(c) (b) - (cid:1) (cid:1) - (cid:1) (cid:1) k x k y - (cid:1) (cid:1) - (cid:1) (cid:1) k x k y - (cid:1) (cid:1) - (cid:1) (cid:1) k x k y - (cid:1) (cid:1) - (cid:1) (cid:1) k x k y 𝑘 ! 𝑘 " FIG. 8. Fermi surfaces of FeSe at (a) T = 0 . , (b) T = 0 . , (c) T = 0 . , and (d) T = 0 . . (a) is the normal state, while (b)-(d)are the nematic states with finite ∆( T ) . are shown in Fig. 8. The Fermi surfaces are distorted withgrowing the nematic order. In the low temperature region[Fig. 8(d)], one Fermi surface near the Γ point disappears ow-ing to the orbital polarization, consistent with experiments .The disappearance of the Fermi surface is related to thechange in the Fermi-sea term of the EQM that will be shownbelow. Note that the shape of the remaining Fermi surfacesin Fig. 8(d) is slightly different from what observed in the ex-periment since we do not take into account a weak spin-orbitcoupling . For a remark, we need to calculate the chemicalpotential at each temperature T to keep the particle numberand reproduce the disappearance of the Fermi surface. -0.06-0.05-0.04-0.03-0.02-0.01 0 0.01 0.02 0 0.02 0.04 0.06 0.08 0.1 Q x - y T fermi-sea termfermi-surface term Q x -y FIG. 9. Thermodynamic EQM Q x − y in FeSe. The EQM Q x − y in FeSe is shown in Fig. 9. At low tem-peratures, the Fermi-sea term with a geometric origin is alsodominant in FeSe. In contrast to LaFeAsO, the Fermi-seaterm is negative. This contribution mainly comes from theelectronic states near the Γ point as we show Fig. 10. Wesee the negative and dominant contribution to the Fermi-seaterm from near the Γ point and the positive contribution fromnear the M point. As for the Fermi-surface term, Fig. 10 alsoshows that the dominant contribution comes from near the Γ -0.06-0.05-0.04-0.03-0.02-0.01 0 0.01 0.02 0.03 0 0.02 0.04 0.06 0.08 0.1 -0.03-0.010.010-0.020.02-0.04-0.05-0.06 T 𝑄 ! ! " ! 𝑄 ! ! " ! (a) -0.01-0.008-0.006-0.004-0.002 0 0.002 0.004 0.006 0.008 0.01 0.012 0 0.02 0.04 0.06 0.08 0.1 Q x - y T fermi-sea termfermi-surface term Q x -y T 𝑄 ! ! " ! 𝑄 ! ! " ! (b) FIG. 10. Contributions to the thermodynamic EQM Q x − y in FeSefrom (a) − π/ < k x , k y < π/ and (b) π/ < k x , k y < π/ .Both Fermi-sea and Fermi-surface terms mainly originate from themomentum near the Fermi surfaces (see also Fig. 11). point. FIG. 11. Fermi-sea term of the EQM in FeSe from each k points for(a) T = 0 . and (b) T = 0 . . (c) shows an enlarged illustration of(b) near the Γ point. The contribution from near the Fermi surfacesdominates the Fermi-sea term. For more details, the Fermi-sea term contributions fromeach k points at T = 0 . and T = 0 . are shown inFigs. 11(a) and 11(b,c), respectively. The opposite contri-bution from the Γ and M points is revealed, consistent withFig. 10. Furthermore, we see a change in the Fermi-seaterm arising from near the Γ point between T = 0 . and T = 0 . , while that from the M point is almost temperatureindependent in this region. This change is caused by the dis-appearance of a Fermi surface discussed above. In Fig. 11(c),we see a large contribution, which is illustrated by white color,from the momentum around which the Lifshitz transition oc-curs.Finally, we discuss the similarities and differences betweenLaFeAsO and FeSe. From the results, we find that the Fermi-sea term with geometric origin is dominant in the thermody-namic EQM of FeSe as well as of LaFeAsO. This findingimplies that the geometrically nontrivial properties of wavefunctions are ubiquitous in iron-based superconductors. Be-cause the geometric properties are owing to the multi-orbitaland multi-band structure, the band degeneracy naturally playsimportant roles for the EQM as well as for the nematic orderand superconductivity. On the other hand, when we look atthe details, the sign of the EQM is opposite between FeSe andLaFeAsO, and the momentum-resolved EQM shows different structures. Thus we need a precise model taking account ofthe realistic electronic structure for quantifying the EQMs ofnematic phases.We would like to stress the usefulness of the thermody-namic formulation for the EQMs. By the thermodynamicEQMs, the nematic order can be quantified in a unified way,even when not only the electronic structures but also the ne-matic order parameters are different between the materials asin the cases of LaFeAsO and FeSe. C. Cuprate superconductors
To illuminate the unique properties of iron-based supercon-ductors, that is, multi-band structure and resulting geomet-ric origin of the EQM, we here calculate the thermodynamicEQM of cuprate superconductors for a comparison. For thenematic order in cuprate superconductors, we consider the d x − y -wave bond order studied extensively . For a com-parison, the orbital order of O2 p x and O2 p y orbitals is alsostudied later. Evaluation of translation-symmetry-breaking or-der, such as the CDW and PDW order, is left for future studies.For the study of cuprate superconductors, we construct the17-orbital tight-binding model, which consists of the d or-bitals of coppers and the p orbitals of oxygens in a unit cell,using the WIEN
2k and Wannier90 code. The space group is I /mmm and lattice parameters are adopted from Ref. 86at T = 295 K. The tight-binding parameters are derived forthe representative mother compound La CuO . AlthoughLa CuO is a Mott insulator and the nematic order occurs byhole doping , we adopt the half-filling model with n = 5 . ,since the dependence on the carrier density is negligible. Themodel takes into account four oxygen ions and one copper ionin the unit cell. Thus, creation operators of the p electronshave index for the sublattice. The two oxygens are located onthe CuO plane, while the other two are apical oxygens. Weset the unit of energy so that the largest nearest-neighbor d - p hopping t d x − y p x is unity.As for the nematic order parameter, the molecular field ofthe d x − y -wave bond order is given by Γ = (cid:88) k (cid:88) σ (cos k x − cos k y ) c † d x − y σ ( k ) c d x − y σ ( k ) . (24)For later comparison, we also examine the p -orbital orderwhose molecular field is written as Γ = (cid:88) k (cid:88) i =1 , (cid:104) c † p x iσ ( k ) c p x iσ ( k ) − c † p y iσ ( k ) c p y iσ ( k ) (cid:105) . (25)The index i = 1 , indicates the oxygens on the CuO plane.Figure 12 shows distortion of the Fermi surface due to the ne-matic order. It is significant in the bond-ordered state, becausethe electronic states near the Fermi level mainly consist of the d x − y orbital, although it is hybridyzed with the p orbitals.The thermodynamic EQM induced by the bond order isshown in Fig. 13. We see that the Fermi-surface term is domi- - (cid:1) (cid:1) - (cid:1) (cid:1) k x k y - (cid:1) (cid:1) - (cid:1) (cid:1) k x k y - (cid:1) (cid:1) - (cid:1) (cid:1) k x k y (a) (b)(c) FIG. 12. Fermi surface of the -orbital model for La CuO . (a)The normal state at T = 0 . [ ∆( T ) = 0 ]. (b) The d x − y -wavebond-ordered state at T = 0 . . (c) The p -orbital-ordered state at T = 0 . . -0.005 0 0.005 0.01 0.015 0.02 0.025 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 Q x - y T fermi-sea termfermi-surface term Q x -y FIG. 13. Thermodynamic EQM in the -orbital d - p model forLa CuO with the d x − y -wave bond order. nant in contrast to the results for the iron-based superconduc-tors. Unlike the iron-based superconductors, the band near theFermi level is isolated from others, although the hybridized d - p orbital forms the Fermi surface. Comparison betweenthe iron-based and cuprate superconductors implies a uniqueproperty of the former from the viewpoint of the EQM; theband degeneracy near the Fermi level gives rise to geometri-cally nontrivial properties that lead to the dominant Fermi-seaterm of the EQM. We would like to stress that the differencemainly comes from the underlying electronic structure and notfrom the character of nematic order parameters. Indeed, theorbital order also induces the dominant Fermi-surface term incuprate superconductors, as we see in Fig. 14.In Appendix E, we show the qualitatively same results forthe three-orbital d - p model, which has been extensively an-alyzed in the previous studies of cuprate superconductors.Thus, just increasing the number of orbitals does not enhancethe Fermi-sea term. The band degeneracy near the Fermi sur-face is an essential condition for a large Fermi-sea term. Com-parison between the d x − y -wave bond order and the p -orbitalorder shows a larger EQM in the former, although we assume -0.005-0.0045-0.004-0.0035-0.003-0.0025-0.002-0.0015-0.001-0.0005 0 0.0005 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 Q x - y T fermi-sea termfermi-surface term Q x -y FIG. 14. Thermodynamic EQM in the -orbital d - p model forLa CuO with the orbital order. the same energy scale of the order parameters. This is sim-ply because the Fermi-surface term is dominant in both cases,and because the distortion of the Fermi surfaces is small in the p -orbital-ordered state. V. SUMMARY AND DISCUSSION
In this paper, after showing the failure of the EQMs givenby the electromagnetism and Wannier function methods, weevaluated the thermodynamic EQMs in LaFeAsO, FeSe, andLa CuO using the first-principles calculation and assumingthe candidate nematic order parameters. The thermodynamicEQMs have been proposed as one of the fundamental quan-tities characterizing the C -symmetry breaking in the variousnematic phases. From the results, we found that the EQMsin iron-based superconductors have a geometric origin. Thisis due to the highly degenerate band structure near the Fermilevel, unique to iron-based superconductors. In contrast, theEQMs of cuprate superconductors mainly originate from thedistortion of Fermi surfaces. In this case, the magnitude andsign of the EQMs can be derived from the band structure,which can be observed by ARPES for instance: The thirdterm of Eq. (11) gives the EQMs. Thus, differences betweeniron-based superconductors and cuprate superconductors inthe nematic phases were elucidated from the perspective ofthe EQMs.In addition to the conceptual meaning characterizing the C -symmetry breaking, the EQMs are related to some elec-tric responses caused by the symmetry breaking. A thermo-dynamic relation between the EQMs and the electric suscep-tibility ∂Q ij ( µ ) ∂µ = − χ eij , (26)has been proved for insulators at T = 0 . This relation impliesthat the geometric contribution plays an essential role for theelectric susceptibility in the insulating ground state of iron-based superconductors’ mother compounds.We also see indirect relations of the EQMs with some opti-cal responses. For example, the optical attenuation coefficientis given by , ε ij att ( ω ) = iπ (cid:88) n (cid:54) = m (cid:90) d d k (2 π ) d g ijnm ( k ) [ f ( (cid:15) n ( k )) − f ( (cid:15) m ( k ))] × δ ( (cid:15) m ( k ) − (cid:15) n ( k ) − (cid:126) ω ) . (27)Here, g ijnm ( k ) = ( A inm ( k ) A jmn ( k ) + c.c ) is the band-resolved quantum metric, which also appeared in the thermo-dynamic EQMs. The anisotropic optical attenuation causedby the nematic C -symmetry breaking may be related to thegeometric term of the thermodynamic EQMs. As for therelation to the nonlinear optics, the photocurrent responsesin time-reversal-symmetric and P T -symmetric systems havebeen recently classified, and the results reveal that the quan-tum metric is an essential quantity for some photocurrent re-sponses, such as shift current, magnetic injection current, andgyration current . Thus, elucidation of linear and nonlin-ear optical responses in iron-based superconductors may be anintriguing future issue. For the photocurrent generation, thespace inversion symmetry must be broken. Indeed, the inver- sion symmetry is broken in some iron-based superconductors,such as FeSe/SrTiO and heavily-doped LaFeAsO . ACKNOWLEDGMENTS
We thank K. Kimura, A. Shitade, and T. Yamashita forfruitful discussions. This work was supported by KAK-ENHI (Grants No. JP18H05227, No. JP18H01178, and No.JP20H05159) from the Japan Society for the Promotion ofScience (JSPS). This work was supported by SPIRITS 2020of Kyoto University.
Appendix A: Geometric contribution to thermodynamic EQMsenhanced by band degeneracy
In the main text, we have shown that the band degeneracyenhances the geometric contribution to the thermodynamicEQMs. To show this explicitly, we discuss an alternative ex-pression of Eq. (11): Q ij = 12 (cid:90) BZ d d k (2 π ) d (cid:88) mn (cid:2) V inm ( k ) V jmn ( k ) + c.c. (cid:3) (cid:20) − δ nm { (cid:15) n ( k ) − (cid:15) m ( k ) } F nm ( k ) + δ nm f (cid:48)(cid:48) ( (cid:15) m ( k )) (cid:21) , (A1) F nm ( k ) = f ( (cid:15) n ( k )) + f ( (cid:15) m ( k ))2 − (cid:15) n ( k ) − (cid:15) m ( k ) (cid:90) (cid:15) n ( k ) (cid:15) m ( k ) d(cid:15)f ( (cid:15) ) , (A2)where V inm ( k ) = (cid:104) u n ( k ) | ∂ k i H ( k ) | u m ( k ) (cid:105) . Since F nm ( k ) = O (1) , the contribution from the non-degenerate bands is suppressed by the factor { (cid:15) n ( k ) − (cid:15) m ( k ) } − . Fornearly degenerate bands, | (cid:15) n ( k ) − (cid:15) m ( k ) | (cid:28) T , Eq. (A2) isapproximated as F nm ( k ) = f ( (cid:15) n ( k )) + f ( (cid:15) m ( k ))2 − β ln(1 + e β(cid:15) n ( k ) ) − ln(1 + e β(cid:15) m ( k ) ) (cid:15) n ( k ) − (cid:15) m ( k ) ≈ f (cid:48)(cid:48) ( (cid:15) m ( k )) { (cid:15) n ( k ) − (cid:15) m ( k ) } . (A3)Because the expression contains the factor f (cid:48)(cid:48) ( (cid:15) m ( k )) , we no-tice that a large contribution is given by the momentum nearthe Fermi surface. From these discussions, we understand thatthe geometric contributions to the EQMs are enhanced by theband degeneracy near the Fermi level. When all the bands arenearly degenerate, we have Q ij = 12 (cid:90) BZ d d k (2 π ) d (cid:88) mn (cid:2) V inm ( k ) V jmn ( k ) + c.c. (cid:3) × (cid:20) (1 − δ nm ) f (cid:48)(cid:48) ( (cid:15) m ( k ))12 + δ nm f (cid:48)(cid:48) ( (cid:15) m ( k )) (cid:21) . (A4)In this case, the Fermi-sea term naturally has a similar form tothe Fermi-surface term. In the realistic situation, the magni- tude of the Fermi-sea term depends on the details of the elec-tronic structure. Appendix B: Chemical potential dependence of the EQM in the10-orbital model for LaFeAsO
Here we show the thermodynamic EQMs for various chem-ical potentials in the model of tetragonal LaFeAsO. Althoughin the main text the chemical potential is determined so thatthe particle number is 6.0, we change the chemical poten-tial with keeping the other parameters. Figure 15 is the re-sult at T = 0 . . It is shown that the geometric contri-bution is dominant only in a small part of the parameterrange. Thus, we consider that the dominant geometric con-0 -0.04-0.03-0.02-0.01 0 0.01 0.02-6 -4 -2 0 2 4 6 8 Q x - y µ fermi-sea termfermi-surface term Q x -y FIG. 15. Chemical potential dependence of the thermodynamicEQM in the model for the tetragonal LaFeAsO. We set T = 0 . . tribution is a unique property of the iron-based superconduc-tors for particle numbers around 6.0. Actually, the nematicorder of LaFeAsO − x F x has been observed in the region, . ≤ n ≤ . . We confirmed that the geometric term isdominant in this region. Appendix C: Tight-binding model reproducing Fermi surfacesof FeSe
To reproduce the Fermi surfaces of FeSe observed in exper-iments, we slightly modify the hopping parameters given bythe first-principles calculation . For this purpose, the en-ergies of the d xy -orbital band and the d xz/yz -orbital band areshifted by ( − . , , . ) and ( − . , , . ) at ( Γ , X , M )points in the folded Brillouin zone, respectively. For this en-ergy shift, the hopping parameters are changed so as to satisfy δE l (Γ) = δt on − site ll + 4 δt nn ll + 4 δt nnn ll , (C1) δE l ( X ) = δt on − site ll , (C2) δE l ( M ) = δt on − site ll − δt nnn ll , (C3)where we represent the energy shifts of the l -orbital band at Γ , X and M points as δE l (Γ) , δE l ( X ) and δE l ( M ) , respec-tively. The modification in the intra-orbital hopping integralis represented by δt ll , and ”on-sine”, ”nn”, and ”nnn” denotethe on-sine, first nearest neighbour, and second nearest neigh-bour hoppings, respectively. In the 10-orbital model with twosublattices in the unit cell, δt nn ll ( δt nnn ll ) is the inter-sublattice(intra-sublattice) hopping. We also tune the chemical poten-tial to keep the filling n = 6 . Using these parameters, weobtain the Fermi surfaces in Fig. 8(a). Appendix D: Sign-reversing orbital polarization in FeSe
ARPES measurements clarified sign-reversing orbital po-larization in FeSe , different from LaFeAsO. Thus,we introduce the molecular field, Eqs. (21)-(23), yielding theorbital polarization with opposite sign between the Γ and M points. To understand the sign reversal, we here consider theunfolded BZ, for simplicity. By Eqs. (21)-(23), the molecular field gives the momentum-dependent energy shift of the d xz and d yz orbitals as δE nem d xz ( k ) = 2∆( T ) (cid:18) cos k y − cos k x + 12 (cid:19) , (D1) δE nem d yz ( k ) = 2∆( T ) (cid:18) cos k y − cos k x − (cid:19) . (D2)The momentum dependence in the energy shift is shown inFig. 16, which resembles a theoretical result for the sign-reversing orbital polarization . In the folded Brillouin zone, k = ( π, and (0 , π ) are equivalent ( M point). Therefore,the energy splitting between the orbitals is δE nem d xz (0 , π ) − δE nem d yz ( π, at the M point, while it is δE nem d xz (0 , − δE nem d yz (0 , at the Γ point. As shown in Fig. 16, the signis opposite between the Γ and M points, consistent with thesing-reversing orbital polarization in FeSe . FIG. 16. Momentum dependence of the energy shift due to themolecular field for FeSe, Eqs. (21)-(23), in the unfolded BZ. (a) and(b) show δE nem d xz ( k ) and δE nem d yz ( k ) , respectivley. Red, blue, andwhite represent positive, negative, and zero value, respectively. Theblack lines show the BZ of the two-sublattice model. Γ and M arethe points of the two-sublattice model. Appendix E: EQM in three-orbital d - p model for cupratesuperconductors Here, we show the EQM in the 3-orbital d - p model whichhas been studied for cuprate superconductors . Themodel takes into account the d x − y -orbital of coppers andthe p x and p y orbitals of oxygens. For comparison with the -orbital model studied in the main text, we assume the halffiling. The hopping parameters and the molecular field arethe same as those in the 17-orbital model. In the cuprates,there is no band degeneracy near the Fermi surface, and thelow-energy electron states are appropriately described by the3-orbital d - p model. Thus, when the EQMs are mainly givenby the Fermi-surface term, we expect qualitatively the sameresults as the -orbital model.1 -0.004-0.002 0 0.002 0.004 0.006 0.008 0.01 0.012 0 0.02 0.04 0.06 0.08 0.1 Q x - y T fermi-sea termfermi-surface term Q x -y FIG. 17. Thermodynamic EQMs in the 3-orbital d - p model forcuprate superconductors with bond order. -0.0025-0.002-0.0015-0.001-0.0005 0 0.0005 0 0.02 0.04 0.06 0.08 0.1 Q x - y T fermi-sea termfermi-surface term Q x -y FIG. 18. Thermodynamic EQMs in the 3-orbital d - p model forcuprate superconductors with orbital order. Indeed, the thermodynamic EQMs show the similar behav-iors to those in the -orbital model, as shown in Figs. 17and 18. In both cases of the d x − y -wave bond order and the p -orbital order, the Fermi-surface term is dominant. The mag-nitude of the EQMs is larger in the bond-ordered state than theorbital-ordered state, like in the -orbital model. On the otherhand, we see differences in the Fermi-sea term between the -orbital and -orbital d - p models: even the sign is opposite inthe orbital-ordered state. This implies the importance of therealistic multi-orbital model for the evaluation of Fermi-seaterms. ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] C. de la Cruz, Q. Huang, J. W. Lynn, J. Li, W. R. II, J. L. Zarestky,H. A. Mook, G. F. Chen, J. L. Luo, N. L. Wang, and P. Dai,Nature , 899 (2008). T. Nomura, S. W. Kim, Y. Kamihara, M. Hirano, P. V. Sushko,K. Kato, M. Takata, A. L. Shluger, and H. Hosono, Supercon-ductor Science and Technology , 125028 (2008). F.-C. Hsu, J.-Y. Luo, K.-W. Yeh, T.-K. Chen, T.-W. Huang, P. M. Wu, Y.-C. Lee, Y.-L. Huang, Y.-Y. Chu, D.-C. Yan, and M.-K. Wu, Proceedings ofthe National Academy of Sciences S. Margadonna, Y. Takabayashi, M. T. McDonald,K. Kasperkiewicz, Y. Mizuguchi, Y. Takano, A. N. Fitch,E. Suard, and K. Prassides, Chem. Commun. , 5607 (2008). A. E. B¨ohmer, F. Hardy, F. Eilers, D. Ernst, P. Adelmann,P. Schweiss, T. Wolf, and C. Meingast, Phys. Rev. B , 180505(2013). M. Rotter, M. Tegel, D. Johrendt, I. Schellenberg, W. Hermes,and R. P¨ottgen, Phys. Rev. B , 020503 (2008). Q. Huang, Y. Qiu, W. Bao, M. A. Green, J. W. Lynn, Y. C. Gas-parovic, T. Wu, G. Wu, and X. H. Chen, Phys. Rev. Lett. ,257003 (2008). K. Ishida, Y. Nakai, and H. Hosono, Journal ofthe Physical Society of Japan , 062001 (2009),https://doi.org/10.1143/JPSJ.78.062001. G. R. Stewart, Rev. Mod. Phys. , 1589 (2011). P. Dai, Rev. Mod. Phys. , 855 (2015). T. Shibauchi, T. Hanaguri, and Y. Matsuda, Journalof the Physical Society of Japan , 102002 (2020),https://doi.org/10.7566/JPSJ.89.102002. I. I. Mazin, D. J. Singh, M. D. Johannes, and M. H. Du, Phys.Rev. Lett. , 057003 (2008). K. Kuroki, S. Onari, R. Arita, H. Usui, Y. Tanaka, H. Kontani,and H. Aoki, Phys. Rev. Lett. , 087004 (2008). H. Ikeda, Journal of the Physical Society of Japan , 123707(2008), https://doi.org/10.1143/JPSJ.77.123707. X. Wang, Q. Liu, Y. Lv, W. Gao, L. Yang, R. Yu, F. Li, andC. Jin, Solid State Communications , 538 (2008). Y. Yanagi, Y. Yamakawa, N. Adachi, and Y. ¯Ono, Jour-nal of the Physical Society of Japan , 123707 (2010),https://doi.org/10.1143/JPSJ.79.123707. R. Thomale, C. Platt, W. Hanke, J. Hu, and B. A. Bernevig, Phys.Rev. Lett. , 117001 (2011). H. Kontani, T. Saito, and S. Onari, Phys. Rev. B , 024528(2011). S. Onari and H. Kontani, Phys. Rev. Lett. , 137001 (2012). R. M. Fernandes, L. H. VanBebber, S. Bhattacharya, P. Chandra,V. Keppens, D. Mandrus, M. A. McGuire, B. C. Sales, A. S.Sefat, and J. Schmalian, Phys. Rev. Lett. , 157003 (2010). R. M. Fernandes and J. Schmalian, Superconductor Science andTechnology , 084005 (2012). T. Yamada, J. Ishizuka, and Y. ¯Ono, J. Phys. Soc. Jpn. , 043704(2014). S. Mukherjee, A. Kreisel, P. J. Hirschfeld, and B. M. Andersen,Phys. Rev. Lett. , 026402 (2015). R. Daou, J. Chang, D. LeBoeuf, O. Cyr-Choini`ere, F. Lalibert´e,N. Doiron-Leyraud, B. J. Ramshaw, R. Liang, D. A. Bonn, W. N.Hardy, and L. Taillefer, Nature , 519 (2010). Y. Sato, S. Kasahara, H. Murayama, Y. Kasahara, E.-G. Moon,T. Nishizaki, T. Loew, J. Porras, B. Keimer, T. Shibauchi, andY. Matsuda, Nature Physics , 1074 (2017). H. Yamase and H. Kohno, Journal of the Physical Society ofJapan , 2151 (2000), https://doi.org/10.1143/JPSJ.69.2151. C. J. Halboth and W. Metzner, Phys. Rev. Lett. , 5162 (2000). C. Honerkamp, M. Salmhofer, N. Furukawa, and T. M. Rice,Phys. Rev. B , 035109 (2001). I. Khavkine, C.-H. Chung, V. Oganesyan, and H.-Y. Kee, Phys.Rev. B , 155110 (2004). E. Berg, E. Fradkin, S. A. Kivelson, and J. M. Tranquada, NewJournal of Physics , 115004 (2009). W. Metzner, M. Salmhofer, C. Honerkamp, V. Meden, andK. Sch¨onhammer, Rev. Mod. Phys. , 299 (2012). S. Bulut, W. A. Atkinson, and A. P. Kampf, Phys. Rev. B ,155132 (2013). S. Sachdev and R. La Placa, Phys. Rev. Lett. , 027202 (2013). Y. Wang and A. Chubukov, Phys. Rev. B , 035149 (2014). Y. Yamakawa and H. Kontani, Phys. Rev. Lett. , 257001(2015). K. Kawaguchi, Y. Yamakawa, M. Tsuchiizu, and H. Kontani,Journal of the Physical Society of Japan , 063707 (2017),https://doi.org/10.7566/JPSJ.86.063707. M. Tsuchiizu, K. Kawaguchi, Y. Yamakawa, and H. Kontani,Phys. Rev. B , 165131 (2018). I. Affleck and J. B. Marston, Phys. Rev. B , 3774 (1988). H. J. Schulz, Phys. Rev. B , 2940 (1989). F. C. Zhang, Phys. Rev. Lett. , 974 (1990). C. M. Varma, Phys. Rev. B , 14554 (1997). P. A. Lee, Phys. Rev. X , 031017 (2014). D. F. Agterberg, J. S. Davis, S. D. Edkins, E. Fradkin, D. J.Van Harlingen, S. A. Kivelson, P. A. Lee, L. Radzihovsky,J. M. Tranquada, and Y. Wang, Annual Review of CondensedMatter Physics , 231 (2020), https://doi.org/10.1146/annurev-conmatphys-031119-050711. J. D. Jackson, Classical Electrodynamics, 3rd Edition (1998). D. Vanderbilt, “Berry phases in electronic structure theory: Elec-tric polarization, orbital magnetization and topological insula-tors,” (Cambridge University Press, 2018). R. D. King-Smith and D. Vanderbilt, Phys. Rev. B , 1651(1993). D. Vanderbilt and R. D. King-Smith, Phys. Rev. B , 4442(1993). N. Marzari and D. Vanderbilt, Phys. Rev. B , 12847 (1997). N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vander-bilt, Rev. Mod. Phys. , 1419 (2012). W. A. Benalcazar, B. A. Bernevig, andT. L. Hughes, Science , 61 (2017),https://science.sciencemag.org/content/357/6346/61.full.pdf. W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Phys. Rev.B , 245115 (2017). L. He, Z. Addison, E. J. Mele, and B. Zhen, Nature Communi-cations , 3119 (2020). M. Ezawa, Phys. Rev. Lett. , 026801 (2018). S. Imhof, C. Berger, F. Bayer, J. Brehm, L. W. Molenkamp,T. Kiessling, F. Schindler, C. H. Lee, M. Greiter, T. Neupert, andR. Thomale, Nature Physics , 925 (2018). M. Serra-Garcia, V. Peri, R. S¨usstrunk, O. R. Bilal, T. Larsen,L. G. Villanueva, and S. D. Huber, Nature , 342 (2018). Z. Song, Z. Fang, and C. Fang, Phys. Rev. Lett. , 246402(2017). S. Franca, J. van den Brink, and I. C. Fulga, Phys. Rev. B ,201114 (2018). M. Hirayama, R. Takahashi, S. Matsuishi, H. Hosono, and S. Murakami, Phys. Rev. Research , 043131 (2020). C. W. Peterson, W. A. Benalcazar, T. L. Hughes, and G. Bahl,Nature , 346 (2018). H. Watanabe and S. Ono, Phys. Rev. B , 165120 (2020). F. Schindler, M. Brzezi´nska, W. A. Benalcazar, M. Iraola,A. Bouhon, S. S. Tsirkin, M. G. Vergniory, and T. Neupert, Phys.Rev. Research , 033074 (2019). W. A. Benalcazar, T. Li, and T. L. Hughes, Phys. Rev. B ,245151 (2019). A. Shitade, H. Watanabe, and Y. Yanase, Phys. Rev. B ,020407 (2018). A. Shitade, A. Daido, and Y. Yanase, Phys. Rev. B , 024404(2019). Y. Gao, D. Vanderbilt, and D. Xiao, Phys. Rev. B , 134423(2018). Y. Gao and D. Xiao, Phys. Rev. B , 060402 (2018). A. Daido, A. Shitade, and Y. Yanase, Phys. Rev. B , 235149(2020). S. Onari, Y. Yamakawa, and H. Kontani, Phys. Rev. Lett. ,227001 (2016). Resta, R., Eur. Phys. J. B , 121 (2011). J. P. Provost and G. Vallee, Communications in MathematicalPhysics , 289 (1980). D. Xiao, M.-C. Chang, and Q. Niu, Rev. Mod. Phys. , 1959(2010). P. Blaha, K. Schwarz, G. Madsen, D. Kvasnicka, J. Luitz,R. Laskowsk, F. Tran, L. Marks, and L. Marks, “WIEN2k: AnAugmented Plane Wave Plus Local Orbitals Program for Calcu-lating Crystal Properties,” (Techn. Universitat, 2019). I. Souza, N. Marzari, and D. Vanderbilt, Phys. Rev. B , 035109(2001). A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt,and N. Marzari, Computer Physics Communications , 685(2008). S. Raghu, X.-L. Qi, C.-X. Liu, D. J. Scalapino, and S.-C. Zhang,Phys. Rev. B , 220503 (2008). J. Maletz, V. B. Zabolotnyy, D. V. Evtushinsky, S. Thirupathaiah,A. U. B. Wolter, L. Harnagea, A. N. Yaresko, A. N. Vasiliev,D. A. Chareev, A. E. B¨ohmer, F. Hardy, T. Wolf, C. Meingast,E. D. L. Rienks, B. B¨uchner, and S. V. Borisenko, Phys. Rev. B , 220506 (2014). Y. Zhang, M. Yi, Z.-K. Liu, W. Li, J. J. Lee, R. G. Moore,M. Hashimoto, M. Nakajima, H. Eisaki, S.-K. Mo, Z. Hussain,T. P. Devereaux, Z.-X. Shen, and D. H. Lu, Phys. Rev. B ,115153 (2016). Y. Yamakawa, S. Onari, and H. Kontani, Phys. Rev. X , 021032(2016). J. Ishizuka, T. Yamada, Y. Yanagi, and Y. ¯Ono, Jour-nal of the Physical Society of Japan , 014705 (2018),https://doi.org/10.7566/JPSJ.87.014705. T. Gorni, P. Villar Arribi, M. Casula, and L. de’ Medici, arXive-prints , arXiv:2101.01692 (2021), arXiv:2101.01692 [cond-mat.str-el]. K. Nakayama, Y. Miyata, G. N. Phan, T. Sato, Y. Tanabe,T. Urata, K. Tanigaki, and T. Takahashi, Phys. Rev. Lett. ,237001 (2014). T. Shimojima, Y. Suzuki, T. Sonobe, A. Nakamura, M. Sakano,J. Omachi, K. Yoshioka, M. Kuwata-Gonokami, K. Ono, H. Ku-migashira, A. E. B¨ohmer, F. Hardy, T. Wolf, C. Meingast, H. v.L¨ohneysen, H. Ikeda, and K. Ishizaka, Phys. Rev. B , 121111(2014). Y. Suzuki, T. Shimojima, T. Sonobe, A. Nakamura, M. Sakano,H. Tsuji, J. Omachi, K. Yoshioka, M. Kuwata-Gonokami,T. Watashige, R. Kobayashi, S. Kasahara, T. Shibauchi, Y. Mat- suda, Y. Yamakawa, H. Kontani, and K. Ishizaka, Phys. Rev. B , 205117 (2015). M. D. Watson, T. K. Kim, A. A. Haghighirad, N. R. Davies,A. McCollam, A. Narayanan, S. F. Blake, Y. L. Chen, S. Ghan-nadzadeh, A. J. Schofield, M. Hoesch, C. Meingast, T. Wolf, andA. I. Coldea, Phys. Rev. B , 155106 (2015). P. Zhang, T. Qian, P. Richard, X. P. Wang, H. Miao, B. Q. Lv,B. B. Fu, T. Wolf, C. Meingast, X. X. Wu, Z. Q. Wang, J. P. Hu,and H. Ding, Phys. Rev. B , 214503 (2015). J. D. Jorgensen, H. B. Sch¨uttler, D. G. Hinks, D. W. Capone II,K. Zhang, M. B. Brodsky, and D. J. Scalapino, Phys. Rev. Lett. , 1024 (1987). J. E. Sipe and A. I. Shkrebtii, Phys. Rev. B , 5337 (2000). J. Iba˜nez Azpiroz, S. S. Tsirkin, and I. Souza, Phys. Rev. B ,245143 (2018). F. Nastos, J. Rioux, M. Strimas-Mackey, B. S. Mendoza, andJ. E. Sipe, Phys. Rev. B , 205113 (2007). H. Watanabe and Y. Yanase, Phys. Rev. X , 011001 (2021). J. Ahn, G.-Y. Guo, and N. Nagaosa, Phys. Rev. X , 041041(2020). Q.-Y. Wang, Z. Li, W.-H. Zhang, Z.-C. Zhang, J.-S. Zhang, W. Li,H. Ding, Y.-B. Ou, P. Deng, K. Chang, J. Wen, C.-L. Song, K. He,J.-F. Jia, S.-H. Ji, Y.-Y. Wang, L.-L. Wang, X. Chen, X.-C. Ma,and Q.-K. Xue, Chinese Physics Letters , 037402 (2012). D. Liu, W. Zhang, D. Mou, J. He, Y.-B. Ou, Q.-Y. Wang, Z. Li,L. Wang, L. Zhao, S. He, Y. Peng, X. Liu, C. Chen, L. Yu, G. Liu, X. Dong, J. Zhang, C. Chen, Z. Xu, J. Hu, X. Chen, X. Ma,Q. Xue, and X. J. Zhou, Nature Communications , 931 (2012). S. Tan, Y. Zhang, M. Xia, Z. Ye, F. Chen, X. Xie, R. Peng, D. Xu,Q. Fan, H. Xu, J. Jiang, T. Zhang, X. Lai, T. Xiang, J. Hu, B. Xie,and D. Feng, Nature Materials , 634 (2013). S. He, J. He, W. Zhang, L. Zhao, D. Liu, X. Liu, D. Mou, Y.-B.Ou, Q.-Y. Wang, Z. Li, L. Wang, Y. Peng, Y. Liu, C. Chen, L. Yu,G. Liu, X. Dong, J. Zhang, C. Chen, Z. Xu, X. Chen, X. Ma,Q. Xue, and X. J. Zhou, Nature Materials , 605 (2013). J.-F. Ge, Z.-L. Liu, C. Liu, C.-L. Gao, D. Qian, Q.-K. Xue, Y. Liu,and J.-F. Jia, Nature Materials , 285 (2015). Y. Miyata, K. Nakayama, K. Sugawara, T. Sato, and T. Taka-hashi, Nature Materials , 775 (2015). J. Shiogai, Y. Ito, T. Mitsuhashi, T. Nojima, and A. Tsukazaki,Nature Physics , 42 (2016). M. Hiraishi, S. Iimura, K. M. Kojima, J. Yamaura, H. Hiraka,K. Ikeda, P. Miao, Y. Ishikawa, S. Torii, M. Miyazaki, I. Ya-mauchi, A. Koda, K. Ishii, M. Yoshida, J. Mizuki, R. Kadono,R. Kumai, T. Kamiyama, T. Otomo, Y. Murakami, S. Matsuishi,and H. Hosono, Nature Physics , 300 (2014). J. Luo and N. E. Bickers, Phys. Rev. B , 12153 (1993). S. Koikegami, S. Fujimoto, and K. Yamada, Jour-nal of the Physical Society of Japan , 1438 (1997),https://doi.org/10.1143/JPSJ.66.1438. T. Takimoto and T. Moriya, Journal of the Physical Society ofJapan66