Geometry-independent superfluid weight in multiorbital lattices from the generalized random phase approximation

The superﬂuid weight of a generic lattice model with attractive Hubbard interaction is computed analytically in the isolated band limit within the generalized random phase approximation. Time-reversal symmetry, spin rotational symmetry, and the uniform pairing condition are assumed. It is found that the relation obtained in Huhtinen et al. [Phys. Rev. B 106 , 014518 (2022)] between the superﬂuid weight in the ﬂat band limit and the so-called minimal quantum metric is valid even at the level of the generalized random phase approximation. For an isolated, but not necessarily ﬂat, band it is found that the correction to the superﬂuid weight obtained from the generalized random phase approximation D (1)s = D (1)s , c + D (1)s , g is also the sum of a conventional contribution D (1)s , c and a geometric contribution D (1)s , g , as in the case of the known mean-ﬁeld result D (0)s = D (0)s , c + D (0)s , g , in which the geometric term D (0)s , g is a weighted average of the quantum metric. The conventional contribution is geometry independent, that is, independent of the orbital positions, while it is possible to ﬁnd a preferred, or natural, set of orbital positions such that D (1)s , g = 0. Useful analytic expressions are derived for both the natural orbital positions and the minimal quantum metric, including its extension to bands that are not necessarily ﬂat. Finally, using some simple examples, it is argued that the natural orbital positions may lead to a more reﬁned classiﬁcation of the topological properties of the band structure.

s,g is a weighted average of the quantum metric.The conventional contribution is geometry independent, that is independent of the orbital positions, while it is possible to find a preferred, or natural, set of orbital positions such that D (1) s,g = 0. Useful analytic expressions are derived for both the natural orbital positions and the minimal quantum metric, including its extension to bands that are not necessarily flat.Finally, using some simple examples, it is argued that the natural orbital positions may lead to a more refined classification of the topological properties of the band structure.

I. INTRODUCTION
An important result of Bardeen-Cooper-Schrieffer theory is the prediction for the superconductive critical temperature T c ∝ e − 1 |U |ρ 0 (E F ) [1], where U is the interaction strength of the effective attractive interaction and ρ 0 (E F ) is the electronic density of state at the Fermi energy.However, the mean-field critical temperature T c gives the energy scale for the breaking of Cooper pairs, but in many unconventional superconductors and in two-dimensional fermionic superfluids the transition to the normal state is controlled by the phase fluctuations of the order parameter rather than Cooper pair breaking [2][3][4][5][6].The quantity that measures the phase stiffness of the order parameter phase with respect to perturbations or thermal fluctuations is known as the superfluid weight D s .The superfluid weight can be obtained experimentally from the penetration depth of the magnetic field characterizing the Meissner effect in supercondutors.From a theoretical point of view, the superfluid weight is defined as a specific limit of the current-current response function [7,8].Equivalently, it can be computed from the change of the free energy due to a twist in the boundary conditions [9].In this work, we adopt this second characterization since it is the most practical for our purposes.
In the rather idealized but popular model of a superconductor as a system of charged particles propagating in a homogeneous medium, the superfluid weight is simply given by D s = q 2 n m , where q and m are the electric charge and the mass of the particles, respectively, and n is the number density [10].Due to the crystalline structure, real materials possess discrete rather than continuous translational symmetry and, as consequence, the inverse mass 1/m in the superfluid weight has to be replaced by the inverse effective mass tensor [ 1  m eff ] i,j = 1 ℏ 2 ∂ ki ∂ kj ε nk obtained from the energy dispersion ε nk of the partially filled band (see (55) in the following).For a long time it was assumed that only the band dispersion plays a role in determining the superfluid weight, which would then be necessarily small if the charge carriers have a large effective mass [11].In particular the superfluid weight should be strictly vanishing in the so-called flat band limit, in which the band dispersion is constant ε nk = ε n as a function of quasimomentum k.
The conventional wisdom was challenged by few specific examples of systems with flat bands for which the superfluid weight can be shown to be nonzero, namely exciton condensates in quantum Hall bilayers [12,13] and the surface states in rhombohedral graphite [14,15].Later on, it was shown at a general level that even in the flat band limit the superfluid weight can be nonzero and large due to an additional contribution associated to the band wave functions rather than the band dispersion [16][17][18].In particular, under some symmetry assumptions, it is shown that the superfluid weight in the flat band limit is proportional to the integral over the first Brillouin zone of the quantum metric, a geometric invariant constructed out of the Bloch wave functions and their derivatives with respect to the quasimomentum [16,19,20].More generally, under the same symmetry assumptions but not necessarily in the flat band limit, the superfluid weight is a sum of two contributions D s = D s,c + D s,g .The conventional contribution D s,c depends on the inverse effective mass, while the so-called geometric contribution D s,g is a weight averaged of the quantum metric over the Brillouin zone (see ( 50)- (51) in the following).
It is important to note that the geometric contribution to the superfluid weight can be nonzero only in multiband/multiorbital lattices, since periodic Bloch functions live in orbital space, whose dimension is the number of orbitals per unit cell.For a simple lattice with one orbital per unit cell, such as the square lattice, a periodic Bloch function becomes a simple scalar, which cannot affect any physical property, and the quantum metric vanishes.The geometric contribution is thus an effect genuinely associated to the multiorbital character of lattice fermions.The geometric contribution is sizable for instance in multiorbital lattices with flat bands that can be realized with ultracold gases in optical lattices, such as the Lieb lattice [17], and also in iron-based superconductors such as FeSe [21], in which multiorbital effects are important [22,23].Another example of material with strong multiorbital character is twisted bilayer graphene [24,25], in which quasi-flat bands accompanied by the onset of a superconducting state are obtained by tuning the twist angle between the two graphene layers to a specific value, called the magic angle.Multiple theoretical works have reached the conclusion that the geometric contribution to the superfluid weight is comparable to the conventional one in magic angle-twisted bilayer graphene [26][27][28][29] and, more recently, experimental evidence for this prediction has been provided [30].The interplay between quantum geometry and superconductivity in twisted multilayer systems is the topic of two recent review articles [31,32].The quantum metric has also been shown to affect the superfluid properties of lattice bosons [33][34][35].
The quantum metric is intimately related to the Berry curvature, another band structure invariant that plays a crucial role in the quantum Hall effect.Indeed, the two quantities are respectively the real and imaginary parts of the quantum geometry tensor, a positive semidefinite complex matrix obtained from the Bloch functions and their derivatives [16].Due to positive definiteness, there are relations between them expressed by inequalities.While the Berry curvature and the associated Berry phase has been extensively studied for instance in the context of the quantum Hall effect [36], the semiclassical theory of electronic motion [37] and the modern theory of polarization [38], the role of the quantum metric in determining various observable properties is currently the subject of ongoing research [39][40][41][42][43][44][45][46].It is clear by now that the quantum metric is relevant for many phenomena other than superfluidity in multiorbital lattices.
The relation between superfluid weight and quantum metric within the mean-field approximation is by now a rather established fact, nevertheless in Refs.[47,48] it was pointed out that these quantities have an unphysical dependence on the positions of the lattice sites, more specifically on the basis vectors determining the relative positions of different sublattices within the unit cell [49].Intuitively, one expects that the superfluid weight should not depend on the relative positions of the lattice sites, called in the following orbital positions.The idea that certain physical observables are independent of the geometry of the lattice has been discussed in Ref. [50].In the mathematical physics literature the related concept of unit cell consistency has also been introduced [51,52].The quantum metric is neither geometry independent nor unit cell consistent and this leads to the unphysical result that the superfluid weight can be nonzero even in the case of a trivial flat band realized in a lattice model composed of completely disconnected unit cells.The simplest example of this unphysical phenomenon is probably the Su-Schrieffer-Heeger model discussed in Sec.VI A of the present work.
The reason behind this inconsistency has been investigated already in Ref. [47] (see also the lecture notes [53] for a more pedagogical presentation) and is due to the fact that the dependence of self-consistently calculated quantities on the change of boundary conditions is neglected when applying mean-field theory.This is a widely used approximation, equivalent to the prescription of replacing the many-body Hamiltonian with the mean-field Hamiltonian when calculating response functions [8,9], and has the well-known drawback of breaking gauge invariance [54].More specifically, the superfluid weight is computed as the second derivative of the free energy (or thermodynamic grand potential) with respect to the phase angle parametrizing twisted boundary conditions.To explicitly preserve translational invariance, it is convenient to implement twisted boundary conditions by means of a constant electromagnetic vector potential A, which cannot be gauged away in a finite system with torus geometry [9].The usual approximation is to replace the exact free energy with the mean-field free energy, which depends on quantities that need to be computed self-consistently, such as the Hartree-Fock potential and, in the case of superconducting systems, the pairing potential.The mean-field potentials depend on the boundary conditions, which means that they need to be computed self-consistently for each different value of A. It is shown in Ref. [47] that the geometry independence of the superfluid weight is restored by taking into account the A-dependence of the pairing potential alone.Indeed, the crucial new result of Ref. [47] is that the superfluid weight in the flat band limit is proportional to the first Brillouin zone integral of the quantum metric minimized with respect to all possible orbital positions.This so-called minimal quantum metric is a geometry independent quantity and it singles out a specific choice of the orbital positions providing a physically sensible result for the superfluid weight.
The purpose of the present work is to extend the results of Ref. [47] in many ways.In Sec.II, we prove rigorously that the superfluid weight is a geometry independent quantity as a consequence of gauge invariance.More specifically, it is shown that a shift in the orbitals positions amounts to a gauge transformation, which does not affect the value of the thermodynamic grand potential.Following Ref. [9], in Sec.III we introduce mean-field theory as a variational approximation for the grand potential.The mean-field grand potential depends on A either directly, through the Peierls phase, or indirectly through the variational parameters that enter in the mean-field Hamiltonian, namely the Hartree-Fock potential and the pairing potential, which are calculated self-consistently for each value of A. A gauge-symmetry preserving approximation for the superfluid weight is obtained by retaining the A-dependence of the pairing potential, as suggested in Ref. [47].An important difference with this latter work is that here we also retain the Adependence of the Hartree-Fock potential.As shown in Ref. [9], this is equivalent to computing the superfluid weight within the generalized random phase approximation (GRPA).The variational approach used here and in Ref. [9] has the advantage to make it clear that the GRPA is simply mean-field theory applied in a fully conserving fashion by taking into account the fluctuations of the mean-field potentials, and is not, strictly speaking, a beyond mean-field approximation.In Sec.III and in Appendix A, we derive the expression for the superfluid weight within the GRPA in the terms of various correlation functions between quadratic operators evaluated at A = 0.The derivation is carried out using an alternative method compared to Ref. [9] and more similar to Ref. [47], first in the case of a general interaction term and later specialized to the Hubbard interaction term.In Sec.IV, we compute analytically the superfluid weight within the GRPA in the case of a generic lattice model with an Hubbard interaction term.We adopt the same symmetry assumptions that allow to derive the relation with the quantum metric in an isolated, but not necessary flat, band [18].To make this work selfcontained and easier to read, we first repeat with our notation the derivation of the conventional and geometric contributions to the superfluid weight within the simplest non-gauge invariant mean-field theory.This derivation is found originally in Ref. [18].Then in Sec.IV B, we extend this derivation to provide a closed-form expression for the correction term to the superfluid weight that is obtained by taking into account the dependence of both the Hartree-Fock and pairing potentials on the constant vector potential A. An important result of Sec.IV B is that this GRPA correction term is also the sum of a conventional part and a geometric part, which depend on the derivatives of the band dispersion and the periodic Bloch functions, respectively.The conventional GRPA correction is geometry independent and always leads to an increase of the superfluid weight compared to the mean-field result.The expression for the conventional part of the GRPA correction is an original result of our work.The conventional GRPA term appears because we take into account the Hartree potential in our case, contrary to Ref. [47] (the Fock potential vanishes for the Hubbard interaction in the absence of magnetic order).
On the other hand, the geometric part of the GRPA correction is essential to cure the problem found in Ref. [47] of the unphysical dependence of the superfluid weight on the orbital positions.Compared to Ref. [47], we provide an explicit expression for the geometric GRPA correction even for an isolated band that is not necessarily flat.
In Sec.V the geometric part of the GRPA correction is analyzed in more detail.More precisely, we show that the sum of the geometric GRPA correction and meanfield geometric contribution gives the minimal quantum metric in the flat band case, as found in Ref. [47].Simple analytical expressions are also provided for both the minimal quantum metric and the associated natural orbital positions that minimize the integrated quantum metric, see Eqs. ( 84)-( 88) and (90).At the end of Sec.V, we explain how these results can be straightforwardly extended to an isolated band that is not necessarily flat.Since we expect the minimal quantum metric and the natural orbital positions to find applications also beyond the context of superfluidity in multiorbital lattices, in Sec.VI we provide some examples, that is we compute these band structure invariants for three different representative lattice models: the Su-Schrieffer-Heeger model, the Creutz ladder and the dice lattice.Finally, in Sec.VII we summarize and discuss our results and single out possible directions for further work.
Appendices A, B and C contain all the necessary computational details that should allow the interested reader to understand, reproduce and ultimately extend and apply our results to other related problems.Appendix A provides a derivation of the GRPA result for the superfluid weight, which in the case of the Hubbard interaction is given by Eqs. ( 30)- (35) in Sec.III.This derivation is different but equivalent to the one of Ref. [9].Appendix B collects useful results for computing the correlation functions that enter in the GRPA expression for the superfluid weight.Finally, in Appendix C the self-consistency equations of mean-field theory are derived in the case of the Hubbard interaction.

