Conserving Approximation of Pairing theories in Fermionic superfluid phase
aa r X i v : . [ c ond - m a t . qu a n t - g a s ] N ov Conserving Approximation of Pairing theories in Fermionic superfluid phase
Yan He and Hao Guo College of Physical Science and Technology, Sichuan University, Chengdu, Sichuan 610064, China and Department of Physics, Southeast University, Nanjing 211189, China (Dated: April 19, 2018)Respecting the conservation laws of momentum and energy in a many body theory is very impor-tant for understanding the transport phenomena. The previous conserving approximation requiresthat the self-energy of a single particle can be written as a functional derivative of a full dressedGreen’s function. This condition can not be satisfied in the G G t-matrix or pair fluctuation theorywhich emphasizes the fermion pairing with a stronger than the Bardeen-Cooper-Schrieffer (BCS)attraction. In the previous work ? , we have shown that when the temperature is above the su-perfluid transition temperature T c , the G G t-matrix theory can be put into a form that satisfiesthe stress tensor Ward identity (WI) or local form of conservation laws by introducing a new typeof vertex correction. In this paper, we will extend the above conservation approximation to thesuperfluid phase in the BCS mean field level. To establish the stress tensor WI, we have to includethe fluctuation of the order parameter or the contribution from the Goldstone mode. The resultwill be useful for understanding the transport properties such as the behavior of the viscosity ofFermionic gases in the superfluid phases. PACS numbers:
I. INTRODUCTION
For strongly correlated systems, such as ultra-cold Fermi gases in the unitary limit with divergent scattering length,perturbation calculations are not reliable because of the lack of small parameters. To capture the strong fluctuations,various approximating methods have been invented, one important example is the t-matrix theory. The t-matrixtheory emphasizes the pairing effect between the fermions due to the stronger-than-BCS attractions. It provides anatural explanation of the preformed pairs before condensation and the formation of the pseudogap. The t-matrixtheory has many applications to the ultra-cold Fermi gas because of the experimental tunable interactions.In recent years a lot of experimental ? and theoretical works ? ? ? ? in ultra-cold Fermi gases focus on the transportphenomena, such as the behavior of the shear viscosity. In order to give reliable calculations of it, the chosen manybody theory should respect the momentum and energy conservation laws. This is the reason that Baym and Kadanoffproposed the conserving approximation long time ago ? ? . The conserving approximation actually sets up a morestringent consistency requirement to the many body theory. Based on the similar consideration, Haussmann andZwerger studied a t-matrix theory with two full dressed propagator in the ladder which is also known as GG theory ? .There are also two other important types of t-matrix theories. The simplest version of t-matrix theory known as G G theory with two bare fermion propagator in the ladder was pioneered by Nozieres and Schmitt-Rink ? ? . Theother one is the G G theory first introduced by Kadanoff and Martin ? , which is designed to be compatible with theBCS-Leggett generalization ? ? . The conserving approximation condition is usually considered not satisfying by thesetwo types of t-matrix theories. However, in our previous work ? we have shown that certain vertex corrections can beintroduced such that the Ward identities of stress tensor or energy-momentum tensor are satisfied both in the G G and G G t-matrix theories above T c , which also means that the momentum and energy are conserved locally. Thisresult opens a way to improve these two t-matrix theories. In this paper, we will show how the local conservationlaws are satisfied in the superfluid phase.Since the way that the conservation law is satisfied in our theory is quite different from that in the traditionalconserving approximation, it is worth to make a detailed comparison of both approaches. In general, the local formsof conservation laws are expressed by the operator equations such as ∂ µ j µ = 0 and ∂ µ T µν = 0. Inserting theseoperators into a n -particle Green’s function, one can find that the divergence of n -particle Green’s function equals tosome contact terms. These are the most original forms of the Ward identities. They must be respected in the exacttheory, however they are usually violated in various levels in the approximated theories.The key point of traditional conserving approximation is that the self-energy of the one-particle Green’s functionis a functional derivative Σ = δ Φ /δG . Here Φ is a functional of the dressed Green’s function G , which also appearsas an interacting part of the free energy in Luttinger-Ward formalism ? . In the field theory language, Φ can alsobe represented as the 2-particle-irreducible skeleton vacuum diagrams. While in real practical calculations, one hasto choose a particular class of diagrams. Due to this truncation, not all the WIs are satisfied in the conservingapproximation. What has been proved is that the one-particle Green’s functions satisfy the conservation law and thetwo-particle Green’s functions satisfy the thermodynamical consistency. It must be noted that the conservation ofthe momentum and energy are proved for the whole system rather than a local form. Moreover, certain importantsymmetries such as the crossing symmetry determined by Pauli principle are also violated. In another approachproposed by de Dominicis and Martin ? , the full dressed two-particle scattering vertex has been introduced as afundamental quantitywhich is determined by the parquet equations. In this scheme, the crossing symmetry is nowrespected but the conservation law is not always satisfied. The FLEX approximation proposed by Bickers ? is roughlyequivalent to the combination of the two approaches mentioned above, however there is still some inconsistency suchthat the vertex derived from the functional derivative of self-energy is not the same as the one used to calculate theself-energy. In summary, the conserving approximation does not mean that all WIs are automatically satisfied. Inorder to respect crossing symmetry, one has to treat the full vertex on the same footing as the self-energy.On the other hand, since the G G or G G self-energy can not be written as a functional derivative, it is believedthat in these theories the conservation laws are not respected. In our theory, we focus on the satisfaction of stresstensor WI or local conservation instead of pursuing the integral form of conservation law. Moreover, the canonicalstress tensor is treated as an external classical field, and we can introduce an external vertex associating with it.The interaction term in the stress tensor serves as a new bare stress tensor vertex. Although the stress tensor WIhas a very complicated momentum dependence, we can rearrange the momentum dependence by introducing a newvertex correction corresponding to the self-energy and show that the WI is satisfied by the full vertex. This alsoimplies that the WIs are satisfied for the one and two-particle Green’s functions. When we apply the same schemeto the superfluid phase, we have encountered huge difficulties because the effects of Goldstone modes and the pairfluctuation are mixed with each other. In this paper, we will stick to the BCS mean field level as an easy startingpoint, the work that the pair fluctuation is included will be built on the top of the BCS mean field theory in thefuture. As we will see later that, even in this level the problem of satisfying stress tensor WI is non-trivial. Althoughthe translational symmetry is not broken, the Goldstone modes of the broken U (1) symmetry plays a very importantrole in constructing the stress tensor correlation function which satisfies the conservation laws. II. STRESS TENSOR WARD IDENTITY (WI) IN NORMAL PHASE
Since the WI associated with the stress tensor or momentum current is not very familiar to the condensed matterphysics community, we first establish these identities for the general scalar fields in this section, and then apply thesegeneral formula to the simplest case, the non-interacting Fermi gas. WIs generally represent the continuity equationsexpressed in terms of Green’s functions or correlation functions. They are obtained by sandwiching the operatorcontinuity equations in various time ordered operator products. These WIs must be satisfied in any theory no matterit is exact or approximate. As a beginning, we start with a simple example about the particle number or the U (1)current conservation which is given by ∂ρ∂t + ∇ · J = 0. The Ward identity of this example is given by q µ Γ µ ( K + Q, K ) = G − ( K + Q ) − G − ( K ) , (1)which is familiar in QED. Here K ≡ k µ = ( ω, k ) and Q ≡ q µ = ( q , q ) are the 4-momenta, G is the dressed Green’sfunction and Γ is the dressed vertex defined by the equation h ψ ( x ) J µ ( z ) ψ † ( y ) i ≡ Z G ( x, x ′ )Γ µ . ( x ′ , y ′ , z ) G ( y ′ , y ) d x ′ d y ′ , (2)where x ≡ x µ = ( t, x ) is the 4-coordinate (We work in the unites where ~ = c = 1 throughout this paper). Moregenerally, inserting the current conservation continuity equation into a n -particle Green’s function of some scalarfields, one can find that the divergence of n -particle Green’s function is not simply zero. ∂ µ h J µ ( x ) φ ( x ) · · · φ ( x n ) i = − X i δ ( x − x i ) h φ ( x ) · · · φ ( x n ) i . (3)The right hand side is often referred as the “contact terms” which comes from the fact that the time derivative hitson the time order products. We will see these terms also appear in other cases.The Ward identity associated with the momentum conservation is more complicated but less familiar, and it is whywe focus on it in this paper. According to the Noether’s theorem, the canonical stress tensor is given by T µν = X a ∂ L ∂ ( ∂ µ φ a ) ∂ ν φ a − g µν L , (4)where L is the Lagrangian density, g µν = (1 , − , − , −
1) is the metric tensor and the index a labels different speciesof fields.The stress tensor satisfies the conservation law ∂ µ T µν = 0. The general WI associated with the correlation functionof the stress tensor and other scalar fields in the coordinate space is given by ? ∂ µ h T µν ( x ) φ ( x ) · · · φ ( x n ) i = − X i δ ( x − x i ) ∂∂x νi h φ ( x ) · · · φ ( x n ) i , (5)where the terms on the right-hand-side are also the contact terms as in the current conservation case. The derivativein the contact terms reflect that the momentum is associate with spatial translation. Applying the above expressionto the 3-point correlation function, we get ∂ µ h T µν ( x ) ψ † ( y ) ψ ( z ) i = − δ ( x − y ) ∂∂y ν h ψ † ( y ) ψ ( z ) i − δ ( x − z ) ∂∂z ν h ψ † ( y ) ψ ( z ) i . (6)Transforming to the momentum space, we find q µ Γ µν ( K + Q, K ) = k ν G − ( K + Q ) − ( k + q ) ν G − ( K ) . (7)This is the general form of the WI associated with the stress tensor for both non-interacting and interacting systems,which will be the focus of this paper. Comparing it to the WI associated with the U (1) current, i.e. Eq. (1), we seethat this WI has a subtle momentum dependence, which makes the establishment of it for the interacting Fermi gasesmuch more difficult than that associated with the U (1) current. To get a warm-up experience of the stress tensor WI,we begin with the non-interacting Fermi gas with the following Lagrangian density L = i ψ † σ ∂ t ψ σ − i ∂ t ψ † σ ψ σ − m ∂ i ψ † σ ∂ i ψ σ + µψ † σ ψ σ . (8)Here ψ σ is the Fermionic field with σ = ↑ , ↓ is the spin index, and we have taken the Einstein summation conventionand a symmetric form for the time derivative. The equation of motion is then given by ∂ L ∂ψ σ − ∂ t ∂ L ∂ ( ∂ t ψ σ ) − ∂ i ∂ L ∂ ( ∂ i ψ σ ) = − i∂ t ψ † σ + 12 m ∂ i ψ † σ + µψ † σ = 0 . (9)which is just the Schrodinger equation for a non-interacting Fermionic field.From Eq. (4), the components of the canonical stress tensor involving momentum density and momentum currentare given by T j = − ( i ψ † σ ∂ j ψ σ − i ∂ j ψ † σ ψ σ ) ,T ij = 12 m ( ∂ i ψ † σ ∂ j ψ σ + ∂ j ψ † σ ∂ i ψ σ ) + δ ij L , (10)which satisfy momentum current conservation ∂ t T j + ∂ i T ij = 0. One can see that in a non-relativistic theory, themomentum density is essentially the same as the U (1) current J j = T j /m .In general, the stress tensor is not uniquely defined. Different forms of the stress tensor lead to different forms ofWI. The component T ij contains a time derivative which makes the frequency summation quite complicated. We canmake use of the equation of motion to get rid of the time derivative and get a more convenient expression of T ij .Hence one finds T j = − ( i ψ † σ ∂ j ψ σ − i ∂ j ψ † σ ψ σ ) ,T ij = 12 m ( ∂ i ψ † σ ∂ j ψ σ + ∂ j ψ † σ ∂ i ψ σ ) − δ ij ∂ i ( ψ † σ ψ σ )4 m . (11)The corresponding (bare) vertices are given by γ j ( K + Q, K ) = k j + q j ,γ ij ( K + Q, K ) = ( k + q ) i k j + ( k + q ) j k i m + δ ij q m . (12)We will stick to this simplified version of vertex in the rest of this paper. We can also verify the simplified WIassociated with these vertices as q µ γ µj ( K + Q, K ) = ( k j + q j G − ( K + Q ) − G − ( K )] . (13)It is simplified a little bit comparing to the previous version of WI in Eq.(7). We are interested in the physical propertyof the 2-point stress tensor response functions which is useful to the understanding of transport property such as theviscosity. The stress tensor-stress tensor correlation function for any Fermi gas is given by Q µj,ab ( x − y ) = − iθ ( x − y ) h [ T µj ( x ) , T ab ( y )] i . (14)By implementing the conservation law of the stress tensor, the divergence of the stress tensor-stress tensor correlationfunction in the coordinate space is evaluated as ∂ µ Q µj,ab ( x − y ) = − iδ ( x − y ) h [ T j ( x ) , T ab ( y )] i . (15)After performing the Fourier transformation, we find the corresponding expression in the momentum space of theabove equation q µ Q µj,ab ( Q ) = h [ T j ( q , t ) , T ab ( − q , t )] i = X p , k (cid:16) p + q (cid:17) j γ ab ( K + Q, K ) h [ ψ † σ p ψ σ p + q , ψ † σ k + q ψ σ k ] i = X k (cid:16) k + q (cid:17) j γ ab ( K + Q, K )( n k − n k + q ) . (16)In the specific case of a non-interacting Fermi gas, the divergence of the stress tensor-stress tensor correlation canalso be obtained by the diagrammatic method and making use of the vertex WI as q µ Q µj,ab ( Q )= X K q µ λ µj ( K, K + Q ) G ( K + Q ) γ ab ( K + Q, K ) G ( K )= X K (cid:16) k + q (cid:17) j h G ( K ) − G ( K + Q ) i γ ab ( K + Q, K ) , (17)which is consistent with the above generic result. III. STRESS TENSOR WI IN SUPERFLUID PHASE
Now we move on to the main focus of this paper, we are going to generalize the above results to the BCS superfluid,i.e. to construct the stress tensor-stress tensor correlation function or the full stress tensor vertex which satisfies WI.From the expression of the bare vertex γ ij , one might guess that a bubble diagram with two vertex insertions maygive rise to the correct stress tensor correlation function. However, one can easily verify that this correlation functionviolates the WI in the superfluid phase. The reason is similar to what happens to the current-current correlationfunction. In the superfluid pase, the broken U (1) symmetry generates the Goldstone mode, whose contribution cannotbe ignored in general gauge choices. Although the translational invariance is not broken, the Goldstone mode stillmakes important contributions to the stress tensor response theory, as we will see in the following discussions.For convenience, we use a more compact expression in the Nambu space to discuss our theory. We introduce thetwo-component Nambu fermion Ψ † k = ( ψ ↑ k , ψ †↓− k ), then the Hamiltonian is H = X k Ψ † k ( ξ k σ − σ ∆)Ψ k , (18)where ξ k = k m − µ is the free fermion dispersion, and ∆ is the order parameter given by ∆( x ) = − g h ψ ↑ ( x ) ψ ↓ ( x ) i with g being the coupling constant. It is easy to see that the inverse bare and full matrix Green’s functions areˆ G − ( P ) = iω n − ξ p σ , ˆ G − ( P ) = iω n − ξ p σ + ∆ σ . (19)By direct evaluation, we get the explicit form of the matrix Green’s function in the Nambu spaceˆ G ( K ) = (cid:18) G ( K ) F ( K ) F ( K ) − G ( − K ) (cid:19) . (20)Here G ( K ) and F ( K ) are the well-known BCS Green’s function and anomalous Green’s function G ( K ) = u k iω n − E k + v k iω n + E k , (21) F ( K ) = − u k v k (cid:16) iω n − E k − iω n + E k (cid:17) , (22)where the quasi-particle dispersion E k = p ξ k + ∆ and coherent factor u k , u k = (1 ± ξ k /E k ). In the Nambu space,the gap equation and number equation can be written as∆ = g X K Tr[ σ ˆ G ( K )] , n = X K Tr[ σ ˆ G ( K )] . (23)Recall that in the Nambu space the bare external electromagnetic vertex in the Nambu space is given by ˆ γ µ ( K + Q, K ) = ( σ , k + q m ). Similarly, one can deduce that the bare vertex in the stress tensor linear response theory in theNambu space becomes ˆ γ j ( K + Q, K ) = ( k j + q j , (24)ˆ γ ij ( K + Q, K ) = h ( k + q ) i k j + ( k + q ) j k i m + δ ij q m i σ . (25)It is straightforward to verify that this bare vertex satisfies the bare WI as follows q µ ˆ γ µj ( K + Q, K ) = ( k + q j h σ ˆ G − ( K + Q ) − ˆ G − ( K ) σ i . (26)To obtain the full stress tensor-stress tensor correlation function which satisfies the WI, we must include the fluctuationof ∆ on the same footing as that of the external metric field which couples to the vertex γ µj . Hence we introduce thegeneralized external disturbing field ˆ Φ = (cid:0) ∆ , ∆ , g νl (cid:1) T (27)and the corresponding vertex ˆ Λ ( K + Q, K ) = (cid:16) σ , σ , ˆ γ νl ( K + Q, K ) (cid:17) T . (28)Then the perturbing Hamiltonian H ′ can be cast in the scalar product as H ′ = X pq Ψ † p + q ˆ Φ · ˆ Λ ( K + Q, K )Ψ p = X pq Ψ † p + q ˆ Φ a ˆ Λ a ( K + Q, K )Ψ p . (29)Introducing the imaginary time formalism, the generalized perturbation η of the stress tensor current due to theperturbative Hamiltonian H ′ is given by η ( τ, q ) = X p h Ψ † p ( τ ) ˆ Λ ( P + Q, P )Ψ p + q ( τ ) i , with ˆ Λ ( K + Q, K ) = (cid:16) σ ( k + q j , σ ( k + q j , ˆ γ µj ( K + Q, K ) (cid:17) T . (30)Here η µ = h T µj i denotes the perturbed expectation value of the stress tensor and η , denotes the perturbed real andimaginary parts of the order parameter multiplied by external momentum. The generalized linear response theory inthe momentum space can be written as η a ( Q ) = X b Q ab ( Q )Φ b ( Q ) . (31)More explicitly in the matrix form, we have η j η j η µj = Q j Q j Q j,νl Q j Q j Q j,νl Q µj Q µj Q µj,νl ∆ ∆ g νl (32)where the response functions are defined as follows Q ab (Ω , q ) = Tr T X iω n X p (cid:0) ˆΛ a ( P + Q, P ) ˆ G ( P + Q )ˆΛ b ( P, P + Q ) ˆ G ( P ) (cid:1) . (33)Detailed expression of these response functions can be found in the appendix. The component Q µj,νl is the gener-alization of the “bare” counterpart given by Eq. (14). However, it does not satisfy the conservation law since thecontribution of the Goldstone mode is not included. Our goal is to find out the “full” response kernel which cancorrectly reflects the spacetime symmetry of the theory, i.e., satisfies the WI. It can be constructed from the gener-alized response function Q ab introduced above. Before we do that, we first prove three identities, which are in factthe WIs associated with the response functions since these identities impose the spacetime symmetry on the quantumcorrelation functions. q µ Q µj = − i ∆ Q j ,q µ Q µj = − i ∆ Q j ,q µ Q µj,νl = − i ∆ Q j,νl + h [ T j , T νl ] i . (34)To prove them, we need to verify that the matrix BCS Green’s function satisfies( p + q j h σ ˆ G − ( P + Q ) − ˆ G − ( P ) σ i = q µ ˆ γ µj ( P + Q, P ) + 2 i ∆ σ ( p + q j . (35)This is in fact the WI associated with the Green’s function and the stress tensor vertex in the Nambu space. Theproof of it is quite straightforward.For the first identity of Eqs.(34), we have q µ Q µj + 2 i ∆ Q j = Tr X P h(cid:0) q µ ˆ γ µj ( P + Q, P ) + 2 i ∆ σ ( p + q j (cid:1) ˆ G ( P + Q ) σ ˆ G ( P ) i = Tr X P h ( p + q j (cid:0) σ ˆ G − ( P + Q ) − ˆ G − ( P ) σ (cid:1) ˆ G ( P + Q ) σ ˆ G ( P ) i = Tr X P ( p + q j (cid:2) iσ ˆ G ( P ) + ˆ G ( P + Q ) iσ (cid:3) = 0 (36)For the second identity, we have q µ Q µj + 2 i ∆ Q j = Tr X P h(cid:0) q µ ˆ γ µj ( P + Q, P ) + 2 i ∆ σ ( p + q j (cid:1) ˆ G ( P + Q ) σ ˆ G ( P ) i = Tr X P h ( p + q j (cid:0) σ ˆ G − ( P + Q ) − ˆ G − ( P ) σ (cid:1) ˆ G ( P + Q ) σ ˆ G ( P ) i = − Tr X P ( p + q j (cid:2) iσ ˆ G ( P ) + ˆ G ( P + Q ) iσ (cid:3) = 0 (37)For the last one, we have q µ Q µj,νl + 2 i ∆ Q j,νl = Tr X P h ( p + q j (cid:0) σ ˆ G − ( P + Q ) − ˆ G − ( P ) σ (cid:1) ˆ G ( P + Q )ˆ γ νl ( P, P + Q ) ˆ G ( P ) i = Tr X P ( p + q j h σ ˆ γ νl ( P + Q, P ) ˆ G ( P ) − ˆ G ( P + Q )ˆ γ νl ( P, P + Q ) σ i = X P ( p + q j Tr (cid:0) [ ˆ G ( P ) σ − σ ˆ G ( P + Q )]ˆ γ νl ( P + Q, P ) (cid:1) = h [ T j ( q , t ) , T νl ( − q , t )] i (38)Now we are ready to construct the full stress tensor response function which satisfies the WI. The key point is thatthe order parameter is not an arbitrary external field but self-consistently determined. The fluctuation of the phaseof the order parameter is the Goldstone mode that plays an important role in the gauge invariant current responsetheory. It also gives important contribution to the stress tensor response theory. Therefore, the fluctuation of theorder parameter must be included in our linear response theory on the same footing as the fluctuation of the externalmetric field. After imposing the self-consistent condition, the perturbed order parameter can be solved and correctlyinclude the contribution of the Goldstone mode.The gap equation gives the self-consistent condition η , = 0. Applying this relation to Eq.(32), we find theperturbed order parameter∆ = − Q j,νl Q j − Q j,νl Q j Q j Q j − Q j Q j g νl , ∆ = − Q j,νl Q j − Q j,νl Q j Q j Q j − Q j Q j g νl . (39)After inserting these results into T µj = Q µj ∆ + Q µj ∆ + Q µj,νl g νl , (40)we get the usual Kubo expression T µi = K µi,νl g νl where the full stress tensor response kernel K µj,νl is given by K µj,νl = Q µj,νl − ( Q j,νl Q j − Q j,νl Q j ) Q µj + ( Q j,νl Q j − Q j,νl Q j ) Q µj Q j Q j − Q j Q j . (41)where the poles of the second term give the excitations of the collective mode. Hence it only appears in the superfluidphase or more generally the broken-symmetry phase. Now we show that the full stress tensor response kernel doessatisfy the conservation law. With the help of the WIs (34), we find that q µ K µj,νl = − i ∆ Q j,νl + h [ T j , T νl ] i + 2 i ∆ ( Q j,νl Q j − Q j,νl Q j ) Q µj + ( Q j,νl Q j − Q j,νl Q j ) Q µj Q j Q j − Q j Q j = h [ T j , T νl ] i . (42)which is also the WI associated with the two-point stress tensor correlation function.From the expression of K µj,νl g νl , we can deduce the full stress vertex asˆΓ µj ( P + Q, P ) = ˆ γ µj ( P + Q, P ) − ( p + q j σ Π µ ( Q ) − ( p + q j σ Π µ ( Q ) (43)where Π , are defined as Π µ ( Q ) = Q j Q µj − Q j Q µj Q j Q j − Q j Q j , Π µ ( Q ) = Q j Q µj − Q j Q µj Q j Q j − Q j Q j . (44)Then the full response kernel can be expressed as K µj,νl ( Q ) = X P h ˆΓ µj ( P + Q, P ) ˆ G ( P + Q )ˆ γ νl ( P, P + Q ) ˆ G ( P ) i (45)By applying the WIs (34), one can see that Π µ , satisfies q µ Π µ ( Q ) = 0 and q µ Π µ ( Q ) = − i ∆. Hence, the full vertexagain satisfies the WI q µ ˆΓ µ ( P + Q, P ) = q µ ˆ γ µ ( P + Q, P ) + 2 i ( p + q j ∆ σ = ( p + q j h σ ˆ G − ( P + Q ) − ˆ G − ( P ) σ i . (46)This is consistent with the conservation law (42) of the response kernel. IV. CONCLUSION
In this paper, we construct a stress tensor linear response theory for BCS superfluid which satisfies the WI dueto the spacetime symmetry. Therefore we have established the local conservation laws of the energy and momentumin the superfluid phase at the BCS mean-field level. The pairing fluctuation effect is not included for now, whileone can treat our theory as a first step toward a the construction of a fully consistent stress tensor linear responsetheory including the fluctuations from the condensed pairs and the non-condensed pairs below T c . We emphasizethat our conserving approximation is quite different from the Φ-derivable theory in the mathematical form becausewe are focusing on the local form of the conservation laws or WI. The key point in constructing the full vertex is totreat the phase fluctuation of the order parameter as an independent field which is only constrained by the consistentgap equation. This method is closely parallel to what is used in establishing the U (1) current WI. The result ofthis paper can be applied to calculate the transport properties such as the viscosity since it can be derived eitherfrom the current-current correlations ? or the stress tensor-stress tensor correlations ? . To satisfy the WI is crucialto connecting these two different approaches.Yan He thanks the support by National Natural Science Foundation of China (Grants No. 11404228). Hao Guothanks the support by National Natural Science Foundation of China (Grants No. 11204032) and Natural ScienceFoundation of Jiangsu Province, China (SBK201241926). Appendix A: Correlation Functions for BCS Superfluids
In this appendix, we list all the detailed expressions of the response functions we have used in the main text. Q j ( ω, q ) = X k k j h (1 + ξ + k ξ − k − ∆ E + k E − k ) ( E + k + E − k )[1 − f ( E + k ) − f ( E − k )] ω − ( E + k + E − k ) − (1 − ξ + k ξ − k − ∆ E + k E − k ) ( E + k − E − k )[ f ( E + k ) − f ( E − k )] ω − ( E + k − E − k ) i (A1) Q j ( ω, q ) = X k k j h (1 + ξ + k ξ − k + ∆ E + k E − k ) ( E + k + E − k )[1 − f ( E + k ) − f ( E − k )] ω − ( E + k + E − k ) − (1 − ξ + k ξ − k + ∆ E + k E − k ) ( E + k − E − k )[ f ( E + k ) − f ( E − k )] ω − ( E + k − E − k ) i (A2) Q j ( ω, q ) = − Q j ( ω, q ) = − iω X k k j h ( ξ + k E + k + ξ − k E − k ) [1 − f ( E + k ) − f ( E − k )] ω − ( E + k + E − k ) − ( ξ + k E + k − ξ − k E − k ) [ f ( E + k ) − f ( E − k )] ω − ( E + k − E − k ) i (A3) Q j,nl ( ω, q ) = ∆ X k k j γ nl ξ + k + ξ − k E + k E − k h ( E + k + E − k )[1 − f ( E + k ) − f ( E − k )] ω − ( E + k + E − k ) + ( E + k − E − k )[ f ( E + k ) − f ( E − k )] ω − ( E + k − E − k ) i (A4) Q mj ( ω, q ) = ∆ X k γ mj ξ + k + ξ − k E + k E − k h ( E + k + E − k )[1 − f ( E + k ) − f ( E − k )] ω − ( E + k + E − k ) + ( E + k − E − k )[ f ( E + k ) − f ( E − k )] ω − ( E + k − E − k ) i (A5) Q j, l ( ω, q ) = X k k j k l ∆ ωE + k E − k h ( E + k − E − k )[1 − f ( E + k ) − f ( E − k )] ω − ( E + k + E − k ) + ( E + k + E − k )[ f ( E + k ) − f ( E − k )] ω − ( E + k − E − k ) i (A6) Q j ( ω, q ) = X k k j ∆ ωE + k E − k h ( E + k − E − k )[1 − f ( E + k ) − f ( E − k )] ω − ( E + k + E − k ) + ( E + k + E − k )[ f ( E + k ) − f ( E − k )] ω − ( E + k − E − k ) i (A7) Q j,nl ( ω, q ) = i X k k j γ nl ∆ ωE + k E − k h ( E + k + E − k )[1 − f ( E + k ) − f ( E − k )] ω − ( E + k + E − k ) + ( E + k − E − k )[ f ( E + k ) − f ( E − k )] ω − ( E + k − E − k ) i (A8) Q mj ( ω, q ) = − i X k γ mj ∆ ωE + k E − k h ( E + k + E − k )[1 − f ( E + k ) − f ( E − k )] ω − ( E + k + E − k ) + ( E + k − E − k )[ f ( E + k ) − f ( E − k )] ω − ( E + k − E − k ) i (A9) Q j, l ( ω, q ) = i ∆ X k k j k l ξ + k − ξ − k E + k E − k h ( E + k + E − k )[1 − f ( E + k ) − f ( E − k )] ω − ( E + k + E − k ) + ( E + k − E − k )[ f ( E + k ) − f ( E − k )] ω − ( E + k − E − k ) i (A10) Q j ( ω, q ) = − i ∆ X k k j ξ + k − ξ − k E + k E − k h ( E + k + E − k )[1 − f ( E + k ) − f ( E − k )] ω − ( E + k + E − k ) + ( E + k − E − k )[ f ( E + k ) − f ( E − k )] ω − ( E + k − E − k ) i (A11) Q mj,nl ( ω, q ) = X k γ mj γ nl h (1 − ξ + k ξ − k − ∆ E + k E − k ) ( E + k + E − k )[1 − f ( E + k ) − f ( E − k )] ω − ( E + k + E − k ) − (1 + ξ + k ξ − k − ∆ E + k E − k ) ( E + k − E − k )[ f ( E + k ) − f ( E − k )] ω − ( E + k − E − k ) i (A12) Q j, l ( ω, q ) = X k k j k l h (1 − ξ + k ξ − k + ∆ E + k E − k ) ( E + k + E − k )[1 − f ( E + k ) − f ( E − k )] ω − ( E + k + E − k ) − (1 + ξ + k ξ − k + ∆ E + k E − k ) ( E + k − E − k )[ f ( E + k ) − f ( E − k )] ω − ( E + k − E − k ) i (A13) Q j,nl ( ω, q ) = ω X k k j γ nl h ( ξ + k E + k − ξ − k E − k ) [1 − f ( E + k ) − f ( E − k )] ω − ( E + k + E − k ) − ( ξ + k E + k + ξ − k E − k ) [ f ( E + k ) − f ( E − k )] ω − ( E + k − E − k ) i (A14) Q mj, l ( ω, q ) = ω X k γ mj k l h ( ξ + k E + k − ξ − k E − k ) [1 − f ( E + k ) − f ( E − k )] ω − ( E + k + E − k ) − ( ξ + k E + k + ξ − k E − k ) [ f ( E + k ) − f ( E − k )] ω − ( E + k − E − k ) ii