Superfluid density, Josephson relation and pairing fluctuations in a multi-component fermion superfluid

In this work, a Josephson relation is generalized to a multi-component fermion superfluid. Superfluid density is expressed through a two-particle Green function for pairing states. When the system has only one gapless collective excitation mode, the Josephson relation is simplified, which is given in terms of the superfluid order parameters and the trace of two-particle normal Green function. In addition, it is found that the matrix elements of two-particle Green function is directly related to the matrix elements of the pairing fluctuations of superfluid order parameters. Furthermore, in the presence of inversion symmetry, the superfluid density is given in terms of the pairing fluctuation matrix. The results of the superfluid density in Haldane model show that the generalized Josephson relation can be also applied to a multi-band fermion superfluid in lattice.

where ρ s is superfluid density (particle number per unit volume), n 0 is order parameter (condensate density) in liquid Helium-4. G(q, 0) is normal single-particle Green function at zero frequency, m is atomic mass of Helium-4. The above equation indicates that the Green function diverges as wave vector q → 0 . Such a divergence of 1/q 2 in Green function is quite a universal phenomenon which can occur in many systems, e.g., superfluid Helium-4, superconductor and ferromagnet 8 . The above Josephson relation in superfluid system can be viewed as a manifestation of Bogoliubov's " 1/q 2 " theorem in spontaneously symmetry breaking system 7 . It also has close connections with the absence of long ranged order (e.g., condensation) at finite temperature in one and two dimensions 9,10 .
Its possible generalization in two-component fermion superfluid has been firstly investigated by Taylor 11 . Using auxiliary-field approach, Dawson et al. also derived a Josephson relation which holds for both bosonic and fermion superfluids 12 . Here is superfluid order parameter (pairing gap in superconductor) and G II (q, ω) is two-particle Green function for pairing states in conventional two-component fermion superfluid. It is found that the superfluid density is determined by superfluid order parameter and the behaviors of two-particle Green function at long wave length limit. In comparison with bosonic superfluid, the above formula shows that the two-particle Green function replace the corresponding single-particle Green function of bosonic superfluid. In addition, pairing gap plays the roles of order parameter in fermion superfluid.
The superfluid properties in multi-component (or multi-band) fermion system have attracted a great interests [13][14][15][16][17][18][19][20][21] . The exotic pairing mechanism in Fermi gas with SU(N) invariant interaction has been proposed [22][23][24][25][26][27][28][29] , dependent on interaction and chemical potential, which can show coexistence of superfluidity and magnetism 30 . Another interesting example of multi-band fermion system is twisted bilayer graphene 31,32 . It is shown that, there exists superconductivity 33 in this system. Its superfluid weight (which is superfluid density up to a constant) have where m is particle mass, µ is chemical potential, V is volume of system and g is interaction strength parameter between two different (spin) components. In the rest of paper, we set = V = 1 for simplifications. Similarly as bosonic superfluid 39 , here we outline how to get the above relation Eq. (2) in the two-component fermion system. First of all, it is well known that the superfluid order parameter (or pairing gap in superconductor) in two-component fermions is In addition, when the superfluid order parameter undergoes a small phase twist e 2iδθ(r)11,40 , the variation (fluctuation) of superfluid order parameter is where δθ(r) is a real function which denotes the small phase twist.
On the other hand, a phase gradient of superfluid order parameter would induce a superfluid current (superflow) 41 , namely where we define superfluid velocity as the gradient of phase, i.e.,v s ≡ ∇δθ(r)/m and superfluid density ρ s as the coefficient before v s in the current δj(r) . We will see the connection between the above two equations (Eqs. (6), (5)) would result in the Josephson relation.
In order to get the relationship between superfluid density ρ s and order parameter , similarly as bosonic case 9,39,42 , here we need a perturbation which include operator of order parameter � (r) = gψ ↓ (r)ψ ↑ (r) and its where ξ is a small complex number, � q = k gψ ↓q+k ψ ↑−k is fluctuation operator of order parameter in momentum space and we use relation ψ σ (r) = k ψ σ k e ik·r . Here we add an infinitesimal positive number ǫ → 0 + in the above exponential which corresponds to choose a boundary condition that the perturbation is very slowly added to the system 43 .
and ω n0 (ρ q ) n0 = −q · (j q ) n0 , we can obtain where density fluctuation operator ρ q = σ k ψ † σ k ψ σ k+q and we use the fact that �ˆ † q=0 � = * = . So For isotropic system, further assuming q � B ∝ δj and using Eq. (14), so we get Here we further use Eq. where the operator of order parameter � αβ (r) = g αβ ψ α (r)ψ β (r) , ψ α(β) is the field operator for α(β)-th component, and g αβ is interaction parameter between α-th and β-th components. The above equation indicates that the number of superfluid order parameter may be an arbitrary integer m ≥ 1 ( m ∈ Z ) in a multi-component fermion superfluid. In such a case, we need a general perturbation Hamiltonian where we relabel the operator of order parameter with � i (i = 1, 2, ..., m) and introduce column vectors .., ξ m } t and {...} t denotes matrix transpose and dot · is the multiplication of matrices.
In the above subsection, the definition of superfluid velocity involves particle mass m. However, the superfluid system we considered here can be a general many-body system with complex energy spectra. For an arbitrarily general multi-component (or multi-band) system, it may be difficult to introduce notion of mass. Furthermore, the superfluid velocity may be unable to be defined unambiguously. In the following manuscript, the superfluid density is simply defined as the coefficient before the gradient of phase in current, i.e., For a general multi-component superfluid system, the phase variation of order parameter δθ is well-defined. Consequently, the physical meaning of superfluid density ρ s is also clear. In fact, the above definition of superfluid density is also consistent with the phase twist method in lattice system 46 .
In addition, even though Eq. (37) is obtained in a translational invariant system (the momentum is a good quantum number), the above formula can be applied equally to lattice system as long as the superfluid state has lattice translation symmetry (see also "Summary" section). This is because, we know for lattice system, the engenstates can be classified by lattice momentum. In the above derivation, on the one hand, we need to identify the (canonical) momentum P = k,σ kψ † σ k ψ σ k with the lattice momentum. On the other hand, the superfluid density may show some anisotropy in lattice system. We can use a similar method as done in bosonic system 39 to deal with the anisotropy in lattice. In such a case, the superfluid density is usually a rank-two tensor and depends on direction of q.
In isotropic case, the superfluid density is a scalar, i.e., ρ s = diag{ρ s , ρ s , ρ s } which does not depend on direction of q . For generally anisotropic system (for example, Lieb lattice system 36 , the spin-orbital coupled cold atoms 47 , etc), the induced current δj can be expressed in terms of vector q and a rank-two tensor m (generally δj is not parallel to q any more), namely where indices i(j) = x, y, z in three dimensional space. In order to get Josephson relation in anisotropic system, we introduce the superfluid density ρ s (q) along q direction, i.e., q · δj(r) ≡ ρ s (q)(q · ∇δθ) , where q ≡ q/q is unit vector along q direction. From Eqs. (31) and (38), we get which is exactly similar to Eq. (34) of isotropic case. Taking qx(r) = −iq · ∇x , qx * (r) = iq · ∇x * and Eq. (33) into account, the following discussions for anisotropic case are same with that in the isotropic case. So the superfluid density along direction q is which usually depends on the direction of q . Furthermore, from the superfluid density along an arbitrary direction q , i.e., it is not difficult to construct a rank-two superfluid density tensor ρ s;ij .
Discussions. In the above derivations, we get the superfluid density in terms of two-particle Green function [the Josephson relation Eq. (37)]. In this section, we will show that for some cases, the above equation can be further simplified. To be specific, if the system has only one gapless collective mode, the superfluid density can be give by the trace of two-particle Green function. When the inversion symmetry is present, the superfluid density can be expressed in terms of the fluctuation matrix of superfluid order parameters.
Only one gapless collective mode. When a system has unique gapless excitation near ground state, e.g., phonon, the Josephson relation can be generalized to a multi-component system with a phase operator method as done in bosonic case 39 . Here we know near the ground state, the phonon corresponds to total density oscillation. Due to the presence of superfluid order parameter, the density oscillation would couple phase oscillation of order parameter 48 . Furthermore, all the superfluid order parameters should share a common phase variation, i.e., δθ i (r) = δθ(r) . On the other hand, near the ground state, similarly as Eq. (5), the fluctuation operator of superfluid order parameters may be expressed in terms of phase operator θ 40 (36) δj(r) = − 4∇δθ(r) q 2 (� t * , � t ) .

