Pair

Bloch electrons in multiorbital systems carry quantum geometric information characteristic of their wavevector-dependent interorbital mixing. The geometric nature impacts electromagnetic responses, and this eﬀect carries over to the superconducting state, which receives a geometric contribution to the superﬂuid weight. In this paper, we show that this contribution could become negative under certain appropriate circumstances. This may facilitate the stabilization of Cooper pairings with real space phase modulation, i.e

Introduction.-In systems with multiple orbital degrees of freedom, electrons on the individual Bloch bands are not featureless independent particles. Instead, the motion of a Bloch electron on one band is inherently connected to that of other Bloch bands at the same wavevector, despite them being distinct energy eigenstates. More specifically, expressed in the band basis under which the kinetic part of the Hamiltonian is diagonalized, the velocity operator possesses finite off-diagonal terms, i.e. the interband velocity, V mn µk = (ǫ nk − ǫ mk ) ψ mk |∂ kµ ψ nk , m = n, where ǫ mk is the energy dispersion of the m-th band with eigenvector |ψ mk . The object i ψ mk |k∂ kµ ψ nk depicts a non-Abelian Berry connection between the Bloch states and therefore characterizes their unique quantum geometric properties. Some aspects of the geometric effect on quantum transport have been known for long. The most classic example is the quantum Hall or Chern insulator, in which the quantized Hall conductance is intimately tied to the Berry curvature of the Bloch bands [1][2][3]. Recently, more aspects of the geometry-induced electromagnetic responses have been discussed at length [4][5][6][7]. The geometric nature carries over to the superconducting state. Hence its footprint must also be found in the superconducting electromagnetic responses. The past several years have witnessed considerable attention on the geometry-induced finite superfluid weight in flat band superconductors [8][9][10][11][12][13][14][15][16][17], where conventional theory would have otherwise predicted vanishing superfluid density and hence unsustainable superconductivity. These studies are of particular relevance to the putative (near) flat band superconductivity reported in twisted bilayer graphene [18,19]. Nonetheless, the geometric effect is not unique to flat band systems. Recently, the discus-sion has been extended to geometry-induced effects, including some peculiar optical anomalies, in a more broad spectrum of multiorbital superconductors [20][21][22][23].
In this paper, we demonstrate that the Bloch quantum geometry could also facilitate the formation of novel phases of matter in multiorbital superconductors. Our discussion is motivated by the observation that the geometry-related contribution to the superfluid weight is not necessarily positive definite. To be more concrete, we show that this contribution contains a part that relates to an effective interband Josephson coupling between the bands. This coupling could become negative if the superconducting order parameters on the multiple bands condense into an appropriate configuration, resulting in a suppressed phase stiffness. More strikingly, in the narrow-or flat-band limit where the pairing gap could become as large as a significant fraction of the bandwidth, the geometric contribution could increase significantly, even to the point of returning an overall negative superfluid weight.
Since the superfluid weight characterizes the realspace superconducting phase stiffness, an intriguing question then arises: does the geometry-induced suppression of superfluid weight provide a natural mechanism for the formation of superconducting states with real-space phase mudulation, such as the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) [24,25] or other more general forms of pair density wave (PDW) order? We answer this question in the affirmative, on the basis of model mean-field Bogoliubov de-Gennes (BdG) calculations which explicitly demonstrate the preference for a PDW rather than the uniform phase. Our study therefore provides an intrinsic mechanism for PDW phases in the absence of external magnetic field.
Negative interband Josephson coupling.-We set the stage of our discussions by first deriving the superfluid weight of a weak-coupling two-dimensional two-orbital model with uniform intraorbital spin-singlet Cooper pairing and without spin-orbit coupling. The BdG Hamiltonian matrix of the model, expressed in appropriate orbital basis, can be written as, Here,k = −k,Ĥ 0k and∆ k respectively describe the kinetic and pairing parts of the Hamiltonian. Following the standard linear response theory, the superfluid weight tensor can be expressed as [10], where µ, ν = x, y, ω n = (2n+1)π/β denotes the fermionic Matsubara frequency, andĜ k (iω n ) = (iω n −Ĥ BdG k ) −1 is the Gor'kov Green's function. The velocity operatorV µk takes the following form, where τ z is the third component of the Pauli matrix operating in the particle-hole space of the BdG Hamiltonian (2). The second term in (3) represents the paramagnetic contribution, which vanishes in single-orbital models at zero temperature. The first term, T µν , is the diamagnetic contribution in which the inverse mass tensor is given by, Utilizing (5) and the periodic boundary condition in the Brillouin zone of a lattice model, the diamagnetic term can be written in a different form, In the second line of the above equation, we have taken the approximation, ∂ kµĤ BdG k ∆→0 = τ zVµk , which is valid in the weak-coupling limit as terms with ∂ kµ∆k are expected to be much smaller.
To facilitate our subsequent derivations, we turn to the band basis representation, using a unitary transformation that diagonalizes the kinetic part of (2), i.e.
With a further assumption that the weak Cooper pairing takes place only between electrons on the same band, i.e. intraband pairing, the BdG