II. GEOMETRY INDEPENDENCE OF SUPERFLUID WEIGHT AND GAUGE INVARIANCE
In this section we introduce the superfluid weight as the second derivative of the thermodynamic grand potential Ω(A) with respect to a constant electromagnetic vector potential A. The vector potential parametrizes twisted boundary conditions in a finite size lattice model with a torus geometry.Moreover, we provide a simple argument based on gauge invariance showing that the superfluid weight is a geometry independent quantity.This section contains only a short summary of the concepts that are needed for the present work.For a more extensive and rigorous presentation, the reader should consult Refs.[9] and [47], which are the basis for the results presented here.The notation used here is essentially the same as the one of Ref. [9].
We consider a generic lattice model described by a noninteracting, or free, Hamiltonian that is quadratic in the fermionic field operators ĉiασ , ĉ † iασ , and depends parametrically on a constant electromagnetic vector potential A through the usual Peierls phase Here H σ free is the hopping matrix, which is translationally invariant since its matrix elements [H σ free ] iα,jβ = [H σ free (i−j)] α,β depend only on the difference between the unit cell indices i = (i 1 , i 2 ) T and j = (j 1 , j 2 ) T , and r iα,jβ is the displacement vector from site jβ to site iα, where α, β = 1, 2, . . ., N orb label the orbitals inside the unit cell.The number of orbitals is denoted by N orb , while N c is the number of unit cells.Without loss of generality, only two-dimensional lattice models are considered in the present work.Note also that the free Hamiltonian commutes by construction with the z-axis spin component operator Periodic boundary conditions are imposed by starting with an infinitely extended and translational invariant lattice, that is a collection of lattice sites located at positions r iα .As a consequence of translational invariance, the site positions transform as follows under shifts of the unit cell index The fundamental vectors a 1 and a The procedure just presented for obtaining a finite size lattice model with periodic boundary conditions from an infinite extended one has been introduced and explained in more detail in Ref. [55].
It is crucial to note that in a finite size lattice with periodic boundary condition, it is not possible to write the displacement vectors as the difference of the site positions, namely r iα,jβ ̸ = r iα − r jβ , otherwise it would be possible to eliminate the Peierls phase e iA•r iα,jβ in (1) by means of a gauge transformation and the constant vector potential A would have no observable effects.Instead, one has to define the displacement vectors more carefully, as explained in Ref. [9].For periodic boundary conditions, the vector potential A affects the eigenvalues and eigenfunctions of the free Hamiltonian (1) since it correspond to magnetic fluxes through the holes of the torus.It can be shown using a gauge transformation, that a nonzero A is equivalent to twisted boundary conditions [55].
The free Hamiltonian (1) is translationally invariant, thus it is convenient to expand the field operators in their Fourier components where the wave vectors k are discretized according to the relation k • R j = 2πn j for integers n j .Inserting the Fourier expansion (4) in (1) leads to In the following we need the average current density operator (current operator for short) defined as Note that here we have not normalized the current density operator by the area of the system as done in Ref. [9].In order to formulate the concept of geometry independence, introduced in Ref. [50], we define an operator encoding a shift in the orbital positions b = The is the shift of the position of the orbital labeled by α since the new orbital positions r ′ iα and displacement vectors r ′ iα,jβ are given The hopping Hamiltonian H σ free (k) in momentum space depends on the choice of the position vectors and the displacement vectors, in fact one has where in the last equation b is a single-particle operator in orbital space, whose eigenvalues are the orbital position shifts b α , namely b |α⟩ = b α |α⟩.The operator b generates gauge transformations parameterized by A Using the shifted displacement vectors, one obtains a new noninteracting Hamiltonian The current operator also changes accordingly In this work, we consider the case of an attractive Hubbard interaction term Thus, the full many-body Hamiltonian is with N = σ iα niασ the particle number operator.Note that the interaction term is invariant under the gauge transformations (10), namely [ Ĥint , Û (A)] = 0, as a consequence only the noninteracting part of the manybody Hamiltonian is modified by the gauge transformation ( 10) Due to the property of similarity-invariance of the trace Tr[H] = Tr U HU −1 , it follows that the thermodynamic grand potential is invariant under gauge transformations (β = 1/(k B T ) is the inverse temperature).The superfluid weight is defined as the second derivative of the grand potential with respect to the uniform vector potential A where is the area of the finite size system with periodic boundary conditions.From (16), it is evident that the superfluid weight is a geometry independent quantity [50], i.e. it is invariant under shifts of the orbital positions inside the unit cell.Notice how this result is a consequence of the fact that, under an orbital position shift, the displacement vectors in (8) are modified by a difference term b α −b β , while the displacement vector r iα,jβ cannot be written as a difference, as discussed above.
The concept of geometry independence, or dependence, has been discussed in Ref. [50] and its utility has been pointed out in relation to several observable quantities.However, the superfluid weight is the first observable that has been shown to be geometry independent as a consequence of gauge invariance.The same approach may be useful to prove the geometry independence of other quantities in the future.