From the above equation, we get the fluctuation operators of order parameters in momentum space
where we use θ † q =θ −q for real phase field θ(r) ( q = 0 ) . From definitions of the G II and F II in Eq. (28), we get ] > 0 is a real number. When the number of superfluid order parameters is one and is real, the relation F II (q, 0) = −G II (q, 0) holds ( q → 0 ) in conventional two-component fermion superfluid.
From Eqs. (33) and (34), we get where we set � * i ξ i ≡ α i e iφ i with amplitude α i , phase φ i . So the superfluid density is On the other hand, we know So finally we get the Josephson relation where G II (q, 0) is two-particle normal Green function (matrix) at zero-frequency. It is found that the above Eq. (44) can be applied to the case of three superfluid parameters in dice lattice, where there exist one gapless phonon, and two gapped Leggett modes arising from the relative phase oscillations between different superfluid order parameters 21 . When the number of order parameters m = 1 , the above equation is reduced to the Josephson relation for conventional two component fermion superfluid 11,12 .
Superfluid density and pairing fluctuations. In this section, we would give a connection between the two-particle Green function and the pairing fluctuation matrix based on BCS mean field theory and Gaussian fluctuation approximation 11 . Furthermore, if the system has spatial inversion symmetry, the coefficient matrix G II is directly proportional to the inverse of pairing fluctuation matrix.
Fist of all, we discuss the results for conventional two-component fermion superfluid. In the following, we assume the pairing gap can be decomposed as the mean field value and small fluctuation δ� ( �(r) = � + δ� ). Expanding action S to second order of δ� , the partition function [49][50][51] where S 0 is the mean-field contribution and Gaussian fluctuation part (40) δ� i (r) ≃ 2i� i δθ (r), where pairing fluctuation field η † q = [� * q , � −q ] and q = (q, iω n ) . ω n = 2nπ/β ( n ∈ Z ) is Matsubara frequency, β = 1/T is inverse temperature. The fluctuation matrix M is a 2 × 2 matrix where G ij (k) is matrix element of Nambu-Gorkov Green function. The collective modes are given by zeros of determinant Det|M(q, iω n → ω + i0 + )| = 0 . As q → 0 , the collective mode is the Anderson-Bogoliubov phonon, which characterizes the density oscillations of superfluid. With the action δS (Gaussian weight), the correlation function (average values of quadratic terms) can be calculated 52 , i.e., On the other hand, we know that the above correlation function has one extra minus sign in comparison with Green function. So the two-particle Green function can be obtained by In addition, if the system has inversion symmetry, i.e., Furthermore, according to Eq. (28), when ω ± iǫ = 0 , the Green function satisfy then the coefficient matrix of Green function G II can be expressed in terms of the inverse of M, i.e.,