Hamiltonian in the basis
Here, ǫ mk and ∆ mk give, respectively, the normal state dispersion and the intraband pairing function of band-m, where m = 1, 2 label the band indices. For future reference, we write down the resultant Bogoliubov dispersion E mk = ǫ 2 mk + |∆ mk | 2 with corresponding eigenstate |mk = (u mk , v mk ) T = (∆ mk , E mk − ǫ mk ) T /N mk and its particle-hole symmetric state |mk = (−v * mk , u * mk ) T , where N mk is a normalization factor. The velocity operator in the same basis can be obtained following the same unitary transformation, where the diagonal elements are simply the band group velocities V mm while the off-diagonal elements, i.e. the interband velocities, as given in (1) [10,20], satisfy V 12 µk = (V 21 µk ) * . Finally, the superfluid weight (3), having the diamagnetic term substituted by (6), is obtained by replacing the velocity operator with (8) and the Green's function withG k (iω n ) = (iω n −H BdG The procedure is equivalent to inserting identity operatorŝ I =Û −1 kÛ k into (3) and then dropping the interband pairings. Taking the xx-component as an explicit example, the zero-temperature superfluid weight reads as follows, As is shown in more detail in the Supplementary, there are multiple terms from the diamagnetic and paramagnetic responses that exactly cancel each other. Hence the two terms on the rhs of this equation are not in one to one correspondence with the two respective responses. Nonetheless, the first term, generated by the virtual intraband transitions |mk ↔ |mk , originates solely from the diamagnetic contribution. In the weak-coupling and isolated-band limit, it returns the conventional expression for the zero-temperature superfluid weight. On the other hand, the second term receives equal contributions from the diamagnetic and paramagnetic responses. Absent in traditional theories, this term is purely of quantum geometric origin. In particular, it is associated exclusively with virtual interband transitions such as |2k ↔ |1k , thereby acquiring an explicit dependence on the interband velocity. Note that, in this approximation the geometric term is finite only when both bands develop Cooper pairing. In fact, it can be cast in a more suggestive form where · · · denotes the expectation value of the ground state. It is thus indicative of an effective interband Josephson coupling, where the interband velocity serves to 'tunnel' Cooper pairs from one band to another. Written more explicitly, this term reads, An interesting observation is that, unlike the first term in (3) which is positive definite, this term could become negative in certain scenarios. In a later illustrative lattice model, V 12 µk is odd in k, and it could be real under a proper gauge choice. Taking a simple example where ∆ mk ≡ ∆ m , the geometric correction (10) is thus negative if sign[∆ 1 ] = −sign[∆ 2 ]. In this case, the total superfluid weight is reduced from the approximation based on conventional isolated-band considerations. Conversely, the superfluid weight is enhanced if the two gaps are of the same sign. Note that the sign of the geometric term must be analyzed on a case by case basis, because V 12 xk could be even in k and ∆ mk may acquire wavevector dependence in different models.
The nontrivial geometric contribution to the superfluid weight has been discussed in quite a number of previous studies [8][9][10][11][12][13][14][15][16][17]. In general, the correction increases with increasing pairing strength [9]. However, in weakcoupling limit the coherence factor |u mk v mk | is peaked around the Fermi wavevectors of band-m and vanishes elsewhere. Hence the term (10) is negligible, except in rare scenarios where the two Fermi surfaces overlap. A bold conjecture is that a strong geometric contribution could be obtained in a model with large ∆/E F , which would see a significantly broadened distribution of |u mk v mk |. However, in this case the approximation made in (6) is no longer valid, and interband pairing may also come into play. Below, we turn to an orbital-basis analysis and extend to the regime where ∆ is comparable to the Fermi energy or the bandwidth. We shall still employ the BdG formalism and the mean-field linearresponse theory for a first-stage study in search for some preliminary qualitative understanding.
Strong geometry-induced suppression of D.-The possibility of negative superfluid weight has in fact been noted in a previous study on twisted bilayer graphene [17], although the implication for PDW order was not mentioned there. Here, we explicitly demonstrate negative superfluid weight in a simpler and more mundane two-orbital model consisting of an s-and a d xyorbital on a square lattice, although the conclusion applies to more general models. Assuming only onsite intraorbital spin-singlet pairing on the two orbitals, the BdG Hamiltonian matrix can be written in the basis (c sk↑ , c dk↑ , c † sk↓ , c † dk↓ ) T as, where ξ ak = −2t a (cos k x + cos k y ) − µ a (a = s, d) represent the intraorbital dispersion and λ k = 4t ′ sin k x sin k y denotes the interorbital mixing. In practice, we keep a balance between t ′ and t s − t d , µ s − µ d , so as to ensure a sizable interband velocity. At this stage, we assume |∆ s | = |∆ d | = ∆ which, when transformed into the band basis, happens to determine the intraband pairing order parameter configuration on the two bands (see the Supplementary). Specifically, the ++ configuration witĥ ∆ ++ = (∆ s , ∆ d ) = (∆, ∆) leads to ∆ 1k = ∆ 2k , while the +− configuration with∆ +− = (∆, −∆) gives rise to ∆ 1k = −∆ 2k . Hence, according to (10), these two distinct uniform phases likely carry opposite geometric corrections to the superfluid weight. The superfluid weight obtained from (3) for a representative given band structure is shown in Fig. 1. Consistent with (10), a dichotomy between the∆ ++ and∆ +− models is indeed observed, with the latter exhibiting a monotonically decreasing D xx as a function of the gap amplitude. The most striking scenario with negative D in the +− model is achieved when ∆ exceeds a large modeldependent value. Since the conventional contribution is positive definite, the negative superfluid weight must be attributed to the quantum geometric effects alone. Showing a negative D for the uniform phase is insufficient to judge the plausibility of any PDW order. To this end, we introduce a new superfluid weight defined on the basis of the corresponding enlarged unitcell of the PDW order. Similar to its original definition, such a quantity measures the phase stiffness among the enlarged unitcells. We take the∆ +− configuration and, for illustrative purpose, consider what we call a π-PDW -one that has modulation wavevector (π, 0) and sees the pairing potential change sign every other site in the x-direction (illustrated in the inset of Fig. 2). The calculation will adopt a new velocity operator expressed in the enlarged basis with two lattice sites per unitcell. Details are provided in the Supplementary.
The superfluid weight D xx and D yy defined for π-PDW as a function of ∆ are also plotted in Fig. 1 (b) for comparison. The two components differ because the π-PDW breaks the four-fold rotational symmetry. Note that the xx-component exhibits an apparent cusp at around ∆ = 2t. In close inspection, this coincides with a qualitative change in the quasiparticle spectrum, i.e. going from having a Bogoliubov Fermi surface [26] to being fully-gapped. Intriguingly enough, both of the two components remain positive even when the uniform superfluid weight turns negative. This suggests, at the very least, that a π-PDW phase has the potential to materialize as real physical order as the uniform phase becomes unstable towards phase decoherence.
Stabilization of PDW phase.-The above analyses is based on the assumption that the PDW pairing potential develops the same magnitude as the uniform phase.
It is so far unclear whether the PDW could become more favorable under a given set of pairing interactions. In the following, we self-consistently determine the relative stability between the uniform and π-PDW phases. Note that the actual ground state could well be a more complex form of PDW with a different modulation wavevector. However, finding the ground state is not our primary objective here. It suffices to demonstrate the preference of at least one certain PDW order over the uniform phase.
We consider onsite pairing interactions H int = i H int,i , where on each site i, Here, U ss , U dd , U sd > 0 designate the strength of the intraorbital (ss and dd) interactions and the interorbital (sd) pair hopping. Among them, U ss and U dd promote onsite intraorbital spin-singlet pairings ∆ s and ∆ d within the two respective orbital manifolds. A sizable interorbital repulsion U sd is chosen to ensure that ∆ s and ∆ d condense into the∆ +− configuration. The order parameters ∆ s and ∆ d are obtained by self-consistently solving the gap equations at zero-temperature. The calculations are performed using the same tight-binding model as in Fig. 1, but with a unitcell stretching two neighboring lattice sites in the x-direction to accommodate the π-PDW (see Fig. 2). Throughout the calculation, the electron filling is kept roughly constant by properly adjusting the chemical potential. Owing to the large interorbital mixing t ′ , the pairings of the two orbitals are strongly coupled, which results in comparable |∆ s | and |∆ d | (not shown in Fig. 2). Thus the superfluid weight of a uniform or PDW phase with a self-consistently obtained pairing amplitude could also be roughly estimated by referring back to Fig. 1 (b).
Keeping the relative strength of the interactions at (U ss , U dd , U sd ) = (1, 0.95, 1.2)U , Fig. 2 plots the phase diagram of the model as a function of U . Generally speaking, the pairing becomes progressively stronger as the interaction increases. In the mean time, the uniform phase features increasingly reduced superfluid weight, in agreement with the preceding analyses. The π-PDW phase sets in beyond a critical interaction strength. Notably, the transition occurs before the uniform phase reaches negative superfluid stiffness. This can be attributed to the much stronger pairing developed in the PDW than in the uniform phase, as is shown in Fig. 2. We note that the PDW phase features positive superfluid weight throughout the concerned range of interaction strength.
On the other hand, no transition to the PDW phase has been found in the∆ ++ configuration (favored for small or negative U sd ). This further testifies how the negative geometric superfluid weight facilitates the formation of PDW order. Following the derivation in the maintext, the diamagnetic and paramagnetic contributions to the superfluid weight can be approximated by replacing the velocity operator and the Green's function in Eq. (3) in the maintext by their band basis counterparts. It is straightforward to obtain the following, ,(S1) Here we make use of the relations V mm µk = −V mm µk and (V mn µk ) * = V nm µk . Note that there are terms from the two contributions that are opposite to each other. Taken together, the total superfluid weight follows as, Here we use the fact in our s-d xy model that V mn µk = −V mn µk for m = n.
Here k x ranges from −π/2 to π/2 due to the doubling of lattice constant in this direction, and the spin indices have been suppressed. Note that in the above construction, we have chosen the gauge under which the relative position between the sublattice sites are manifested. Such a gauge is important to obtaining the correct forms of the velocity operators, as electron hoppings on all of the bonds must be accounted for on equal footing. In this gauge, the (normal state) velocity operator is straightforwardly obtained by taking the derivative of (S7), i.e.V 0,µk = ∂ kµĤ0,k . The BdG Hamiltonian of the PDW phase can be written in the Nambu spinor basis (c Ask↑ , c Bsk↑ , c Adk↑ , c Bdk↑ , c † Ask↓ , c † Bsk↓ , c † Adk↓ , c † Bdk↓ ) T as,Ĥ In our study, the pairing matrix with on-site spin-singlet intra-orbital pairing on the two respective orbitals is given by,∆ With the interactions given in Eq. (12) of the maintext, the matrix elements in (S9) can be self-consistently determined by solving the following gap equations, ∆ P s = k −U ss c P sk↓ c P sk↑ + U sd c P dk↓ c P dk↑ , ∆ P d = k −U dd c P dk↓ c P dk↑ + U sd c P sk↓ c P sk↑ , P=A, B. (S10) Note that the summation over k x ranges from −π/2 to π/2 as mentioned above.