Superfluid density and collective modes of fermion superfluid in dice lattice

The superfluid properties of attractive Hubbard model in dice lattice are investigated. It is found that three superfluid order parameters increase as the interaction increases. When the filling factor falls into the flat band, due to the infinite large density of states, the resultant superfluid order parameters are proportional to interaction strength, which is in striking contrast with the exponentially small counterparts in usual superfluid (or superconductor). When the interaction is weak, and the filling factor is near the bottom of the lowest band (or the top of highest band), the superfluid density is determined by the effective mass of the lowest (or highest) single-particle band. When the interaction is strong and filling factor is small, the superfluid density is inversely proportional to interaction strength, which is related to effective mass of tightly bound pairs. In the strong interaction limit and finite filling, the asymptotic behaviors of superfluid density can be captured by a parabolic function of filling factor. Furthermore, when the filling is in flat band, the superfluid density shows a logarithmic singularity as the interaction approaches zero. In addition, there exist three undamped collective modes for strong interactions. The lowest excitation is gapless phonon, which is characterized by the total density oscillations. The two others are gapped Leggett modes, which correspond relative density fluctuations between sublattices. The collective modes are also reflected in the two-particle spectral functions by sharp peaks. Furthermore, it is found that the two-particle spectral functions satisfy an exact sum-rule, which is directly related to the filling factor (or density of particle). The sum-rule of the spectral functions may be useful to distinguish between the hole-doped and particle-doped superfluid (superconductor) in experiments.

For a single flat band system, due to the divergence of effective mass 17 , it is believed that the flat system usually has no superfluidity, e.g., the superfluid density (which is superfluid weight up to a constant) ρ s = 0 . However, for a multi-band system, the situation may be quite different 18,19 . In contrast to the single-band system, the superfluid density in a multi-band system can be divided into two part contributions. One of that consists of the diagonal matrix elements of current operator (conventional part), and the other one is the off-diagonal terms (the so called geometric part [20][21][22] ). Even through the conventional part may vanish in flat band limit, the geometric one may be finite. Furthermore, it is found that under some conditions, the geometric part of superfluid density is related to the geometric quantity of energy band, namely, the geometric metric tensor of Bloch states 23 .
A flat band can appear in dice ( T 3 ) lattices, similar as graphene, which also has the honeycomb lattice structure (see Fig. 1). In comparing with the graphene, there is an extra lattice site (per unit cell) in the center of hexagon of honeycomb lattice [24][25][26] . The low-energy physics of the dice lattice are also described by the Dirac equation, but its pseudospin S = 1 instead of 1