III. GENERALIZED RANDOM PHASE APPROXIMATION
An attractive interaction, such as the one in (13), generally leads to the emergence of superfluid phases characterized by a nonzero superfluid weight [7,8].In the context of superfluidity and superconductivity, mean-field theory is the simplest and most common approximation.This is a variational approximation for the grand potential based on the Bogoliubov inequality [9,56,57] The auxiliary grand potential Ω 0 = −β −1 ln Tr e −β Ĥ0 is obtained from a quadratic variational Hamiltonian Ĥ0 , which, in the case of the Hubbard interaction (13), takes the form The expectation value on the right hand side of ( 18) is evaluated with respect to the statistical ensemble associated to the variational Hamiltonian Ĥ0 , namely The coefficients Γ σ α and ∆ α that appear in Ĥ0 are variational coefficients that are chosen so as to minimize the mean-field grand potential Ω m.f. on the right hand side of (18).When the minimum of Ω m.f. is attained, the Hartree potential Γ σ α and the pairing potential ∆ α satisfy the self-consistency conditions of mean-field theory where the expectation values are evaluated as in (22).
In (23) we define σ =↓ if σ =↑ and viceversa.In the equations above, we have assumed that translational symmetry is not broken, therefore the expectation values ⟨ĉ † iασ ĉiασ ⟩ and ⟨ĉ iα↓ ĉiα↑ ⟩ do not depend on the unit cell index i.It is also assumed that spin rotational symmetry around the z-axis is preserved by the mean-field solution, indeed the variational Hamiltonian (19) commutes with the spin operator (2).This implies that the Fock meanfield potential term, proportional to ĉ † iα↑ ĉiα↓ in the case of the Hubbard interaction, does not appear in Ĥ0 .This assumption is justified since spin rotational symmetry breaking is associated to magnetic order, which generally occurs for repulsive interactions.
The mean-field grand potential is minimized with respect to ∆ α and Γ σ α for each separate value of the vector potential A, therefore the Hartree and pairing potentials become themselves function of A. It has been shown in Ref. [47] that one cannot ignore this dependence in the pairing potential ∆(A) when the mean-field approximation for the superfluid weight is evaluated as We denote by dΩ m.f.A, ∆(A), Γ(A))/dA l the full derivative of the mean-field grand potential, including also the A-dependence of the mean-field potentials Γ σ α (A) and ∆ α (A), while the partial derivative ∂Ω m.f.A, Γ(A), ∆(A))/∂A l denotes the derivative with respect to the first argument only.Replacing the full derivatives with the partial derivatives in ( 25) is a commonly used approximation [8,16,18].However, it has the disadvantage of breaking gauge invariance and thus geometry independence [47].On the other hand, it has been shown [9] that the superfluid weight computed from (25) (with the full derivatives) is in fact equivalent to the generalized random phase approximation [58,59], which is an approximation that preserves gauge invariance.
As shown in Appendix A, it is possible to express the full second derivatives of the mean-field grand potential in (25) in terms of correlation functions evaluated on the mean-field statistical ensemble for A = 0 only.This is advantageous since it is not necessary to solve the meanfield problem for several different values of A in order to evaluate the superfluid weight.To present this result, we introduce the following convenient notation where Â and B are arbitrary operators and the notation Â(τ ) = e τ Ĥ0 Âe −τ Ĥ0 for the time evolution of an operator in imaginary time τ has been used.In fact, the symbol introduced in (26) is simply the imaginary time Green's function evaluated at the Matsubara frequency iω n = 0 [60] and has some useful properties that are easy to prove In Appendix B, we introduce a compact method to evaluate the correlation function in (26) between two translationally invariant quadratic operators Â and B.
To express the result for the superfluid weight in the generalized random phase approximation, we need correlations functions of the form (26) of pairs of operators taken from the set { Ĵl , Nασ , Dα }.The correlation functions that involve the components of the current operator Ĵl (6) are organized into a vector . . .
Instead, all the remaining correlation functions are collected into a matrix Finally, we need also the following matrix As shown in Appendix A, the full derivative of the meanfield grand potential can then be expressed as In the last equality we have used the geometric series expansion As mentioned above, the superfluid weight is often computed by retaining only the second partial derivatives of the mean-field grand potential in (35), given by The first term on the right hand side is known as the diamagnetic part of the current-current response function, while the second term is the paramagnetic one [8,9].It can be shown [18] that the diamagnetic part is equal to a correlation function of the form (26) (see (B43)).

IV. SUPERFLUID WEIGHT IN THE ISOLATED BAND LIMIT
The aim of this section is to compute analytically the superfluid weight in the isolated band limit within the generalized random phase approximation, which means that ( 35) is reduced to integrals over the first Brillouin zone of certain combinations of the band dispersions and band wave functions, and of their derivatives with respect to quasimomentum.Ultimately, these integrals have to be evaluated numerically, however, in Sec.VI, we provide also some examples in which fully analytical results can be obtained.The final expression presented below is valid for generic lattice models under few assumptions, the most important being time-reversal symmetry and the uniform pairing condition to be introduced in the following.For completeness, we first evaluate the terms corresponding to the second partial derivatives of Ω 0 (37) and reobtain the known result that the superfluid weight in a multiband/multiorbital lattice can be separated into two contributions, called conventional and the geometric, respectively [16,18].In particular, the quantum metric enters in the geometric contribution to the superfluid weight.However, as pointed in Refs.[47,48], the quantum metric depends on the orbital positions, therefore it is not a geometry independent quantity.In order to restore gauge invariance and thus geometry independence, we evaluate the remaining terms in (35).This amounts to computing the full derivatives with respect to the vector potential A rather than just the partial derivatives.The same approach has been used in Ref. [47], with the only difference that the Hartree potential Γ σ α in the variational Hamiltonian Ĥ0 (A) is neglected in this latter work.By slightly modifying the derivation in Appendix A, it is shown that neglecting the Hartree potential amounts to setting to zero in (35) all the correlation functions in which the number operators Nασ appear.One of the main result of this work is to take into account the Hartree potential and show that in this way one obtains a new correction to the superfluid weight that is geometry independent and is proportional to the derivatives of the band dispersion.Therefore, in the language of Refs.[17,18], this is a conventional contribution since it vanishes in the flat band limit.
We start by rewriting the variational quadratic Hamiltonian Ĥ0 in Nambu form where the column (row) vector ĉk (ĉ † k ) is defined in (B4) and H 0 (k) is the single-particle Bogoliubov-de Gennes (BdG) Hamiltonian given by The Nambu form for translational invariant quadratic operators is discussed in Appendix B. From now on time-reversal symmetry is assumed, which implies The second assumption, which enables the analytic evaluation of the superfluid weight for generic lattices, is called the uniform pairing condition, expressed by The pairing potential ∆ can also be taken real and positive.With a slight abuse of notation, we indicate with ∆ both the scalar value of the uniform pairing potential in (42) and the matrix in (41), which becomes proportional to the identity.The uniform pairing conditions is justified in Appendix C starting from the self-consistency equations of mean-field theory ( 23)- (24).Under the above assumptions, the BdG Hamiltonian can be diagonalized as follows Here, ) the diagonal matrix of the band dispersions.Finally, the BdG wave functions are given by The quantities u nk and v nk are the usual BCS coherence factors [54,61].Note that in a multiband lattice model there is a pair (u nk , v nk ) of coherence factors for each band.
The total superfluid weight is separated into two contributions The first contribution D (0) s is the one given by the second partial derivative of Ω 0 , that is (37), while D (1) s includes all the remaining terms in (35).In the following, D (0) s is called the "mean-field theory (MFT) superfluid weight", because in all previous works, with the exception of Refs.[47,48], this is the only term that is evaluated when the superfluid weight is computed within the mean-field approximation.Instead, D s is referred to as the "GRPA correction"to the superfluid weight.In fact, it would not be incorrect to consider the sum of the two contributions in (48) the actual mean-field theory result for the superfluid weight since it is obtained by taking the full derivative of the mean-field free energy Ω m.f., see (35).However, we avoid this nomenclature in order not to create confusion when referring to previous works.
As shown in the following, both the MFT superfluid weight and the GRPA correction can be written as the sums of a conventional contribution and a geometric contribution, which we indicate as D s,g , with j = 0, 1.As explained in Appendix B, the conventional contribution depends only on the intraband matrix elements of the current operator.According to (B58), this means that only the derivative of the band dispersions ∂ε nk /∂k l with respect to the quasimomentum k, that is the group velocity, enters into D (j) s,c and not the derivatives of the Bloch functions |∂ l g nk ⟩.On the other hand, the geometric contribution D (j) s,g is associated to the interband matrix elements of the current operator, therefore it depends only on scalar products of the form ⟨g mk |∂ l g nk ⟩, but not on the group velocity.It was pointed out for the first time in Refs.[16,18] that the derivatives of the Bloch functions affect the mean-field superfluid weight in the form of an additional geometric contribution, indicated by D (0) s,g in our notation.A major result of the present work is to show that, under the same assumptions, the GRPA correction D s,g to be presented below.