So the superfluid density
Next, for the case of several superfluid order parameters, the proof is similar. Assuming the number of superfluid order parameters is an arbitrary integer m, then Gaussian fluctuation part of action (50) G II (−q, iω n ) = G II (q, iω n ), F II (−q, iω n ) = F II (q, ω n ).
= lim q→0 4 q 2 (� t * , � t ). www.nature.com/scientificreports/ w h e r e M ( q ) i s a 2m × 2m f l u c t u at i o n m at r i x , a n d p a i r i n g f l u c t u at i o n f i e l d η † q = [� * 1q , � * 2q , ..., � * mq , � 1,−q , � 2,−q , ..., � m,−q ] . With the action δS (Gaussian weight), the correlation functions (and the matrix elements of Green function) can be calculated, i.e., Furthermore, in the presence of inverse symmetry, the coefficient matrix can be obtained in terms of the fluctuation matrix, i.e., G II = −M −1 . So the superfluid density is given by During the derivations of Eqs. (37), (44) and (56), we can see that the superfluid density is mainly related to the low energy collective excitations (e.g., phonon for neutral fermion superfluid) and the behaviors of Green functions at long wave length limit. In addition, it is found that the superfluid density in a multi-band fermion superfluid can be divided into two parts, one is the conventional part, which arises from the diagonal matrix elements of current operator, and the other one is the off-diagonal terms (the so called geometric part 36,46 ), which connects with the geometric metric tensor of Bloch states 53 . The geometric part plays an important role in the understanding of superfluidity of flat band, where the conventional part is negligible 36 . In addition, it is expected that some other physical quantities, for example, sound velocity 20 should also have similar geometric part in a multi-component (or multi-band) superfluid. The Josephson relations [Eqs. (37), (44) and (56)] give the total superfluid density, which not only includes conventional part, but also the geometric part of a multiband fermion superfluid.
The non-vanishing superfluid density in superconductor can result in the perfect diamagnetism (the Meissner effect) 54 . The penetration depth of a magnetic field ( ) is related to the superfluid density through relation ρ s ∝ 1/ 2 . So the superfluid density can be experimentally obtained by measuring the penetration depth of magnetic field in superconductor. On the hand, for two-dimensional high-T c superconductors, the superfluid density at zero-temperature is related to the superfluid transition temperature by T c ∝ ρ s (T = 0) (the so called Uemura relation 55 . So the superfluid density can be also evaluated by measuring the superconductor transition temperature T c . For neutral superfluid system, the superfluid density would result in a reduction of moment of inertia of system (from its classical value). So the superfluid density can be obtained through measuring the moment of inertia of superfluid system 56,57 . In order to measure the moment of inertia of atomic gas, an optical method is proposed through imparting non-zero angular momentum into cold atom gas 58 .
An example: superfluid density in Haldane-Hubbard model. As an application of the Josephon relations, e.g., Eqs. (44) and (56), the superfluid density is calculated for Haldane-Hubbard model in twocomponent Fermi gas with on-site attractive interaction −U ( U > 0) 59

. The Hamiltonian for Haldane-Hubbard model is
where µ is chemical potential, and n iσ = c † iσ c iσ is particle number operator. ǫ i = 1(−1) for sublattice A (B) and M is energy offset between sublattice A and B. t ij is hopping amplitude between lattice sites i and j, which is t for nearest neighbor sites, t ′ e −iφ (t ′ e iφ ) for clockwise (anticlockwise) hopping between next-nearest neighbor sites 59 . The distance between nearest neighbor sites is a.
In the following, we focus on the case of φ = π/2 and M = 0 , where the inversion symmetry is not broken. Due to the presence of two sublattices, Haldane model is a two-band fermion system. In addition, there exist two superfluid order parameters for two sublattices, i.g., where i is unit cell index, ↓ (↑) are two spin component indices. Furthermore, when M = 0 , the two order parameters are equal, e.g., i,A = i,B = . When the filling factor n = 2 (2 particles per unit cell in half-filling), for weakly interacting case, the system is a band insulator, and the superfluid order parameter vanishes ( = 0 ). Only when the interaction is strong enough, the system enters a superfluid phase ( = 0 ). Away from the half-filling, the system is usually in a superfluid phase 59 . In addition, it is found that, there exists only one gapless excitation near q = 0 , which corresponds to total density oscillations. The superfluid density can be also calculated with current-current correlation functions 9,34,47 or phase twist method 46,60 . Assuming the superfluid order parameters undergo a phase variation, e.g, △ iA(B) → � iA(B) e 2iq·r iA(B) , the superfluid density (particle number per unit cell) tensor ρ sij can be written as (56) ρ s = lim q→0 4 q 2 (� t * , � t ). where is thermodynamical potential (per unit cell) in grand canonical ensemble. Figure 1 shows the evolutions of superfluid density as the interaction increases. The superfluid density is obtained with three different formulas, i.e. Eqs. (58), (56) and (44). First of all, it is found that the superfluid density is isotropic and behaves as a scalar in two-dimensional space, i.e., ρ sij = diag(ρ s , ρ s ) , which is a consequence of C 3 rotational symmetry of honeycomb lattice 46 . Secondly, the results from the Josephson relations [Eqs. (56), (44)] are consistent with the results from phase twist method [Eq. (58)]. In addition, When filling factor n = 2.1 , the superfluid density increases as the interaction gets strong. However when the filling factor n = 2.8 , the superfluid density gets smaller and smaller when the interaction increases. Such an interesting feature is also reflected in Fig. 2, that when the filling factor is near half-filling ( n = 2 ), the superfluid density increases as the interaction gets strong. However, when the filling is far away from half-filling, the situation is reversed (see Fig. 2). Figure 2 shows the superfluid density as functions of filling factor ( n = 0 → 4 ). From the Fig. 2, we can see that the superfluid density is symmetric with respect to half-filling ( n = 2 ) due to particle-hole symmetry 59 . For weak interaction cases ( U = 2t and U = 3t ) and half filling, the system is a insulator (see Ref. 59 ), the superfluid density vanishes (see Fig. 2). For fully occupied case ( n = 4 particles per unit cell), the system is equivalent to the completely empty case ( n = 0 ) due to the particle-hole symmetry, the system is also a insulator. Consequently, superfluid density is also zero. When the filling factor falls into the middle of upper band ( n ≈ 3 ), the superfluid density reaches its maximum. The appearance of double dome structure for weak interactions in Fig. 2 is in qualitative agreement with the results obtained through dynamical mean-field theory (DMFT) in Ref. 46 .

Summary
In conclusion, we investigate the Josephson relation for a general multi-component fermion superfluid. It is found that the superfluid density is given in terms of two-particle Green functions. When the superfluid has only one gapless collective excitation, the Josephson relation can be simplified, which is given in terms of superfluid order parameters and trace of Green function. Within BCS mean field theory and Gaussian fluctuation approximation, the matrix elements of Green function can be given in terms of pairing fluctuation matrix elements. Furthermore, in the presence of inversion symmetry, it is shown that the two-particle Green function is directly proportional the inverse of pairing fluctuation matrix. The formulas for superfluid density are quite universal for generic multi-component fermion superfluids, which can be also applied to a multi-band superfluid with complex energy spectra in lattice.
Josephson relation for multi-component fermion superfluid provides a general method for calculations on superfluid densities in terms of two-particle Green functions and fluctuation matrix. Our work would be useful for investigations on the superfluid properties of multi-component (or multi-band) superfluid system with complex pairing structures.