Results
The structure of dice lattices is depicted in Fig. 1, every unit cell consists of three different lattice sites (A, B and C), in which sites A and B are located at the vertices of hexagon to form honeycomb lattice, and site C is located at the center of the hexagon. In this work, we only consider the nearest neighbor hopping and equal hopping amplitude, i.e., t 1 = t 2 = t (see Fig. 1). The Hubbard model in which two-component Fermi gases with spin σ =↑, ↓ interacting through on-site attractive interaction −U ( U > 0 ) is where µ describes the chemical potential. In the following, we would use i to label the unit cell. The creation operators for three sub-lattices are represented by ψ † iA , ψ † iB and ψ † iC , respectively. n iσ = ψ † iσ ψ iσ is the particle number operator. Under periodic boundary condition , using a Fourier transform, it is convenient to obtain the non-interacting Hamiltonian with ψ † (k) = [ψ † A (k), ψ † B (k), ψ † C (k)] and www.nature.com/scientificreports/ where h 12 (k) = t δ cos(k · a δ ) + it δ sin(k · a δ ) , h 21 (k) = h * 12 (k) , h 23 (k) = h 12 (k) , h 32 (k) = h 21 (k) , and k = (k x , k y ) . The two dimensional vectors a δ are shown in Fig. 1. We note that it is a multi-band (threeband) system and there exist a flat band [ ε 0 (k) = 0 ] in between the highest [ ε + (k) ] and the lowest bands [ ε − (k) = −ε + (k) ] (see Fig. 1).
It is worth noting that an important role played by particle-hole symmetry in this model. Similarly as Haldane-Hubbard model 40 , applying a particle-hole transformation to the system, where ε i = 1 for the A and C sublattices, and for the B sublattices ε i = −1 , the hopping term in Hamiltonian Eq. (1) is invariant under the above transformation. After adding some constant terms, the interaction term remains invariant and chemical potential term changes a sign. At this point, we only need to consider the case where the number of particles n > 3 ( n = 3 is half-filling, i.e. there are three particles in a unit cell). For filling n ≡ n − 3 ≥ 0 (particle-doped case) and − n = 3 − n ≤ 0 (hole-doped case), the chemical potentials satisfy a relationship It indicates that, at half-filling ( n = 3 ), the chemical potential is exactly zero. Furthermore, we will see that all the other physical quantities, e.g., superfluid order parameters, superfluid density, etc., are also symmetrical with respect to the half-filling n = 3 (see Fig. 2).
The order parameters ( � i = −U�ψ i↓ ψ i↑ � ) can be calculated with mean field theory (see "Methods"). The evolutions of the superfluid order parameters are depicted in panel (a) of Fig. 2. It is found that i is always real. The order parameters are symmetrical with respect to the half-filling due to the particle-hole symmetry. The Hamiltonian is subject to the inversion symmetry, A is the same as C but differ from B . As one increases the interaction U, i always increases. The order parameters can increase without limit as long as the interaction is strong enough. For example, when filling factor is half-filling ( n = 3 ), � i ≃ U/2 as U → ∞ , which is a half of the energy of two particle bound states 41 .
When the filling factor approaches zero, the order parameters decrease to zero. The fully occupied case ( n = 6 ) is equivalent to completely empty ( n = 0 ) case due to the particle-hole symmetry. With the increasing of filling factor n, the i gets larger. However, when the filling factor is in between the lowest band and the middle flat band, the order parameter gets a dip due to the small density of states near the Dirac points [especially for weakly interacting cases in panel (a) of Fig. 2]. The suppressions of order parameters are also found in Lieb lattices 42 . When the filling factor lies in the middle region (e.g., n ≃ 3 ), the order parameter reaches its maximum due to enhancement of density of states in flat band.   The superfluid densities for interaction U = 1 , U = 2 , U = 3 and U = 4 are also shown in blue, green, red and black solid lines, respectively.  43 and ferromagnetic transition 44 . Similarly, the behaviors of superfluid order parameters are also affected significantly by the flat band 45 . For example, when the filling factor is in the middle of the flat band (e.g., half filling n = 3 ), the infinite large density of states results in order parameter � A = � C ≃ U/4 as the interaction approaches zero, which is very different from that in the usual Bardeen-Cooper-Schrieffer (BCS) superconductor (or superfluid), where pairing gap is exponentially small when interaction U → 0 . In addition, for half-filling n = 3 , it is found that the order parameter for sublattice B is much smaller than � A(C) when the interaction is small ( U/t ≪ 1 ), i.e., � B /� A(C) ≪ 1 (see Fig. 2). This is because the weights of sublattices A and C in flat band wave function is dominant over that of sublattice B. Superfluid density. The superfluid density can be calculated with phase twist method 20,46 . Assuming the superfluid order parameters undergo a phase variation, e.g, △ i → i e 2iq·r i , the superfluid density (particle number per unit cell) tensor ρ sij can be written as where is thermodynamical potential (per unit cell) in grand canonical ensemble. Similarly as the Haldane model 20 , due to the C 3 symmetry of honeycomb lattices 19,47 , the superfluid density tensor can be simplified into a scalar, i.e., ρ sij = diag{ρ s , ρ s }.
The superfluid density ρ s is plotted as a function of filling factor for different interactions in panel (b) of Fig. 2. First of all, we see that the superfluid density is also symmetrical with respect to the half-filling due to the particle-hole symmetry. When the filling factor is completely empty ( n → 0 ) [or fully occupied case ( n → 6)], the superfluid density vanishes linearly. For the small filling factor ( n ≪ 1 ) and weak interaction ( U ≪ t ) case (BCS-limit 48 ), the superfluid density can be written as where m * = √ 2/(3t) is the effective mass near the bottom of lowest band, which is determined by the singleparticle energy band structure, and does not depend on interaction strength U. The superfluid density is proportional to filling factor (or particle density) and inversely proportional to the effective mass of energy band 17 , which is very similar to that of Bose-Einstein condensate with spin-orbital coupling 49 .
Differently from order parameters, when interaction U/t ≫ 1 , the superfluid density remains finite. In such case, the superfluid consists of tightly bound pairs. Two bound pairs could not occupy a same lattice site due to the constraint of exclusive principle of fermions. So the pairs can be viewed as hard-core bosons, which is BEC limit of fermion superfluid. When the number of the effective hard-core bosons is small n ≪ 1 (dilute gas limit), the superfluid density is determined by the effective mass of tightly bound pairs, i.e., where n ′ = n/2 is filling factor of effective hard core bosons and m ′ * = U/(16t 2 ) is its effective mass, which is determined by the dispersion of two-particle bound state energy (Note m ′ * � = 2m * due to the coupling of between the motions of center of mass and relative motion of two particles in lattice 50 ). The above equation indicates that when the interaction gets stronger (in strong interaction limit), the effective mass becomes larger, the resultant superfluid density gets smaller. For a moderate interaction strength (e.g., U/t > 5 ), and small filling ( n ≪ 1 ), it is expected that the superfluid density should satisfy 8t 2 n/U < ρ s < 2.12nt.
In addition, differently from the order parameters, the evolution trends of superfluid density with the interactions can bedifferent for different regimes of filling factors. For example, when the filling is small, superfluid density diminishes with the increasing of interaction, which is opposite to that of order parameters [see panel (b) of Fig. 2]. Only when the filling is sufficiently large (for example, n ≃ 3 ), the superfluid density grows up as the interaction increases (within a certain range of interaction, see the following discussions and the panel (a) of Fig. 3). Due to the vanishingly small density of states near the Dirac points, the superfluid density also develops a dip for weak interacting cases. Similarly as Lieb lattices 18 , a triple dome structure appears in the superfluid density, especially for weakly interacting cases [see panel (b) of Fig. 2].
For intermediate filling factor (e.g, n = 3 in flat band), it is found that the superfluid density approaches zero as U → 0 . Once the interaction is turned on, the superfluid density grows rapidly [see panel (a) of Fig. 3]. Specially here, the non-vanishing superfluid density for weak interactions mainly arises from the contributions of off-diagonal matrix elements of current operator (the geometric part). In the isolated flat band limit 18,22 or under the assumption of uniform pairing 1,20 , the superfluid density in flat band can be related to geometric metric tensor of Bloch sates. However, here we note that the superfluid density approaches zero according to the law of as U → 0 with two constants a ≃ 0.4077 and b = 9.822 for half filling n = 3 [see panel (a) in Fig.3], which is different from the linear dependence in isolated flat band limit 18 . As U → 0 , the appearance of the logarithmic singularity (UlnU) is attributed to the existences of the gapless Dirac points and vanishingly small pairing gap ( ∝ U → 0 ). This is because when filling factor is in the middle of flat band (half-filling of n = 3 ) and U → 0 (or � ≡ � A = � C ≃ U/4 → 0 ), it is found that the superfluid density mainly arises from matrix elements of current operator between the flat band and other bands near the Dirac points (the off-diagonal part). Furthermore, the contribution near the Dirac points can be approximated by With the increasing of interaction, the superfluid density reaches its maximum values at U ≃ 3.7t [see panel (a) in Fig. 3]. After passing over that point, it decreases. For strong interaction limit and finite n, the trend of evolutions of superfluid density can be also understood qualitatively by effective mass of tightly bound pairs. Considering the fermion nature of tightly bound pairs, it is expected that the superfluid density should be smaller than that determined by m ′ * , namely, ρ s < 8t 2 n/U for BEC limit ( U ≫ t ). For example, when n = 3 , it is found that the superfluid density is also inversely proportional to interaction, i.e., ρ s → C/U as U → ∞ [see panel (a) of Fig. 3]. The numerical results give the constant C ≃ 12 for half-filling case ( n = 3 ), which is almost a half of the estimated up bound [ 8t 2 n = 24 (t=1)]. Nevertheless, the up bound would provide much more precise estimation for superfluid density when filling factor n is small [see panel (b) of Fig. 3]. For a general filling n, it is found that the constant C can be given by a parabolic function (see "Methods"), namely Consequently, the superfluid density The asymptotic results Eqs. (11) and (12) in strong interaction limit are very similar to the behaviors of the superfluid density of hard-core bosons in lattices 52,53 . Collective modes. In this section, we would present Gaussian fluctuation method [54][55][56] to calculate the collective modes and spectral functions 40 . The order parameters can be decomposed as mean-field part and fluctuation, i.e., � i = � i + δ� i with i = 1, 2, . . . m , and m is the number of order parameters. The partition function can be written as where S 0 is the mean-field contribution and the Gaussian fluctuation part (see "Methods") where is matrix element of Nambu-Gorkov Green function. Here k = (k, iω n ′ ) and ω n ′ = (2n The collective modes are given by zeros of determinant Det|M(q, iω n → ω + i0 + )| = 0 . As q → 0 , the gapless collective mode is the Anderson-Bogoliubov phonon, which characterizes the density oscillations of superfluid. The gapped ones are Leggett modes, which correspond relative oscillations of densities of different sublattices [5][6][7][8]40 .
The two-particle normal (anomalous) Green function matrix elements are defined as 47 where ω n0 = E n − E 0 , E n and |n� are eigenvalues and eigenstates of Hamiltonian, respectively. iq is pairing fluctuation operator of superfluid order parameters [see Eq. (22) in the next section]. With the action δS (Gaussian weight), the correlation functions (and the matrix elements of two-particle Green functions) can be calculated 47 , i.e., The two-particle spectral functions can be expressed in terms of the Green functions as 57 where 0 + denotes a positive infinitesimal number. In addition, in the presence of inversion symmetry, the superfluid density can be also related to the pairing fluctuation matrix M 47 and (. . .) t denotes matrix transpose. In dice lattice case, the number of order parameters m = 3 and we identify 1 = A , 2 = B and 3 = C . Figures 4 and 5 show the collective modes with the increasing of wave vector q. It is found that there exist three undamped collective modes for strong interactions [see panels (a) and (b) of Fig. 5]. The lowest one is the gapless phonon as q → 0 , which corresponds to total density oscillation. The upper two gapped excitations are the Leggett modes, which corresponds to the relative density oscillations between sublattices. As the interaction strength decreases, the bottom of two-particle continuum lowers, the region of the existence of collective modes shrinks. The collective modes merge into the continuum and becomes damped [see panels (a) and (b) of Fig. 4]. As a result of that, the collective modes are not well defined. Only when the interaction is strong enough, the undamped Leggett modes survive (see Fig. 5). The collective modes also reflected in the two-particle spectral functions as shown in Figs. 4 and 5. The sharp peaks in spectral functions correspond the collective modes [see Two-particle spectral functions. The two-particle spectral functions of a general superfluid state can be also defined as     where ψ αk+q fermion field operator for single particle state α in momentum space, which satisfies anti-commutation relation The two-particle normal (anomalous) Green functions can be expressed in terms of spectral functions, i.e., Similarly as single-particle spectral function 58 , the two-particle spectral functions satisfies a sum-rule, i.e., Taking the diagonal elements ( γ = α; δ = β), where N cell is the number of unit cells, N α is the particle number of α-th component. In our dice lattice case, there exist three order parameters and they can be identified as � A(B/C)↓A(B/C)↑ = � A(B/C) . The sum-rule is where indices i, j = A, B, C , particle number of i-th sublattice N i = N i↑ + N i↓ .
If we take the diagonal elements of A ij and sum over i, then we get a much more useful sum-rule where total particle number N = N A + N B + N c and filling factor n = N/N cell . It should be emphasized that the above sum-rule of spectral functions is exact, which does not depend on what state the system is. One can use the above sum-rule to check the calculations of various theories. The above Eq. (28) shows that the integral of spectral functions is proportional to filling factor n − 3 (relative to the half filling). If the filling factor is half filling ( n − 3 = 0 ), the integral of spectral functions is exactly zero. While when filling factor is larger (smaller) than half-filling [particle (hole)-doped case], the integral of spectral functions is negative (positive). Our numerical calculations verify the above observations. For example, when the filling is half-filling ( n = 3 ), we see that the spectral functions are odd functions of frequency and the integral of spectral functions vanishes exactly [see panels (c, d) of Fig. 4 and panel (c) of Fig. 5]. Panel (d) of Fig. 5 shows the spectral functions for particle-doped case ( n = 5 ), the weight of negative frequency ω < 0 is dominant, and (21) A αβ,γ δ (q, ω) = n �0|� αβq |n��n|� † γ δq |0�δ(ω − ω n0 ) − n �0|� † γ δq |n��n|� αβq |0�δ(ω + ω n0 ), www.nature.com/scientificreports/ the integral of spectral functions is negative. The above findings may be useful to distinguish between the holedoped and particle-doped superfluid or superconductor in experiments. In addition, for a fixed q , when ω approaches some non-degenerate collective modes (e.g., phonon or Leggett modes), from the definition Eq. (21), the spectral functions should take the following forms where |n� = |q� is wave function of the collective mode, ω n0 is its excitation energy. As a result of that, when ω ≃ ω n0 , the heights of peaks of two-particle spectral functions near the collective mode should satisfy a relation, i.e., When q is small ( q → 0 ), the matrix elements of two-particle Green functions of zero frequency ω = 0 are proportional to product of two order parameters, namely 47 where ] ≃ 1/(q 2 ρ s ) > 0 is a real number and θ q the phase operator 59,60 . In such a case, it is expected that, for small q, the two-particle spectral functions near the phonon states [ ω ≃ ω phonon (q) or ω ≃ −ω phonon (−q) ] are also proportional to product of superfluid order parameters, i.e., where Z + = |�0|θ q |q�| 2 and Z − = |�0|θ −q | − q�| 2 are two real number, which are approximately equal as q → 0.
Our numerical calculations indeed verify Eqs. (30) and (32) for the case of q → 0 . For example, some numerical results of spectral functions are listed in Table 1. It is found that, for extremely small wave vector ( q = 0.01 ), both Eqs. (30) and (32) are satisfied approximately. As q increases from q = 0.01 to q = 0.1 and eventually q = 0.5 , Eq. (32) becomes invalid gradually (see the columns 6, 7, 8 in the Table 1). Nevertheless, Eq. (30) still holds (see the last column of the Table 1).

Discussion
In conclusion, we investigate the superfluid properties of attractive Hubbard model in dice lattice. The order parameters and superfluid density as functions of filling factors and interactions are analyzed in great detail. We emphasize the important roles of effective masses in the understanding of the asymptotic behaviors of superfluid density, especially for small filling factor. When both the filling factor and interaction are small, the superfluid density is given by the effective mass of single particle energy band. When interaction is strong, the superfluid density is inversely proportional to interaction, which is also related to effective mass of tightly bound pairs. At strong interaction limit, the asymptotic behavior of superfluid density is captured by a parabola function of filling factor. In addition, the flat band has significant influences on superfluidity. To be specific, the infinite large density of states of flat band results in a linear interaction-dependence of superfluid order parameters. When interaction is weak, the existences of Dirac points and vanishingly small order parameters cause a logarithmic singularity of superfluid density.
Due to the existence of three order parameters, there are three collective modes, i.e., one is the gapless phonon, which corresponds to total density oscillations; the others are gapped Leggett modes, which are characterized by the relative density oscillations between different sublattices. It is found that the undamped Leggett modes can exist in strong interaction cases. In addition, the collective modes can be reflected by sharp peaks in the two-particle spectral functions. The behaviors of two-particle spectral functions near the phonon modes are also analyzed in detail. In addition, an exact sum-rule of spectral functions is derived. For theoretical aspects, the sumrule can be used to check the various theoretical calculations. In experiments, the sum-rule of the spectral functions may be useful to distinguish between the hole-doped and particle-doped superfluid (or superconductor).  (30) and (32) of spectral functions near the phonon modes with U = 4 and filling factor n = 5 [same as that in panel (b) and (d) of Fig. 5].

(29)
The Bogoliubov transformation of above equation is readily carried out to obtain the eigenenergies E 1,2,3+(−) (k) , where E n− (k) = −E n+ (−k) < 0 and 1, 2, 3 represent three branches of energy. The thermodynamic potential per unit cell in the ground state is where N cell is the number of unit cells. The order parameters (or pairing gaps) and chemical potential can be obtained by solving ∂�/∂� A,B,C = 0 , and n = −∂�/∂µ consistently. In the whole paper, we set t = 1 as the energy unit and lattice spacing a = 1.

Asymptotic behaviors of superfluid density.
For dice lattices, the superfluid order parameters for three sublattices are Using the mean-field decoupling, the BdG Hamiltonian is where the particle and hole part Hamiltonian are and respectively. The above BdG matrix can be diagonalized, e.g., where quasi-particle energy E n+ (k) = E n (k)) > 0 and E n− (k) = −E n+ (−k) < 0 due to particle-hole symmetry of the Bogoliubov-de Gennes Hamiltonian. For three band dice lattices, the indices n = 1, 2, 3 . Eigen-states |n+, k� and |n−, k� ("quasi-particle basis") can be also written as www.nature.com/scientificreports/ where |ps, k� and |ht, k� are eigen-states of particle and hole part Hamiltonian, respectively, i.e, They are also the "energy band basis" for single-particle Hamiltonian. Then W + and W − are transformation matrices between the above two sets of basis. Due to the time reversal symmetry in dice lattices, and |pn, k� = |hn, k� ≡ |n, k�. The thermodynamic potential per unit cell is where N cell is the number of unit cells. The constant Cons = −3µ + The superfluid density can be obtained with phase twist method. Assuming the superfluid order parameters undergo a phase variation, e.g, △(r) → �(r)e 2iq·r , the superfluid density tensor ρ sij can be written as where �(q) is thermodynamical potential. In the presence of phase twist, the BdG Hamiltonian (after applying a gauge transformation) becomes