A. MFT superfluid weight
Using the results in Appendix B, more specifically by taking the isolated band limit in (B67), we obtain for the MFT superfluid weight [ [ where we have introduced the quantum metric defined in terms of the projector P (k) = |g nk ⟩⟨g nk | on the only partially filled band labeled by n.The partially filled band in the noninteracting limit is determined by the condition min k ε nk ≤ µ ≤ max k ε nk on the chemical potential µ.Isolated band limit means that the partially filled band is separated from all other bands by a large band gap E gap Moreover, it is assumed that the interaction strength is much smaller than the band gap In this limit, the superfluid weight and all of the observable properties of the system are determined by the band dispersion ε nk and the Bloch functions |g nk ⟩ of the partially filled band.Rather counterintuitively, in order to obtain the correct expression for the geometric contribution to the MFT superfluid weight it is essential to take into account the interband matrix elements of the current operator, as explained in Appendix B. Naïvely, one may neglect the interband matrix elements of the current operator in the isolated band limit, but the result is that the geometric contribution is lost.In general, the geometric contribution is considerably smaller than the conventional one for a dispersive band with a bandwidth much larger than the interaction strength, so the former can be safely neglected [18,32].On the other hand, when the partially dispersive band is quasi-flat, this not anymore an acceptable approximation, since the conventional contribution vanishes in the flat band limit (∂ l ε nk = 0) while the geometric one does not.It is possible to rewrite the conventional contribution in a more suggestive form by using (B59) in the form given by the last term and performing an integration by parts.In this way, one obtains the following expression with The first term inside the parenthesis is the inverse effective mass tensor weighted by the occupation factor n k .It is not difficult to show that n k is precisely the occupation number (including both spins) of the state with energy ε nk at thermal equilibrium.Therefore, the interpretation of the first term is that all of the particles participate to the superfluid flow at zero temperature and the superfluid weight is an average measure of the effective mass of the carriers of the superfluid current.The second term vanishes at zero temperature and can be intepreted as a depletion of the superfluid component due to the thermal excitations of quasiparticles [62].However, this interpretation does not provide the full picture since the geometric contribution ( 51) is not taken into account.

B. GRPA correction to the superfluid weight
In the previous section we have reobtained the general result for the MFT superfluid weight in multiband/multiorbital lattices in the presence of an Hubbard interaction term, which has been discussed in a number previous works.In this section, we move on to consider the GRPA correction to the superfluid weight, consisting of the terms in (35) which depends on the matrix A and the vector v l .Thus, we need to evaluate their these objects in the isolated band approximation.The details of the calculation are again provided in Appendix B. For the matrix elements of the matrix A α,β defined in (31), we obtain with the two 4 × 4 matrices C c k and C g k given by For reason that will become clear in a moment, we distinguish two contributions to the matrix A, a conventional one proportional to C c k and a geometric one proportional to C g k .The same separation applies to the vector v l (30) Ĵl These results are obtained from (B55)-(B56), (B61) and (B64) in Appendix B. It is clear from (65) that v c lα is a conventional term since it depends only on the group velocity ∂ l ε nk .In fact, expressing v c lα in terms of the group velocity alone is a nontrivial result that involves a fair amount of algebra, see (B60).On the other hand, one can see from (66) that v g lα is purely geometrical, since it depends only on the derivatives of the Bloch functions, or equivalently, only on the derivatives of the flat band projector ∂ l P (k).
The key observation that allows to considerably simplify the calculation is to notice that the vectors v c,g lα are eigenvectors of the matrices C (c,g) k and B α (33) As a consequence, we have also with It is then clear that GRPA correction to the superfluid weight can be written as the sum of a conventional and geometric components since the vectors introduced in (64) are orthogonal v c lα • v g mβ = 0.The conventional component takes the form [D (1)  s Here, A c = A/N c denotes the unit cell area, U = diag(U 1 , U 2 , . . ., U α=N orb ) is a diagonal matrix with the coupling constants U α on the main diagonal, λ c is a matrix with components given by ( 71) and c l = (c l,1 , c l,2 , . . ., c l,α=N orb ) T are vectors whose components are given by c lα = A −1 ( Ĵl , Nα↑ ).Recall that, in the thermodynamic limit, the summation over wave vectors in the expression for ( Ĵl , Nα↑ ) in ( 65) is replaced by the Brillouin zone integral k → A (2π) 2 d 2 k.At zero temperature λ c = 0 and only the first term in (73) survives.Since the components of the vector c l are purely real and U α > 0, the matrix M with components [M ] l,m = c T l U c m is positive semidefinite, which means that the conventional part of the GRPA correction leads to an enhancement of the superfluid weight compared to its mean-field value at zero temperature.
It can be also shown that the second term in ( 73) is always negative semidefinite, which is consistent with the discussion below (A20) (see also Ref. [47]).To prove this, one notes that the matrix with c = diag(c 1 , c 2 , . . ., c N orb ) the diagonal matrix obtained from the components of the generic complex vector c = (c 1 , c 2 , . . ., c N orb ) T .In fact, this is just a special case of the Schur product theorem according to which the Hadamard product (entry-wise product) C = A • B of two positive (semi)-definite matrices A and B is positive (semi)-definite [63].In our specific case, we have Q(k) = P (k) • P * (k) with both P (k) and P * (k) positive semidefinite.Thus λ c is negative semidefinite and as a consequence also ((λ c ) −1 − U ) −1 .This concludes the proof.By the same token, is a positive semidefinite matrix, showing that the conventional GRPA correction D s,c always leads to an increase of the mean-field superfluid weight even at finite temperature.It will become clear in the following that the conventional GRPA correction D (1) s,c is also geometry independent.
Using (75) with λ g in place of λ c , we can write the geometric component of the GRPA correction to the superfluid weight in the following form [D (1)  s Similarly to (73), the matrix λ g has matrix elements λ g α,α ′ given by ( 72) and d l = (d l,1 , d l,2 , . . ., d l,α=N orb ) T is a vector with components d lα = A −1 Ĵl , Dα , whose expression in terms of the periodic Bloch functions and their derivatives is given in (66).Note that the components of the d l vector just introduced are purely imaginary since Ĵl , , see ( 27)- (28).Therefore, if it is shown that the matrix U −1 − λ g is positive definite and so is its inverse, then D (1) s,g is negative semidefinite, meaning that superfluid weight is decreased by the geometric GRPA correction compared to the mean-field result D (0) s .This is in contrast to the conventional GRPA correction D (1) s,c discussed previously.A subtle point here is that U −1 − λ g is not invertible, as shown below.This is not a problem since the formula in (76) makes sense and gives the correct result if one denotes by 1/(U −1 − λ g ) the Moore-Penrose inverse (pseudoinverse).
To prove that U −1 − λ g is positive semidefinite, we need the self-consistency equation of mean-field theory for the parameters ∆ α .By assuming the uniform pairing condition (42) and taking the isolated band limit, one obtains from (C7) Therefore, we can write with The matrix R(k) is defined in terms of the Hadamard product Q(k) = P (k) • P * (k) introduced in (74), while diag P (k) is the diagonal part of P (k), that is the matrix obtained by setting to zero all of the matrix elements away from the main diagonal.Since the linear combination of positive semidefinite matrices with positive coefficients is again positive semidefinite, we just need to prove that R(k) is positive semidefinite.This is done as follows where the vector c and the diagonal matrix c are as in (74).This concludes the proof that D s,g is negative semidefinite.By taking c as the identity matrix, one obtains that the right-hand side of ( 80) is zero, showing that R(k) and U −1 − λ g are never invertible.This can be traced back to the fact that the mean-field free energy does not change if all the variational parameters ∆ α are multiplied by the same constant phase e iϕ .The invariance under global phase rotations of the order parameter is a general property of superconducting systems, which is a consequence of gauge invariance.
The expressions ( 73) and ( 76) for the conventional and geometric components of the GRPA correction are the main results of this section, together with the statements regarding their positive or negative semidefiniteness, respectively D s,g ≤ 0. The formulas (73) and (76) allow to compute the full GRPA correction D (1) s of generic lattice models in the isolated flat band limit in terms of the energy dispersion ε nk and the Bloch functions |g nk ⟩ of the only partially filled band.An important difference between the conventional and the geometric components is that the latter is not geometry independent, meaning that it depends on the orbital position vectors r iα that enter in the Fourier transforms of the field operators (4) and of the free Hamiltonian (5).In the next section, it is explained how one can take advantage of this fact and set the geometric GRPA correction to zero by a suitable choice of the orbital positions.

V. MINIMAL QUANTUM METRIC AND NATURAL ORBITAL POSITIONS
In this section, we first assume that the partially filled band is not only isolated, but also flat in order to simplify the presentation.All of the results can be straightforwardly extended to the case of an isolated, but not necessarily flat, band as explained towards the end.In the isolated flat band limit, only the geometric contribution to the superfluid weight s,g survives, while the conventional one vanishes s,c = 0. Since the energy dispersion of the band is just a constant ε nk = ε n, the superfluid weight depends only on an invariant built out of the flat band Bloch functions.The quantum metric ( 52), which appears in the expression (51) for the geometric contribution to the MF superfluid weight, is invariant under multiplication of the Bloch functions by an arbitrary k-dependent phase factor |g nk ⟩ → e iϕ(k) |g nk ⟩, because the projector P (k) = |g nk ⟩⟨g nk | is unaffected by this transformation.In this sense the quantum metric is a band structure invariant.However, if we perform a shift of the orbital positions as in (8), the flat band projector transforms as (see (9)) The components b l of the vector b are operators acting in orbital space and encode the position shifts b α for each orbital, that is b l |α⟩ = [b α ] l |α⟩.The presence of the commutator term in (82) implies that the quantum metric is not geometry independent.It follows that it is essential to include the GRPA correction in order to restore the geometry independence of the total superfluid weight, which has been established in Sec.II using the principle of gauge invariance.On the other hand, it is clear that the quantities ⟨α|P (k)|α⟩ and | ⟨α|P (k)|α ′ ⟩| 2 , appearing in (65) and (71), respectively, are invariant under orbital position shifts (81), therefore the conventional component of the GRPA correction ( 73) is geometry independent.
The purpose of this section is to show that it is possible to set to zero the geometric component of the GRPA correction, that is D s,g = 0, by a suitable choice of the orbital positions.It is also shown that these orbital positions, called in the following the natural orbital positions, are the ones that minimize the trace of the quantum metric integrated over the Brillouin zone, in agreement with the results of Ref. [47].Using the nomenclature of this last reference, the integrated quantum metric computed using the natural orbital positions is called the minimal quantum metric.
Using ( 51) and ( 76), we can write the superfluid weight in the isolated flat band limit in the following way where which is also flat, and the minimal quantum metric M l,m is defined as We call M l,m the minimal quantum metric as in Ref. [47], but a more proper name would be minimal integrated quantum metric, since M l,m is the quantum metric (52) integrated over the whole Brillouin zone.The vector s l is related to the vector d l introduced in (76) and has purely real components.As before, in (88) we have used the Hadamard product As in (76), R −1 denotes the pseudoinverse since R is not an invertible matrix.
It is possible to show that the result in ( 83)-( 88) for the superfluid weight in the isolated flat band limit is applicable also in the case of several degenerate flat bands.The only modification is that the projector reads P (k) = n∈F |g nk ⟩⟨g nk | in this case, where the sum runs over the set F of degenerate partially filled flat bands.An example of a lattice with degenerate flat bands is the dice lattice presented in Sec.VI C.
We now show that, if the orbital positions are chosen so as to minimize the trace of the integrated quantum metric Tr M = l M l,l , then the vectors s l vanish and so does the geometric GRPA correction D (1) s,g proportional to the quantity s T l R −1 s m /2 in (84).To this end, we need to compute how the integrated quantum metric changes under a shift of the orbital positions.From (81), we obtain By setting to zero the derivatives of the trace Tr M ′ = l M ′ l,l with respect to the shift components [b α ] l , it is found that the integrated quantum metric is minimized when the following condition is satisfied The linear system (90) necessarily has a solution because Tr M ′ is a positive quantity and has a minimum as a function of the shifts [b α ] l .However, the solution is not unique since Tr M ′ is invariant if all orbitals are shifted by the same amount ([b α ] l = [b β ] l for all α, β) as one expect from translational invariance.Indeed, if e = (1, 1, ..., 1) T is the vector with all the components equal to one, it is easy to verify that s • e = 0 and Re = 0. Thus, a solution of (90) can be expressed in terms of the pseudoinverse R −1 of the matrix R, namely b l = R −1 s l /2.According to the properties of the pseudoinverse, this solution has minimum norm On the other hand, we have from ( 82) and ( 87) This result together with (89) implies that the minimal quantum metric M (84) is geometry independent, as expected from the discussion in Sec.II.Another consequence is that, if the orbital shifts are chosen so as to satisfy (90), then s ′ l = 0 and the minimal quantum metric M (84) coincides with the integrated quantum metric M, thus it is a positive semidefinite matrix.This also means that the correlation functions (J l , D α ) (66) and the GRPA correction in ( 83)-( 84) vanish if calculated with the orbitals positions defined by (90).Importantly, it is always possible to find such a preferred set of orbital positions, which we call natural orbital positions.In general, we have observed that the natural orbital positions are always unique up to arbitrary translations.This is indeed the case of the examples presented in Sec.V, however a general proof of this fact is lacking at present.
Even though only the flat band case has been considered so far in this section, it is easy to show that natural orbital positions, for which D (1) s,g = 0, exist also in the case of an isolated band that is not necessarily flat.In order to do this, one simply needs to repeat the arguments presented in this section with the inclusion of the weight factor E −1 nk tanh(βE nk /2) under the integral sign d 2 k in ( 85), ( 87) and (88).Indeed, it is immediate to see that the derivations of the transformation properties in (89) and (91) are unaffected by the weight factor.
It is a remarkable fact that the expression ( 84)-( 88) for the minimal quantum metric and the linear system (90) defining the natural orbital positions are obtained also from the analysis of the two-body excitation spectrum.Here, we quickly recall how this has been done in Ref. [47].The dispersion of propagating two-body bound states is given by the effective Hamiltonian first introduced in Ref. [64], under the assumption that for all α and β.This condition guarantees that the uniform pairing condition is satisfied for U α = U β = U (see (C8)) and implies that the two-body bound state with lowest energy for q = 0 is represented by the effective state vector |Ψ 0 ⟩ with ⟨α|Ψ 0 ⟩ = ⟨β|Ψ 0 ⟩ = N −1/2 orb for all α, β.The effective Hamiltonian transforms as h(q) → h ′ (q) = e −iq•b h(q)e iq•b under an orbital position shift, therefore its eigenvalues are geometry independent quantities.Using second order perturbation theory, one finds that the inverse effective mass of the bound state represented by |Ψ 0 ⟩ is proportional to the minimal quantum metric (84), which coincides with the integrated quantum metric M l,m = N orb ⟨Ψ 0 |∂ l ∂ m h(q)|Ψ 0 ⟩, if the natural orbital positions are used in (92).In particular, the linear system determining the natural orbital positions obtained in Ref. [47] (see Eq. (38) in this reference) is a special case of (90) when (93) holds.Indeed, the components of s l can be written in the following equivalent form which can be identified with the right-hand side of Eq. ( 38) in Ref. [47], while on the left-hand side it is also easy to identify the matrix R (88).The second equality in the above equation is a consequence of the fact that the projector P (k) is in general not periodic, but rather ϱ-equivariant [52], namely P (k + g j ) = e −igj •b P (k)e igj •b for some given shift operator b [see ( 9)], while g j are reciprocal lattice vectors defined by a i • g j = 2πδ i,j , therefore ⟨α|P 2 (k)|α⟩ = ⟨α|P (k)|α⟩ is periodic.Our derivation of ( 84)-( 88) and ( 90) is equivalent to the one of Ref. [47] but is based on the direct evaluation of the GRPA expression for the superfluid weight.It is rather instructive to see how the same results can be obtained with completely different methods.This provides strong evidence that the relation between minimal quantum metric and superfluid weight in a flat band is an accurate, or even exact, result.Our approach has the advantage that it is straightforward to relax the flat band assumption, moreover, we have shown that the condition (93) is unnecessary.This means that the minimal quantum metric and the associated natural orbital positions can be defined for arbitrary bands or composite bands, therefore they are likely to find applications in other contexts.With this perspective in mind, we illustrate these band structure invariants by considering some representative examples in the following section.

VI. EXAMPLES
In this section, we illustrate the ideas introduced in the previous section by considering three examples: Su-Schrieffer-Heeger (SSH) model, Creutz ladder and dice lattice.In all of them, the Brillouin zone integrals that enter in the definitions of the minimal quantum metric and natural orbital positions can be worked out analytically.This gives us the opportunity to understand these concepts in the simplest possible setting.See also Ref. [47] for an analogous discussion of the Lieb lattice.

A. Su-Schrieffer-Heeger model
The SSH model is a one-dimensional lattice model widely used to illustrate the topological properties of the band structure and the concept of bulk-edge correspondence [65].As illustrated in Fig. 1, the unit cell of the SSH model contains two orbitals, labeled by α = 1, 2. The single-particle free Hamiltonian H free (1) for the SSH model is given by The two real and positive parameters v and w denote the intra-cell and inter-cell hopping matrix elements, respectively.The lattice constant is fixed to a = 1, thus the crystal momentum k takes values in the range [−π, π].
The orbital positions are chosen to be the same for the 1 2 v w w=0 v=0 FIG.1: The SSH model is a simple linear chain with alternating hoppings.The unit cell, shown in the top row as a green rectangle, consists of two orbitals labeled by α = 1, 2, while v and w denote the hopping amplitudes within a unit cell (thick line) and between neighboring unit cells (thin line), respectively.By convention, the two orbitals in the unit cell have the same position, see (98).For illustration purposes, they are displaced from each other in the figure.The two energy bands of the SSH model are both flat when either w = 0 (middle row) or v = 0 (bottom row).In both cases, the lattice model reduces to a collection of disconnected dimers, therefore transport is not possible in any form.However, the quantum metric and thus the MF superfluid weight is zero in the first case (w = 0) and nonzero if w ̸ = 0. To resolve this inconsistency, it is necessary to take into account the geometric GRPA correction to the superfluid weight, which amounts to computing the quantum metric using the natural orbital positions.The red rectangles in the middle and bottom rows contain orbitals whose natural positions are identical.
two orbitals inside the same unit cell (see top row in Fig. 1), that is r j,(α=1) = r j,(α=2) = j . (98) With this convention for the orbital positions, the Fourier transform of the free Hamiltonian (5) reads It is straightforward to obtain the energy dispersions ε n (k) and projection operators P n (k) for the upper (n = +) and lower (n = −) band The spin index σ has been dropped since since the projector operator P σ ± (k) does not depend on the spin.We need also the derivative of the projector given by Since the diagonal matrix elements of the projector (101) are constant and equal to ⟨α|P ± (k)|α⟩ = 1/2, the SSH model satisfies the uniform pairing condition if U 1 = U 2 = U .Under this condition, the geometric contribution to the MF superfluid weight is given by the formula in (51) in the isolated band limit (adapted to onedimension), which depends on the quantum metric G(k).
We obtain from (102) that the quantum metric is The quantum metric is the same for the two bands and diverges at k = π when |w − v| → 0 concomitantly with the closure of the energy gap between the two bands at the same k point.The integral of the quantum metric over the Brillouin zone M = π −π dk 2π G(k) can be performed by using the change of variable e ik → z and applying the residue theorem to the resulting contour integral for v > w , for v < w .

(104)
We have computed the integrated quantum metric away from the flat band limit, which in the SSH model is obtained when v = 0 or w = 0, since it can be of interest also for other applications besides the computation of the superfluid weight.In the case w = 0, the unit cells are completely disconnected, as shown in Fig. 1, and one expects the absence of all forms of transport.In fact, the superfluid weight vanishes since M = 0. On the other hand, if v = 0, then the quantum metric is nonzero (M = 1/2) even though the lattice is again composed of disconnected two-orbital dimers and one would again expect vanishing superfluid weight.As explained below, this apparent paradox is due to the fact that we have ignored the geometric GRPA correction D (1) s,g or, in other words, we have calculated the quantum metric using orbitals positions that are different from the natural ones.
To compute the natural orbital positions and the GRPA correction, we need the vector s = (s α=1 , s α=2 ) T (86)- (87), in which the spatial index l has been dropped since we are dealing with a one dimensional model, and the pseudoinverse of the matrix R (88).It is straightforward to obtain the latter from (101) Instead, the components of the vector s are calculated again using a contour integral as in ( 104) The position shifts, given by b = (b α=1 , b α=2 ) T = R −1 s/2 vanish in the case v > w, which confirms that the choice in (98) corresponds to the natural orbital positions.On the other hand, for v < w the natural orbital positions are given by As a consequence, we have r j,2 = r j+1,1 .This means that in the SSH model the natural orbital positions are given by assigning the same positions to the orbitals connected by the hopping with largest magnitude, as illustrated in Fig. 1.We can also compute the minimal quantum metric for arbitrary value of the hopping amplitudes.From ( 105) and ( 106), one obtains s T R −1 s = 1 in the case v < w, therefore the minimal quantum metric is given by for v < w . (108) In particular M = 0 for v = 0 or w = 0, thus the superfluid weight vanishes in the flat band limit in the case of the SSH model after taking into account the geometric GRPA correction, as expected.Note that the expressions of the minimal quantum metric in the two cases are related by to each other through the interchange v ↔ w.This is consistent with the fact that one can interchange the two hopping hopping amplitudes v and w by a simple redefinition of the unit cell, which amounts to a different specification of the orbitals positions.Thus, the minimal quantum metric captures an intrinsic property of the band wave functions, which is not affected by the unit cell choice or the specific assignment of the orbital positions.Note also that the natural orbital positions respect the symmetries of the lattice, more specifically the reflection symmetry with respect to the middle point of the line connecting two nearest-neighbor orbitals.We will see more examples of this general phenomenon in the following.
Besides the superfluid weight, another potentially interesting application of the natural orbital positions is related to topological invariants.The SSH model possess a topological invariant, the winding number W [65][66][67], which takes value W = 0 for v > w and |W| = 1 for v < w when computed using the projector operator given in (101).In fact, one can show that the quantum metric is bounded from below by the winding number [68], namely M > W 2 /2.However, it is well-known that the winding number depends on the choice of unit cell used to compute it, which amounts to a choice of the orbitals positions [69,70].Interestingly, our results show that the winding number of the SSH model computed with the natural orbital positions is always zero.Our next example, the Creutz ladder, possesses bands with nonzero winding number, when computed using the natural orbital positions.This suggests that the natural orbital positions can be used to provide a more refined classifications of the topological properties of the band structure.

B. Creutz ladder
The Creutz ladder shown in Fig. 2 is a one-dimensional lattice model with two orbitals per unit cell introduced in Ref. [71].It has the peculiar property that both of its two bands are perfectly flat for a specific choice of the hopping matrix elements.This model has been studied extensively both in the bosonic and the fermionic case with the inclusion of different types of interaction terms [48,55,68,[72][73][74][75][76][77][78].
Adopting the convention for the phases of hopping matrix elements introduced in Ref. [55], the free Hamiltonian of the Creutz ladder reads with t > 0 the energy scale of the hopping amplitudes.This Hamiltonian is represented graphically in Fig. 2. Choosing the orbital positions as for the SSH model (98), i.e. same position for orbitals in the same unit cell, the Fourier transform of the free Hamiltonian of the Creutz ladder takes the form The dispersions of the two bands, labeled by n = ±, and the associated projection operators are given by As in the case of the SSH model, the spin index has been dropped since the projection operator is the same for spin up and spin down.Again, the uniform pairing conditions is satisfied for From (113) it is immediate to obtain the result thus the components of the vector s in ( 87) are all zero.It follows that the orbital positions given by ( 98) are in fact natural ones for the Creutz ladder and the superfluid weight is obtained simply from the integrated quantum metric M, while the GRPA correction vanishes.As noted in Ref. [68], the quantum metric of the Creutz ladder is constant, independent of k, For completeness, we compute also the matrix R (88) and its pseudoinverse Note that the natural orbital positions respect the symmetries of the Creutz ladder, in this specific case they are preserved under the interchange of the two orbitals inside the unit cell.The unitary operator R implementing this transformation on the field operators is Rĉ The sign factors (−1) j in the definition of R correspond a gauge transformation and are necessary in order to preserve the signs of the hopping matrix elements of the Creutz ladder Hamiltonian (see Fig. 2).
As mentioned in the previous section, the winding number computed with the natural orbital positions is W = 1 [68], in contrast to the SSH model, for which it is always zero.It is an interesting open question for the future is to understand whether the different winding numbers of the two models manifest in some observable properties, for instance in the edge states that occur in topologically nontrivial lattice models due to the bulkedge correspondence.

C. Dice lattice
In this final example, we explicitly calculate the set of natural orbital positions and the minimal quantum metric for the two-dimensional dice lattice, also known as T 3 lattice [55, 79,80].The graphical representation of the free Hamiltonian of the dice lattice is given in Fig. 3, while its Fourier transform is provided in Ref. [55] and is not repeated here for brevity.The labeling of the orbitals shown in the figure is the same as in Ref. [81] and different from Ref. [55].Note that in Ref. [55] a more general model is considered, whereas here we specialize to the case in which all the hopping matrix elements are equal up to a sign.The hopping sign is denoted by the color of the bonds in Fig. 3.As discussed in Ref. [79], with this specific choice of hopping signs, the band structure of the dice lattice is composed of six doubly degenerated flat bands.We focus on the lowest pair of degenerate flat bands (n = 1, 2) with energy The parameter ε h is the on-site energy of the hub sites, while the rim sites have zero on-site energy (see Fig. 3 for the definition of hub and rim sites).The unit of energy in the above expression is the absolute value of the hopping amplitude between any two sites.The periodic Bloch functions for this pair of flat bands are where the normalizing constant is given by c = (ε 2 + 6) −1/2 and we have introduced the quantities k i = k • a i , with a i=1,2 the fundamental vectors of the Bravais lattice The periodic Bloch functions in ( 119) and ( 120) are obtained with the following choice of the orbital positions i.e. the same position vector is assigned to all of the orbitals inside the rectangular unit cell shown in Fig. 3 since r iα is independent of the orbital index α.Note that these are different from the positions of the lattice sites actually used in Fig. 3.
A convenient feature of the dice lattice is that the components of the periodic Bloch functions and therefore also of the projection operator P (k) = |g k,1 ⟩⟨g k,1 |+|g k,2 ⟩⟨g k,2 | are polynomials in e ik1 and e ik2 , which allows to obtaining analytical results.Indeed, the calculation of the matrix R (88) and the vector s l (86)-( 87) is straightfoward but tedious and is best done with the help of a computer algebra system.The explicit expressions of these quantities are not particularly illuminating and are not provided here.We simply note that the following vectors give a solution of the linear system (90) for arbitrary values of the parameter ε h .In fact, they give the unique solution satisfying the constrain e • b l = 0, thus this is the solution obtained by using the pseudoinverse of the matrix R, as explained in Sec.V. Due to this constrain, the vectors b l give the natural orbital positions for the dice lattice in a coordinate system in which the origin is the baricenter of the orbitals inside a unit cell.The positions of the lattice sites shown in Fig. 3 coincide with the natural orbital positions given by ( 123)-(124).One can see once again that the natural orbital positions provide a maximally symmetric arrangement of the lattice sites.It also straightforward to compute the minimal quantum metric (84), which is proportional to the identity and whose diagonal elements are (l = x, y) The superfluid weight is proportional to the minimal quantum metric according to (83), provided the uniform pairing condition (42), or equivalently (C8), is satisfied.For a discussion of the uniform pairing condition in the case of the dice lattice see Ref. [81].Note that the integrated quantum metric M (85) of the dice lattice is in general not proportional to the identity matrix when computed using orbitals positions that are different from the natural ones.This means that using the integrated quantum metric alone in ( 83) can lead to an unphysical anisotropy of the superfluid weight.Thus, it is important to include the geometric GRPA correction, which amounts to using the natural orbitals positions, in order to restore the proper symmetry of the superfluid weight tensor.

VII. DISCUSSION AND CONCLUSION
In this work, we have provided the analytic expression for the superfluid weight within the GRPA in the isolated band limit.By analytic we mean that the evaluation of the superfluid weight is reduced to performing Brillouin zone integrals of certain combinations of the band dispersions ε nk and the band wave functions |g nk ⟩ of a generic lattice model and their derivatives with respect to quasimomentum.This allows to relate the superfluid weight, which is an important observable for superfluid systems, to the properties of the band structure, in particular the effective mass and invariants such as the quantum metric (52).Our results hold under specific assumptions: the interaction is of the Hubbard form (13) with negative coupling constants U α < 0 (attractive) that can depend on the orbital α, and the free Hamiltonian (1) is invariant under spin rotations along a given axis and also time-reversal symmetric for A = 0.Moreover, it is assumed that the uniform pairing condition is satisfied, namely that the pairing potential ∆ α is independent of the orbital index α (42).Under the same assumptions, it was shown in Refs.[16] and [18] that the superfluid weight computed within the mean-field approximation (BCS theory), denoted here by D s,c can be written as in (50) or, equivalently, as in (55) in terms of the effective mass, while the geometric contribution D (0) s,g is a weighted averaged of the quantum metric over the Brillouin zone (51).The geometric contribution becomes important for bands with small bandwidth compared to the interaction coupling constants U α since the conventional contribution vanishes in the flat band limit.
A major result of this work is that the same separation holds also at the level of the GRPA.Indeed, we find that the correction term to the superfluid obtained from the GRPA D s,g is also the sum of a conventional contribution D s,c and a geometric contribution D (1) s,g .The expression for the conventional part of the GRPA correction is given by combining ( 65), ( 71) and (73).One can see from the expression (65) for correlation functions of the form ( Ĵl , Nασ ) that D s,c vanishes in the flat band limit, as in the case of D (0) s,c .We also show that the conventional GRPA correction is a positive semidefinite tensor, meaning that this correction term leads to an increase of the superfluid weight compared to the mean-field result.
On the other hand, the geometric part of the GRPA correction, given by ( 66), ( 72) and ( 76), is not necessarily zero in the flat band limit and, as discussed in Sec.V and through examples in Sec.VI, it is important to include it in order to restore the geometry independence of the superfluid weight.Indeed, we find that, within the GRPA, the superfluid weight in the flat band limit is proportional to the minimal quantum metric M (84)-( 88), the integral of the quantum metric over the Brillouin zone minimized with respect to the orbital positions (see Sec. V).The relation between minimal quantum metric and superfluid weight was pointed out in Ref. [47] for the first time, where it was shown that the geometry independence of the superfluid weight is restored by not neglecting the A-dependence of the pairing potential ∆(A) when taking derivatives of the mean-field free energy (25).
In this work we extend the results of Ref. [47] in several ways.First, in (25) we take also into account the A-dependence of the Hartree-Fock potential Γ(A), which amounts to computing the superfluid weight within the full GRPA [9].By doing so, we obtain the conventional part of the GRPA correction to the superfluid weight, which is a new result.This does not affect the results of Ref. [47] since D s,c = 0 in the flat band limit, which means that the relation between superfluid weight and minimal quantum metric holds at the level of the GRPA.This is consistent with numerical investigations performed using quantum Monte Carlo methods [82,83], which find that (83) gives a rather accurate estimate of the superfluid weight in the flat band limit.
Another useful result of the present work is the derivation of the simple analytical expressions for both the minimal quantum metric (84)-( 88) and the natural orbitals positions, the latter given as the solution of the linear system (90).Our approach is based on the direct evaluation of the GRPA formula for the superfluid weight.Compared to the one of Ref. [47], where these band structure invariants were introduced, we are able to remove unnecessary assumptions, namely the flat band requirement and the condition (93).This is very important in view of potential future applications in other contexts.While the quantum metric has found by now many important applications, the issue of its geometry dependence is almost never addressed, therefore we expect the analytical formulation of the minimal quantum metric and the natural orbital positions to become very useful in this sense.Also, our gauge symmetry-based argument for establishing the geometry independence of the superfluid weight (Sec.II) may be extended to other observable quantities.
In Sec.VI, an interesting application is already provided, since we find that the SSH model is topologically trivial when the winding number is computed using the natural orbital positions, while the Creutz ladder is not.For the future, it would be interesting to better understand in what sense these two models, and other ones as well, are topologically distinct.Note that the SSH model is often used as a toy model to illustrate the concepts of bulk-edge correspondence and topological phase transition [65], while at the same time it is debated whether it is a truly topologically lattice model given that its winding number is not unit cell consistent [70].In view of these examples, we speculate that the concept of natural orbital positions may ultimately lead to a more refined classification of the topological properties of the band structure.
As we have seen in Sec.VI (see also Ref. [47]), another useful feature of the natural orbital positions is to provide a set of positions that are maximally symmetric without requiring any input other than the band projector P (k).Thus, they could find applications in the field of electronic structure theory, for instance.In this sense, more work is needed in order to better understand the physical meaning of the natural orbital positions since symmetry alone is not sufficient to determine them.For instance, this is the case of the Lieb lattice with staggered hopping [47], in which many orbital positions are compatible with the lattice symmetries, but the natural ones are uniquely determined up to translations.
Finally, an interesting direction for future work is to extend our results to bosonic superfluids.Some work as already been done in this direction [33][34][35]84], however the superfluid weight has not been computed in the full GRPA approximation since the A-dependence of the Hartree-Fock potential is not taken into account.This might be especially important in cases where translational symmetry is spontaneously broken by interactions, for instance the supersolid phase observed in ultracold gases with dipolar interactions [85][86][87][88].
Then, by minimizing the mean-field grand potential [9], one obtains the self-consistency equations of mean-field theory The variational coefficients Γ σ j,j corresponds to the Hartree potential, Γ σ i,j for i ̸ = j to the non-local Fock potential and ∆ i,j to the pairing potential.
Following the notation of Ref. [9], we define the vector of mean-field potentials as and we use the following variation of the Einstein summation convention ∂g .

(A10)
The vector of expectation values of operators quadratic in ĉjσ , ĉ † jσ is then written as (see (A4)-(A5)) Using this notation we can write the expectation value of the interaction term in a concise form Thus, from this last equation and (A3), we have for the the mean-field grand potential Taking the partial derivative with respect to ∆ c gives The first partial derivatives of the mean-field grand potential vanish when These are just the self-consistency equations of meanfield theory (A6)-(A8) written using our Einstein summation convention (A10)-(A12).After taking another partial derivative of (A14) and imposing that the meanfield potentials are self-consistent (A15), we have Recall that both the mean-field grand potential Ω m.f. and Ω 0 depend either directly or indirectly on the vector potential A, namely Ω m.f.A, ∆ a (A) and Ω 0 A, ∆ a (A) (see (25)).The self-consistency equations are satisfied for any A, which means that This last result gives the derivatives of the mean-field potentials with respect to A Here and in the following, we use the notation [M a,b ] −1 to denote the inverse of the matrix with elements M a,b rather than the inverse of a single matrix element.We ignore here the subtleties occurring when the inverse of the Hessian matrix ∂ 2 Ω m.f.∂∆ a ∂∆ b does not exists.In this case one should express the solution of the linear system (A17) using the pseudoinverse (Moore-Penrose inverse).
We can now compute the full derivatives of the meanfield grand potential and express them in terms of the mean-field solution for A = 0. We start with the first full derivative where in the last equality we have used (A14)-(A15).The second full derivative is then ∂∆ a ∂A l . (A20) where in the last two equalities we have used (A17) and (A18).Since the self-consistent solution minimizes Ω m.f., the Hessian matrix ∂ 2 Ω m.f.∂∆ a ∂∆ b (and also its inverse) is positive semidefinite.This fact together with (A20) implies that the superfluid weight is bounded from above (in the sense of matrix inequalities) by the matrix ∂ 2 Ω m.f.∂A m ∂A l .
It is convenient to express the result in (A20) only in terms of the quantities ∂ 2 Ω 0 ∂∆ a ∂∆ b , which are correlations function relative to the mean-field statistical ensemble, as shown below.We have already done this is in the case of ∂ 2 Ω m.f.∂∆ a ∂∆ b in (A16), thus we only need to do the same for ∂ 2 Ω m.f.∂A m ∂A l and ∂ 2 Ω m.f.∂∆ a ∂A l .Since Ω m.f.= Ω 0 + ⟨ Ĥ − Ĥ0 ⟩ we consider the partial derivatives of the expectation value Again this quantity vanishes if the self-consistency equations (A15) are satisfied, therefore from (A19), we have where (6) has been used in the last equality.Thus, we find the important result that the first full derivative of the grand potential with respect to A is proportional to the current.Taking another partial derivative with respect to A m in (A21) and (A14) and imposing selfconsistency gives Using these results together with (A16) and the last line of (A20) leads to To simplify this expression we use the identity where A and B are matrices defined by Thus, our final result is With further simple manipulations, it is possible to prove the equivalence between this and the result for the superfluid weight in the generalized random phase approximation provided in Ref. [9].We can specialize (A12) to the case of the attractive Hubbard interaction given by Ĥint = − j U j nj↑ nj↓ , with U j > 0 . (A29) The Hubbard interaction in the above equation is more general than (13) since the latter is translationally invariant, while the former is not as we allow for an arbitrary dependence of the coupling constant U j on the site index j.The Hubbard interaction in (A29) is obtained by choosing the interaction coefficients in (A1) as follows Then, the expectation value of the interaction term in (A12) becomes This gives the matrix B in (33), after translational invariance is enforced.Instead, the expressions for the matrices A α,β (31) and the vectors v lα (30) in terms of correlation functions of the form (26) are obtained by using the results [9] and respectively.As an example, we have In the second line we have performed the substitution i → iα and used the fact that the correlation function Ĵl , ĉ † iα↑ ĉ † iα↓ is independent of the unit cell index i since the operator Ĵl is translationally invariant.This completes the derivation of the expression for the superfluid weight in the generalized random phase approximation given in (35) in the main text.

Appendix B: Nambu formalism and evaluation of correlation functions
In this section, it is explained how to evaluate the correlation function (26) between translationally invariant quadratic operators, which are conveniently recast in the following Nambu form Here ĉkσ (ĉ † kσ ) is the column (row) vector whose components are the field operators ĉkασ (ĉ † kασ ) for α = 1, . . ., N orb .The column vector ĉk in (B4), grouping together both creation and annihilation operators, is called a Nambu spinor.Consistently with the Nambu spinor structure, the quadratic operatic Â(k) corresponds to the single-particle operator A(k), a 2N orb × 2N orb matrix consisting of four blocks of dimension N orb denoted by [A(k)] i,j , as shown in (B3).We adopt the convention that the blocks of a single-particle operator in the Nambu representation are labeled by superscripts, while subscripts label the matrix elements in each block.This convention is used already in (B2).Note that, in order to bring a quadratic operator in Nambu form, it is necessary to rearrange the field operators, which produces additional c-number terms due to the fermionic anticommutation relations.However, all of the c-number terms cancel out when the expectation value of the same operator is subtracted, namely Â(k) − ⟨ Â(k)⟩.This combination is precisely the one appearing in (26), implying that, for the purpose of computing correlation functions, we are free to reorder the field operators and represent quadratic operators in Nambu form.The ultimate reason for using the Nambu formalism is that it allows to diagonalize in a convenient way quadratic operators that contain anomalous terms, such as ĉiα↓ ĉiα↑ and ĉ † iα↑ ĉ † iα↓ .In our specific case, we need to diagonalize the variational Hamiltonian Ĥ0 (see Sec. IV).
In the following we denote by ⟨ĉ k (τ )ĉ k ′ ⟩ the set of expectation values obtained by replacing each of the Nambu spinors with any of their components, namely . ., N orb .The same convention is used for ⟨ĉ † k (τ )ĉ k ′ ⟩ and so on.Using this notation, we can express in a concise way the constraints imposed by the conservation of spin Ŝz and momentum on the expectation values of products of two operators, which read These relations are used in the evaluation of the expectation value that appears under the integral sign in ( 26) Here T τ is the time-ordering symbol for imaginary time where â and b are anticommuting operators.Note that 0 < τ < β in (26) and the field operator are ordered accordingly in (B8).In the third equality, we have used Wick's theorem [60], which holds since all expectation values are evaluated with respect to the quadratic Hamiltonian Ĥ0 , see (22).In the last line of (B8) we have introduced the standard definition of the imaginary-time Green's function After expanding the Green's function using Matsubara frequencies ω n = (2n + 1)π/β in (B8) and performing the imaginary time integral in (26), one obtains (B12) The Nambu form H 0 (k) of the variational Hamiltonian is given in (39).It is a standard result that the Matsubara Green's function can be written as [60] G The second line follows from (43).To proceed, we also need the operators Ĵl , Nα and Dα in Nambu form, namely It is important to keep track of the minus sign associated with the down spin in (B17).
In order to perform the Matsubara frequency summation in (B12), it is convenient to introduce the following matrix (B21) This is simply the Green's function G(k, iω n ) (B13) from which the Bloch functions U k have been removed.Note that each of the blocks of L(k, iω) is diagonal, namely [L(k, iω)] i,j m,n = [L(k, iω)] i,j n,n δ m,n , and the diagonal elements are given by where i ̸ = j in the last equation.Using these definitions and the Matsubara frequency sum with n F (E) = (e βE + 1) −1 the Fermi-Dirac distribution, we obtain the following useful results (i ̸ = j) As an example, it is shown how to compute the correlation function ( Dα , D † β ) using the above results.From (B18)-(B20) and (B12) Finally the Matsubara sum is evaluated with (B27).From (B26)-(B29), it is apparent that several relations hold between the correlation functions ( Â, B), where the operators Â and B are taken from the set { Nασ , Dα , D † α }.They are the following In (B32) and (B33) σ denotes the spin opposite to σ.
From the properties ( 27)-( 28), we obtain also the additional relation If n denotes the only partially filled band, then in the isolated band limit E nk ∼ U and E nk ∼ E gap for n ̸ = n and the leading order contribution ∼ U −1 to the matrix A α,β is obtained by retaining only the term n = m = n in (B30) and in all of the other correlation functions appearing in (31).Thus, the Matsubara sum (B25) becomes Observe that, in the isolated band limit, all the matrix elements of A α,β are real since the Matsubara sums (B26)-(B29) are always real and the Bloch functions of the band n appear in the combination , which is positive.As a consequence, the correlation functions that appear in (B33) and (B35) are all equal.With the notation and the results established so far, it is immediate to obtain the matrix elements of A α,β in the isolated flat band limit, which are shown in ( 57)-( 62).The correlation functions that involve the current operator Ĵl require some care when taking the isolated flat band limit.To begin with, time-reversal symmetry implies that the current operator in Nambu form becomes To derive (B39) from (B14), we have used the identity Using this last result, one can show that the diamagnetic term in (37) can written as a correlation function in the same way as the paramagnetic one, that is The fact that the diamagnetic term is nonnegative is an immediate consequence of the general property (29).
In order to compute Ĵl τ z , Ĵm τ z and other correlation functions that involve the current operator, it is useful to introduce the following operator in Nambu form It is not difficult to show that this operator satisfies the following properties with From these relations and the results in (B26)-(B29), one obtains the following expressions for the matrix elements of Using the definitions in ( 45)-( 47) and after some laborious algebra, one finds the following alternative expressions for the same matrix elements (B50)-(B51) These are useful to compute the components of the vector v lα in (30), which read The result for [M − l (k)] 1,1 n,n is also useful and is computed in a similar way The easiest way to obtain the second equality is to take the limit ε mk − ε nk → 0 in (B52).
In the case m ̸ = n, the term ∝ ∆ 2 in (B53) can be ignored from the start, thus [M + (k)] 1,1 m,n ≈ [M − (k)] 1,1 m,n .Indeed, this term gives a contribution of order ∆ 2 /E gap ∼ U 2 α /E gap , which vanishes in the isolated band limit.Then, using again (B58), we derive the following useful result (B61) In the isolated band limit, only the term n = n corresponding to the partially filled band gives a nonzero contribution.Note that, to obtain the correct result for the correlation function (B61), it is important to retain all of the matrix elements of the current operator [ Jl (k)] m,n , even if n, m ̸ = n.Indeed, retaining these matrix elements allowed us to use the completeness relation for the Bloch functions in (B60).The need to take into account also the interband matrix elements of the current operator even in the isolated flat band limit is a general phenomenon, as we will see in the following.For this reason, one has to be particularly careful when evaluating correlation functions that involve the current operator.
To compute the correlation function between the current operator and the pairing operator (B56), one can use the approximation In fact, according to (37), only their difference is required × ⟨∂ l1 g mk |g nk ⟩ ⟨g nk |∂ l2 g mk ⟩ . (B67) In the isolated band limit, the first sum gives the conventional contribution to the superfluid weight in (50), while the second one with m ̸ = n corresponds to the geometric contribution in (51).Indeed, in the case of the geometric contribution, one can proceed as follows Thus, the Bloch function |g nk ⟩ of the partially filled band enters only through the quantum metric (52).Again, it is important to notice that the completeness relation has been used in the above derivation.
Appendix C: Self-consistency equations of mean-field theory for the Hubbard interaction and uniform pairing assumption In this section, we solve the self-consistency equations of mean-field theory in the case of the Hubbard interaction, thus justifying the uniform pairing condition (42).In order to solve the self-consistency equations ( 23) and (24), it is necessary to compute the expectation values that appear on the right hand side.These are obtained from the imaginary-time Green's function (B10) by taking the limit τ − τ ′ → 0 − Using the definition in (B21), we obtain Recall that ε k = diag(ε nk ) is a diagonal matrix containing the band dispersions ε nk , while E > k = diag(E nk ) is also diagonal but contains the quasiparticle dispersions E nk instead, that is the eigenvalues of the BdG Hamiltonian H 0 (k), see (43).Note that the result in (C5) is valid only under the uniform pairing condition (42).Then, the Green's function at τ = 0 − is computed by combining (B13) with (C5).Thus, from (C1)-(C3), we can rewrite the self-consistency equations as The parameters ∆ α , obtained from the second equation for a given value of ∆, do not in general satisfy the uniform pairing condition.However, it is possible to adjust the relative strength of the coupling constants U α so as to ensure that ( 42) is at least approximately satisfied.
In the isolated band limit, it is possible to derived an explicit condition on the coupling constants U α that ensures uniform pairing.First, the self-consistency equation relative to the Hartree potential Γ σ α is neglected for simplicity and only the partially filled band n is retained in the sum over bands in (C7).Indeed, the contribution of the terms n ̸ = n is negligible in the isolated band limit since ∆ is of the same order of U α .If it is assumed that the n-th band is flat (ε nk = ε n), then the quasiparticle dispersion E nk = E n = (ε n − µ) 2 + ∆ 2 is also flat and the uniform pairing conditions is equivalent to the following requirement namely that the quantity on the left hand side is independent of the orbital index α, when it is not zero.(C8) can always be satisfied by a suitable choice of the coupling constants U α .In this case, the self-consistent value of the pairing potential ∆ (called also the order parameter) is obtained from This equation is obtained by combining (C7) and (C8) and is identical to the self-consistency equation of the Weiss mean-field theory of ferromagnetism, in which the quasiparticle energy E n plays the role of the magnetization.If the quantity k |⟨α|g nk ⟩| 2 is independent of the orbital index α, for instance because of some lattice symmetry, then the uniform pairing condition follows from U α = U β = U = N orb Ū for all α, β, where N orb is the number of orbitals in the unit cell.This is the case considered in Ref. [16] and other subsequent works.
is also the sum of a conventional contribution D , as in the case of the known mean-field result D

FIG. 2 :
FIG.2:The Creutz ladder is a one-dimensional lattice model composed of two simple linear chains.The two chains are distinguished by the orbitals index α = 1, 2 and coupled by inter-chain hopping matrix elements, shown as the red and black diagonal lines in the figure.The horizontal and diagonal black lines correspond to the hopping amplitude t > 0 in the free Hamiltonian (109) and the red lines to −t.With this choice of the hopping matrix elements, the two bands of the Creutz ladder are perfectly flat and geometrically non-trivial since the minimal quantum metric is nonzero.As in the case of the SSH, orbitals in the same unit cell (green rectangle) have the same position.As shown in the main text, this choice of the orbital position is natural.The two chains are displaced from each other in the transverse direction only for illustration purposes.

FIG. 3 :
FIG.3: Schematics of the dice lattice.The orange hexagons and the blue triangles denotes lattice sites.The hexagons are called hub sites and are six-fold coordinated while the triangles are called rim sites and are three-fold coordinated.The bonds between sites represent Hamiltonian matrix elements, which are all equal up to a sign encoded in the color of the bond (+ for black and − for orange).With this choice of the hopping matrix elements, the Bravais lattice is rectangular with the fundamental vectors given in (121).The red rectangle denotes the associated unit cell, which contains six orbitals, labelled by the orbital index α = 1, . . ., 6 as shown in the figure.The black cross denotes the baricenter of the orbitals within the chosen unit cell.The natural orbital positions are specified by the vectors bx (123) and by (124) with the baricenter taken as the origin of the coordinate system.The positions of the lattice sites in the figure coincide with the natural orbital positions.