k)] and the BdG matrix
After diagonalizing the above BdG matrix, we can get the eigen-energies E n+ (k, q) and E n− (k, q) . The thermodynamical potential becomes In addition, if a Hamiltonian H, its eigen-states |n� and eigen-energies E n depend on parameters k and q , they should satisfy (42) (k, q) . |m+� ≡ |m+, k, q� ( |n−� ≡ |n−, k, q� ) and E m+ ≡ E m+ (k, q)[E n− ≡ E n− (k, q)] are eigen-states and eigen-energies of Hamiltonian H BdG (k, q) . The first term of superfluid density depends on the second derivative ∂ 2 E n− ∂k i ∂k j , it would vanish. Because quasi-particle energy E n− (k x , k y ) is a periodic function of k x , the summation over k can be transformed into an integral. An integral of a derivative of periodic function over a period is zero. So the superfluid density is When q → 0, and Using the energy-band basis, the superfluid density (49) H|n� = E n |n�, �m|n� = δ mn , www.nature.com/scientificreports/ When interaction U → ∞ , the chemical potential µ ∝ U and the order parameters A ≃ B ≃ C ∝ U , which corresponds to the uniform pairing case 20 . In such a case, the superfluidity takes place independently in the respective energy bands (no inter-band pairings). The transformation matrix takes the following form In addition, the particle-hole symmetry of BdG matrix implies that the transformation matrices satisfy When U → ∞ , the order parameter and chemical potential where n is filling factor (particle number per unit cell). In addition, the quasi-particle energies satisfy where ε n (k) ≃ −µ as U → ∞ . The superfluid density is reduced to where |n� ≡ |n, k� is energy band basis, and h(k) is single-particle Hamiltonian (Eq. 3).
So the superfluid density where is a rank-two tensor, which only depends on single-particle Hamiltonian. For dice lattices in our work, S ij = 6t 2 δ ij . So the asymptotic formula for superfluid density (particle number per unit cell) is obtained ( U → ∞)

Gaussian fluctuation method.
To investigate the collective modes, it is very convenient to formulate the theory with functional integration. The partition function of grand canonical ensemble is written as functional integral of Grassman fields ψ and ψ where we have the action S given by where the inverse temperature β = 1/T . In this paper, we will only focus on the properties at zero temperature and will take β → ∞ at the end. By introducing the Hubbard-Stratanovich fields � iA(B,C) for sub-lattice A(B,C) as usual, we have (55) W + ns (k) ≃ δ ns In the following, we will use the subscribe i 1 , and ] to denote the positions of sub-lattices A [B/C] unless stated otherwise. With the basis ψ = [ψ i 1 A↑ ,ψ i 2 B↑ ,ψ i 3 C↑ ; ψ i 1 A↓ , ψ i 2 B↓ , ψ i 3 C↓ ] , the inverse Nambu-Gorkov Green function G −1 takes following form where The functional integral is quadratic with respect to ψ fields, we can integrate it out Further assuming iA , iB and iC can be written as � iA (τ ) = � A + δ� iA (τ ) , � iB (τ ) = � B + δ� iB (τ ) and � iC (τ ) = � C + δ� iC (τ ) , where � A(B/C) are not dependent on spatial and time variables. It is convenient to write the G −1 in momentum and Matsubara frequency spaces, i.e., |� iAσ (τ )| 2 U + |� iBσ (τ )| 2 U + |� iCσ (τ )| 2 U −ψG −1 (� iA (τ ), � iB (τ ), � iC (τ ))ψ .