Microscopic derivation of density functional theory for superfluid systems based on effective action formalism

Microscopic derivation of density functional theory for superfluid systems based on effective action formalism Takeru Yokota, Haruki Kasuya, Kenichi Yoshida, and Teiji Kunihiro Institute for Solid State Physics, The University of Tokyo, Kashiwa, Chiba 277-8581, Japan Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan Department of Physics, Kyoto University, Kyoto 606-8502, Japan

However it is known to be clumsy for a naive application of the original DFT to describe the systems with SSB such as magnetization or superfluidity solely with the particle density. So the order parameter as given by an anomalous density is included explicitly as an additional density to describe such a system with SSB [9][10][11][12][13][14][15][16][17][18][19][20].
A great challenge in the current study of DFT is to develop a microscopic and systematic framework for constructing the EDF from an inter-particle interaction. The notion of the effective action [21,22] has brought forth the nonperturbative framework for dealing with the quantum many-body problems, and the concepts and recipes developed and accumulated in the study of quantum-field theory (QFT) can be brought into the study of many-body systems. In particular, the effective action for composite fields [23][24][25] links DFT and concepts and recipes of QFT. Indeed the two-particle point-irreducible effective action [26] for the density field is identified with the freeenergy density functional containing all the information of not only the ground state but also the excited states [27,28]. Thus a promising strategy to overcome the perennial challenge has come out. In Ref. [29], what they call the inversion method was developed to calculate the ground-state energy nonperturbatively, and it was further extended to the Fermion system with superfluidity [30]. Along with this line of method, the notion of effective field theory is applied to construct the free-energy density functional including superfluidity [31][32][33][34].
Although the inversion method is certainly a well-founded and powerful nonperturbative scheme, the notions and techniques developed in QFT may not have been fully utilized. An established method is currently available for the calculation of effective actions, which is called the functional renormalization group (FRG) [35][36][37][38]. The FRG provides a nonperturbative and systematic procedure for the analyses of renormalization flows by solving one-parameter functional differential equations in a closed form of effective actions. Recently, the application to the effective action for composite fields aimed at an ab-initio construction of DFT, which we call the functionalrenormalization-group aided density functional theory (FRG-DFT), has been developed: After the proposal of the formalism in Refs. [39,40], several approximation schemes have been proposed, whose performances were tested in (0+0)-or (0+1)-dimensional toy models [41,42]. The application was extended to (1+1)-dimensional systems composed of nucleons in the canonical formalism [43] and the grand-canonical formalism [44]. And the excited states in (1+1)-dimensional systems were investigated [45]. Furthermore, analyses of more realistic systems such as homogeneous electron gases in (2+1) dimensions [46] and (3+1) dimensions [47] were recently achieved.
The purpose of this paper is to extend the FRG-DFT formalism as developed in Refs. [44,45] so as to be applicable to superfluid systems at finite temperature. To this end, we introduce the effective action for the particle-number and nonlocal pairing densities, whereby the Hohenberg-Kohn theorem in terms of the Helmholtz free energy is found to be readily established. We then derive the flow equation for the effective action by introducing a flow parameter λ with which the strength of the inter-particle interaction is adiabatically increased as λ is changed from 0 to 1. The differentiation of the effective action with respect to λ gives the coupled flow equations for the particle-number and pairing densities. Together with the variational equation which the equilibrium densities should satisfy, the exact expression for the KS potential generalized to including the pairing potentials is derived. The resultant KS potential has a nice feature that it expresses the microscopic formulae of the external, Hartree, pairing, and exchange-correlation terms, separately. As a demonstration of the validity of our formalism, we shall show that a lowest-order approximation to our flow equations with respect to the flow parameter reproduces the gap equation, the ground-state energy and pairing 2/26 density given in the BCS theory for the case of short-range interactions in the weak-coupling limit. Although we show the validity of our microscopic formulation by applying it to simple cases, an advantage of our exact FRG-DFT formalism lies in the fact that it provides ways to systematically improve the correlation part. This paper is organized as follows: In Sec. 2, we introduce the effective action for the particlenumber and pairing densities and derive the flow equation to calculate the effective action. From the flow equation, the exact expression of the KS potential is derived. In Sec. 3, we apply the lowestorder approximation to our formalism and compare our results with those given by the BCS theory. Section 4 is devoted to a brief summary and the conclusions.

Formalism
In this section, we develop a formalism to describe microscopically the superfluid systems at finite temperature in the framework of FRG-DFT, where the effective action given in the functionalintegral form is fully utilized. The formalism starts with construction of the effective action Γ[ρ, κ, κ * ] for the nonlocal pairing densities as well as particle-number ones, and then we show that Γ[ρ, κ, κ * ] can be identified with the free-energy density functional for which the Hohenberg-Kohn theorem holds. We then introduce a flow parameter λ with which the inter-particle interaction is switched on gradually as λ is varied from 0 to 1. Differentiating thus constructed effective action Γ λ [ρ, κ, κ * ] with respect to λ, we obtain the coupled flow equations for the particle-number and pairing densities. Thus the problem of deriving the governing equation of the superfluid systems in the context of DFT is reduced to solving the flow equations. Inherently in the way of the introduction of the flow parameter, the governing equation thus obtained actually leads to the exact expression of the KS potential. Moreover, because of the very exact nature of the governing equation thus obtained, the resultant equation can be a starting point for an introduction of various systematic approximations. We consider a system composed of non-relativistic fermions interacting via two-body interactions subject to the external potential at finite temperature 1/β with β being the inverse temperature. The action describing the system in the imaginary-time formalism is given as follows: Here, the short-hand notations ξ = (τ, x, a), where τ, x, and a represent the imaginary time, spatial coordinate, and index for internal degree of freedoms, respectively with ǫ being an infinitesimal constant. The fermion fields are represented by Grassmann variables ψ(ξ) and ψ * (ξ). The static external potential and the instantaneous two-body interaction are denoted by V(ξ) = V a (x) and U(ξ, ξ ′ ) = δ(τ − τ ′ )U aa ′ (x − x ′ ), respectively. The shifted coordinate ξ + ǫ τ = (τ + ǫ, x, a) is introduced so that the corresponding Hamiltonian becomes normal-ordered [44].

Effective action for superfluid systems
Let us introduce the effective action, which is to be identified with the free-energy density functional for superfluid systems. The construction of the effective action [21,22] starts by defining the generating functional of the correlation functions of density fields. In this work, we generalize the previous works [27,44,45] and introduce the generating functional for the correlation functions of 3/26 particle-number and (nonlocal) pairing densities as ) denotes the external source coupled to the respective density fields: with the following abbreviated notations for arguments We emphasize that the pairing fields are introduced as nonlocal ones with which the present formalism is applicable for describing not only the conventional spin-singlet pairing but any types of pairings as in Ref. [48] and furthermore the UV divergence, which appears in the case of local interactions and local pairings, can be avoided. Now we define the generating functional for connected correlation functions W[ J] by Then the effective action Γ[ ρ = (ρ, κ, κ * )] for the particle-number density ρ(ξ) as well as nonlocal pairing densities κ ( * ) (ξ, ξ ′ ) is given by a Legendre transformation of W[ J]: Here J sup [ ρ] denotes the source that gives the supremum of the quantity in the bracket of the first line and gives the expectation value of the density operatorρ i on account of Eqs. (2) and (5).

The Hohenberg-Kohn theorem for superfluid systems
Now that having defined the effective action by Eq. (6) for superfluid systems in a rather natural way, we are in a position to show one of the most important points in this article, i.e., a proof of the Hohenberg-Kohn theorem for superfluid systems with nonzero temperature [11] on the basis of the effective action without recourse to any heuristic arguments.
• Variational principle for Eq. ), which will be, however, taken to be zero in the end of calculation. We shall consider a generic case where the chemical potential µ a may depend on the internal degree of Fig. 1 Schematic pictures of the κ-dependence of I in the case of (a) no external potential and (b) small nonzero external potential, where the ρ-dependence is not depicted. Two different equilibrium pairing densities are denoted by κ (1,2) ave . δ is an infinitesimally small but nonzero number. Without the external potential, κ ave is not determined uniquely in the variational way since I has the same value in the region shown in (a). With a small external field being applied, the equilibrium density [κ (1) ave in (b)] may be obtained variationally since the flat region is slightly tilted. freedom a; however, notice that the possible spatial dependence of µ a can be neglected without loss of generality because it can be absorbed away into the redefinition of V a (x). Now under the condition that the external potential δ iρ V(ξ ρ ) − µ i (ξ i ) is applied to the system, consider the variational problem of the following functional: ). The solution ρ ave of the variational equation and thus one obtains by the use of Eqs. (2) and (7). Equation (9) shows that ρ ave is the average of the density in the equilibrium state in the presence of the external potential A remark is in order here: Although the vanishing limit is formally taken for µ κ ( * ) (ξ κ ( * ) ), this limit should be taken with a care when SSB is concerned; it should be interpreted as 0+ in the mathematical notation meaning an infinitesimally small but nonzero value to obtain ρ ave variationally. Figure 1 shows a schematic picture of I[ ρ] solely as a function of κ, say, with ρ fixed. If we have two solutions κ (1) ave and κ (2) ave as the pairing density in equilibrium, I[ ρ] with µ κ ( * ) = 0 is actually flat in the interpolated densities [κ = ακ (1) ave is a convex functional as a result of the Legendre transformation. This prevents one from determining the equilibrium density uniquely by a variational method. As long as µ κ ( * ) is kept small but nonvanishing, however, the flat region of I[ ρ] is slightly tilted, then one can uniquely identify the equilibrium density in the variational method.
Furthermore, it follows immediately from Eq. (9) that the external potential determines the equilibrium density uniquely. On the other hand, it is easily demonstrated that the equilibrium 5 where we have used J i . This argument is based on the fact that the functional is strictly convex as seen in Fig. 1(b). Therefore, one has a one-to-one mapping between the external potential δ iρ V(ξ ρ ) − µ i (ξ i ) and the equilibrium density ρ ave . In particular, the existence of the map from ρ ave to δ iρ V(ξ ρ ) − µ i (ξ i ) can be seen more clearly through the fact that Γ[ ρ] is linearly dependent on the external field, which we now show below. • One-to-one map between the external fields and the densities.
From Eq. (6), we obtain To show this, let us assume that two sets of the external potential and the chemical potential denoted by (V 1 , µ 1 ) and (V 2 , µ 2 ), respectively, give the same densities ρ ave . By use of the variational equation of I[ ρ] and Eq. (11), however, we obtain Subtracting these equations from each other, we have the equality which proves that Thus the Hohenberg-Kohn theorem for superfluid systems with nonzero temperature, which was first demonstrated in Ref. [11], is established in terms of the effective action. Furthermore it is also 6 In fact, in the limit of µ κ ( * ) → 0, Eq. (13) at ρ = ρ ave becomes because of the relation J sup [ ρ ave ] = µ and Eq. (6). In terms of words, since Z[ µ] is the grand partition function in the presence of the chemical potential µ, where {E n } is the energy eigenvalues of the system and satisfies E gs = E 0 < E 1 < · · · . Therefore, we can identify the EDF with the effective action as In particular, the universal part of the EDF is obtained by substitution of Γ with Γ| V=0 . This concludes the generalization of the correspondence between the EDF and the effective action to the superfluid systems [27,28,44].

Flow equation
The effective action Γ[ ρ] contains the whole contents of physics of the quantum system in a compact form of the functional integral. Although the compact and exact form of the action allows us to define consistent approximations, the problem is developing a computational method to perform the functional integral in a systematic and hopefully exact way. In this subsection, we are going to develop a practical method to calculate the effective action in such a manner. A wise way originally developed by Wegner [35] and Wilson [36] is to convert the functional integral to a differential equation, which is called a renormalization-group (RG) equation or flow equation, and solve the equation with some systematic approximation being adopted. In the case of the effective action solely for the particle-number density and accordingly without SSB, a calculational framework is developed based on an RG/flow equation in Refs. [39,40]; see also Refs. [44,45]. Here, we generalize this framework to the case with SSB for the calculation of Γ[ ρ] containing the pairing density fields.
Following Refs. [39,40,44,45], we introduce a regulated two-body interaction U λ (ξ, ξ ′ ) with a flow parameter λ and define the regulated action as where U λ (ξ, ξ ′ ) is defined by with a constraint so that the action evolves from the free S λ=0 to the fully-interacting one S λ=1 as λ runs from 0 to 1. With this action, the λ-dependent generating functionals and effective action are defined in the same manner as in Sec. 2.1: where Now utilizing the λ dependence of the action (19), the flow equation of Γ λ [ ρ] can be obtained simply by differentiating it with respect to λ [39,40,[43][44][45]. It turns out, however that the most convenient way to obtain the flow equation where with O[ψ, ψ * ] being some functional of ψ and ψ * . Noting that U λ (ξ where an infinitesimal number ǫ ′ τ plays the same role as ǫ τ but its zero limit should be taken after that of ǫ τ ; see Appendix in Ref. [44]. From Eq. (18), one can immediately see .
Using these relations, one finds that Eq. (21) takes the form of the flow equation for W λ [ J] as follows: . (27) By taking the second derivative of Eq. (6) with respect to ρ, we have One finds that Eq. (27) is the inverse of Eq. (28): Here, and By taking the derivative of Eq. (6) with respect to λ, we obtain Substituting Eqs. (20), (29), and (31) into Eq. (26), we eventually arrive at the formal flow equation Now Eqs. (25) and (29) together with Eq. (20) tell us that ( δ 2 Γ λ δ ρδ ρ ) −1 ρρ gives essentially the correlation functions of the densities ρ, which collectively denote the particle-number and pairing densities. To see more explicitly the physical contents of the flow equation (32), in particular, the contributions of the pairing densities, we find it convenient to separate the λ = 0 part of the correlation function 9/26 where the particle operators are free ones: By applying the Wick's theorem, we obtain where we have used the equality ρ i (ξ i ) λ, J sup,λ [ ρ] = ρ i (ξ i ) which follows from Eqs. (24) and (20). One sees that the first term given by a product of the pairing condensates in Eq. (33) is the counter part to that of the particle-number densities given in the first term of Eq. (32), while the terms in the second line and the remaining order-λ terms are identified with the exchange and correlation terms, respectively. Thus one finds that it is natural to represent Eq. (32) in a form where the terms with different physical significance are written separately as, where G (2) xc,λ [ ρ](ξ ρ 1 , ξ ρ 2 ) denotes the exchange-correlation term given by with We call G (2) x [ ρ](ξ ρ 1 , ξ ρ 2 ) and G (2) c,λ [ ρ](ξ ρ 1 , ξ ρ 2 ) the exchange and correlation terms, respectively. We note that the flow equation gives the systematic way to calculate the higher-order contribution for Γ λ [ ρ], or G (2) c,λ [ ρ], since it is not only an exact but also closed equation for Γ λ [ ρ]. Moreover, it also provides us with the basis for defining consistent approximations, say, respecting symmetries, if necessary. In this respect, some approximation schemes developed for the study of the normal phase such as the vertex expansion [39,40] is expected to be utilized for the determination of Γ λ [ ρ] also in the case of superfluid systems. 10/26

Expressions of Helmholtz free energy and ground-state energy
Having derived the flow equation (34) for the effective action Γ λ [ ρ], we here give a reduced formula for the Helmholtz free energy which is given in terms of the effective action.
The Helmholtz free energy given in Eq. (13) now takes the following form where ρ ave is the equilibrium density of the fully-interacting (λ = 1) system in the presence of a given chemical potential µ: which is nothing but Eq. (9) The ground-state energy is obtained by taking the zero temperature limit in Eq. (40): where T [ ρ gs ] := lim β→∞ Γ λ=0 | V=0 [ ρ ave ]/β with ρ gs ≔ lim β→∞ ρ ave being the ground-state density is interpreted as the kinetic energy. Notice that T [ ρ gs ] is neither the kinetic energy of a free gas nor that of the interacting system, but that of the non-interacting system with an appropriate potential, as may be utilized in KS formalism. Indeed, we shall show in the following subsection that the KS theory for superfluid systems at finite temperature emerges naturally in the present formalism.

The Kohn-Sham potential
In our formulation of the FRG-DFT theory for superfluid systems, the starting action with λ = 0 is that for a non-interacting system with the external forces ensuring the densities. Thus it is natural that J sup,λ [ ρ] is related to the KS potential V KS [ ρ] that is generalized so as to incorporate the anomalous or pair potentials besides the normal ones: The KS potential is defined by the derivative of the EDF subtracted by the non-interacting kinetic energy part at zero temperature. In terms of the effective 11/26 action at finite temperature, the kinetic energy part T KS [ ρ] may be identified with Γ λ=0 | V=0 , and then we naturally arrive at the definition of the KS potential as By differentiating Eq. (19) with respect to the densities, we obtain Using this relation and Eq. (11) for Γ = Γ λ=0 , Eq. (42) is rewritten as follows: The KS potential determines ρ ave for given chemical potentials µ: which is equivalent to the variational equation in the KS DFT [49] and determines ρ ave from the KS potential (42). The present formulation does not refer to any single-particle orbitals for establishing the KS DFT. Note that Eq. (45) does possess the property of a self-consistency in terms of the densities. In fact the equilibrium density of the full-interacting system with a given chemical potential µ is also the equilibrium density of the non-interacting system under the external potential µ − V KS [ ρ ave ], as we will now show explicitly: First of all, we notice that Eq. (45) is also expressed as Eq. (44). Then inserting this relation into Eqs. (18) and (20), one finds that Eq. (20) tells us that the equilibrium density ρ ave in Eq. (39) also has the following expression where Z λ | V=0 [ J] is defined by the action S λ | V=0 [ψ * , ψ] with no external potential in Eq. (18). Note that Eq. (46) is a self-consistent equation for ρ ave since V KS [ ρ ave ] in the right-hand side depends on ρ ave in the left-hand side, as noted before. Next, let us reduce the exact expression of the central quantity V KS [ ρ] in our formalism in a more calculable form from the flow equation (34). After differentiating Eq. (34) with respect to ρ and then integrating it with respect to λ, we have 12 /26 where U λ (ξ κ ( * ) := (ξ, ξ ′ )) = U λ (ξ, ξ ′ ) and Eq. (43) has been used. With use of Eq. (44), we finally arrive at the following expression for V KS [ ρ]: It is noteworthy that V i KS [ ρ] is now naturally decomposed into four terms, which denote the external, Hartree, pairing, and exchange-correlation terms, respectively. In particular, a microscopic expression of the exchange-correlation term V i KS,xc [ ρ] is explicitly given in terms of G (2) xc,λ [ ρ]: We emphasize that our FRG-DFT formalism applied to superfluid systems has naturally led to the self-consistent equation for the densities that are kept the same between the non-interacting and interacting systems together with the microscopic expressions of the KS potential and the kinetic energy without recourse to notions of single-particle orbitals. It is also noted here that the effective action formalism with the inversion method was applied to a microscopic construction of the KS potential for a normal system [50] and a superfluid system of the local singlet pairs [33]. However, the exact expression of the KS potential is not derived in the inversion method. Furthermore, the single-particle basis are introduced in developing the KS DFT. 1 In the next section, we are going to show that Eq. (40) together with Eqs. (45) and (48) leads to the ground-state energy in the Hartree-Fock-Bogoliubov approximation, which is furthermore reduced to the well-known gap equations of BCS theory in the case of short-range interaction.

Analysis in the lowest order approximation: neglect of correlations
We demonstrate that our formalism actually reproduces the results of the mean field theory at the lowest-order truncation, implying that G (2) xc,λ does not have λ dependence: G (2) xc,λ ≈ G (2) x .

The case of a generic interaction
Applying the lowest-order truncation to Eq. (48), we have 1 In Ref. [34], it was anticipated that the KS DFT may be constructed without single-particle orbitals in the FRG-DFT.

13/26
for the KS potential. On account of Eq. (36), δG (2) x [ ρ](ξ ρ 1 , ξ ρ 2 )/δρ i 1 (ξ i ) is represented as follows: Here, the following notation has been introduced: , λ=0 are expressed in terms of the density correlation functions of the non-interacting system, which can be further reduced to simpler forms using the Wick's theorem. When the superfluidity is not considered [43,44], Γ (n≥2) λ are also written in terms of the connected density correlation functions. We are going to extend them so as to incorporate the anomalous density-correlation functions as where G (n) λ,i 1 ···i n (ξ i 1 1 , · · · , ξ i n n ) is the connected correlation function given by , Using these relations, the last term in the right-hand side of Eq. (50) is rewritten as follows:

14/26
By use of the Wick's theorem, is evaluated as follows: Here, we have introduced in which the right-hand side is defined by Eq.
where we have used the equal-time commutation relation: Substituting Eq. (56) into Eq. (54), we have 15/26 Although the last term in the right-hand side may not be possible to represent in terms of the density correlation functions in general, it can be further simplified for an arbitrary U λ . In the next subsection, we shall discuss the case with short-range interactions and then show that the last term of Eq. (57) is approximately represented in terms of the density correlation functions. Finally, let us discuss the Helmholtz free energy. Applying the approximation G (2) xc,λ ≈ G (2) x to Eq. (40), we have From Eqs. (33) and (36), Then the integral with respect to ξ ρ 1 and ξ ρ 2 in the third line of Eq. (58) is evaluated as follows: where we have used the equal-time commutation relation: . (60) At the zero temperature limit β → ∞, it gives the ground-state energy: .
The third term in the right-hand side is interpreted as the Hartree-Fock energy; The fourth term is the contribution from the pairing condensate. This formula of the ground-state energy is identical to that given in the Hartree-Fock-Bogoliubov approximation [51][52][53]. 16

The case of short-range interactions in the weak coupling
Let us take the case where the inter-particle interaction is short-range in the weak coupling limit of the paring force as in the BCS theory [54]. Then it will be found that Eq. (57) is reduced to a simple and familiar form.
In the case of short-range interactions, the integral in the last term of Eq. (57) is dominated by the contribution from the region where x 2 is close to x 1 . Therefore, we can make the following approximation for short-range interactions that holds exactly for the zero-range interaction: where we have used and introduced ρ a 1 a 2 (τ, x) = ψ * (τ + ǫ ′ τ , x, a 1 )ψ(τ, x, a 2 ) ρ . We evaluate ρ a 1 a 2 (τ, x) in the weak coupling limit leading to small κ in comparison with the particle density ρ. In the non-interacting case λ = 0, the vanishing external field J κ ( * ) sup,λ=0 [ ρ] = 0 gives κ ( * ) = 0. Conversely, we have lim κ,κ * →0 J κ ( * ) sup,λ=0 [ ρ] = 0 since κ ( * ) uniquely determines the external field. Therefore, by use of Eqs. (55) and (22), ρ a 1 a 2 (τ, x) is approximated in the weak coupling as Notice that there is no term causing a mixing between ψ ( * ) (τ, x, a) with different a's in the exponent S λ=0 [ψ, ψ * ] − ξ J ρ sup,λ=0 [ ρ](ξ)ρ(ξ) any more. Thus we have ρ a 1 a 2 (τ, x) ≃ 0 for a 1 a 2 in the weak coupling limit. In other words, ρ a 1 a 2 (τ, x) with a 1 a 2 is negligible in comparison with the diagonal component ρ a 1 a 1 (τ, x) = ρ(τ, x, a 1 ) in the weak coupling limit or small κ ( * ) . Thus, we can make the approximation for ρ a 1 a 2 (τ, x) as follows: Then, Eq. (61) is reduced as follows, 17/26 Applying this approximation to Eq. (57), we arrive at the following expression of the basic equation for the KS theory for the case of short-range interactions in the weak coupling as

Homogeneous systems
In this subsection, we are going to discuss a homogeneous system with V(ξ) = 0 in Eq. (1), and show that Eq. (65) together with Eq. (45) is reduced to the well-known gap equation in the BCS theory [55]. For simplicity, we consider fermions with spin but not any other internal degrees of freedom: we assume that the spin is saturated in the system, and accordingly the chemical potential µ has no spin dependence. We also assume that the two-body interaction is a spin-independent short-range one U ab (x) = U(x) and consider the weak-coupling limit so that the approximated equation Eq. (64) is applicable. Although we treat the spin-singlet condensate in the present investigation exclusively in the present article, the extension to the spin-triplet case is straightforward by considering the vector pairing densities.
Here, we introduce the following notation: Owing to the translational symmetry, the coordinate dependence of V ρ KS (ξ ρ ) is simplified as follows: Here, ρ ave = ρ ↑ ave + ρ ↓ ave is the total particle-number density. Other components for κ ab ave (τ, x) and V κ KS,ab (τ, x) are assumed to be zero since we focus on the systems with the spin-singlet condensate. For convenience for the description of the homogeneous system, let us move to the momentum representation. Then Eq. (65) is reduced to the following two equations for the particle and pairing 18/26 densities, respectively, Here, we have introduced the notations as P = (p 0 , p) and Q = (q 0 , q) with the Matsubara frequencies p 0 and q 0 , the spacial momenta p and q, andṼ κ,s KS (P),Ũ(p) andκ s ave (P) are the Fourier transforms of V κ,s KS (τ, x), U(x), and κ s ave (τ, x), respectively. The following abbreviation is also used for the Q-'integrations': Notice that the right-hand side of Eq. (66) and henceṼ κ,s KS (P) in the left-hand side are independent of p 0 . Thus we write asṼ κ,s KS (P) =Ṽ κ,s KS (p) from now on. Now we calculate the equilibrium densities in the KS system via the variational equation (45) with the obtained KS potential Eq. (66). However, we have already seen that the solution ρ ave of Eq. (45) satisfies the self-consistent equation Eq. (46), which is now simplified to Furthermore, since the exponent in the integral becomes a bilinear form when Eq. (66) is inserted, we can perform the path integral analytically to give where p = dp/(2π) 3 , ε(p) = p 2 /2 + V ρ KS − µ, and E(p) = ε(p) 2 + |2Ṽ κ,s KS (p)| 2 , as is detailed in Appendix. Substitution of Eq. (68) into the second equation of Eq. (66) in turn leads to Equation (68) implies thatṼ κ,s KS (p) is identified with the energy gap ∆ s (p) in the BCS theory as ∆ s (p) = −2Ṽ κ,s KS (p) * . With this identification, we find that Eq. (69) is nothing else than the celebrated gap equation in the BCS theory [54].
As is also shown in the Appendix, the Helmholtz free energy is calculated from Eq. (60) to yield where V = dx is the volume of the system.

Conclusion
On the basis of the effective action formalism, we have developed a generalized density-functional theory (DFT) for superfluid systems at finite temperature in a rigorous way. The rigorous formulation is combined with the renormalization group method by introducing a scale parameter that is multiplied to the inter-particle potential, which is reminiscent of the adiabatic connection method. A possible difficulty that the spontaneous symmetry breaking (SSB) is not described when symmetrybreaking term is absent is nicely circumvented by a suitable choice of the external fields by making use of the functional renormalization group (FRG). Then we have established the Hohenberg-Kohn theorem with the SSB being taken into account. By introducing flow-parameter-dependent source terms appropriately fixing the ground-state densities during the flow, we have derived the flow equation for the effective action for the particle-number and nonlocal pairing densities where the flow parameter λ runs from 0 to 1, corresponding to the non-interacting and interacting systems, respectively.
Integrating this flow equation and using the variational equation for the thermal equilibrium densities, we have arrived at the exact self-consistent equation to determine the Kohn-Sham potential V i KS [ ρ] generalized to including the pairing potentials at finite temperature without introduction of single-particle orbitals. The resultant Kohn-Sham potential has a nice feature that it expresses the microscopic formulae of the external, Hartree, pairing, and exchange-correlation terms, separately. In particular, the exchange-correlation term V i KS,xc [ ρ] is given explicitly in terms of the corresponding density-density correlation functions G (2) xc,λ [ ρ]. As a demonstration of the validity of our formulation as a microscopic theory of DFT for superfluidity, we have shown that our self-consistent equation with Kohn-Sham potential leads to the ground-state energy of the Hartree-Fock-Bogoliubov theory when the correlations are neglected. And the self-consistent equation is further reduced for short-range interaction and in the weak-coupling limit to reproduce the gap equation of the BCS theory for the spin-singlet pairing.
Although we have shown the validity of our microscopic formulation of DFT based on FRG by applying it to simple cases, an advantage of our exact FRG-DFT formalism lies in the fact that it contains the equation (37) with which the correlation part can be improved systematically. One of the most effective systematic schemes used in the FRG-DFT is the vertex expansion, where the Taylor expansion around the equilibrium densities are applied to the flow equation for the effective action. In fact, it is shown [45] for the system without superfluidity that the vertex expansion only up to the second order gives an approximate scheme superior to the random phase approximation (RPA). Therefore it is naturally expected that a simple vertex expansion applied to the basic equation in the present work would lead to an approximation scheme that incorporates correlations beyond those given by the quasi-particle RPA. We hope to report on such an attempt near future. 20/26 Our formulation was given to general cases, i.e. any spatial profiles and the dependence of the interaction on the internal degrees of freedom are not imposed. Such a general formulation may be helpful when applying to the superfluidity in various systems. Introduction of the nonlocal particle density is an interesting extension of the present formalism, with which we can describe the particlehole and particle-particle correlations on the same footing. Furthermore, our general formulation for nonlocal pairing should provide us with a sound basis for describing pairing phenomena with any symmetry such as spin-triplet superfluidity and/or the paring with spin-orbit coupling incorporated. We hope that there will be a chance to report on such investigations.
To calculate the Matsubara frequency summation in Eq. (A13), let us introduce a complex function and a weighting function It is observed that the product h(−iz)g(z) decays sufficiently fast as |z| → ∞, and h(−iz) has simple poles at z = ± 2 ǫ sinh −1 ǫE In the second to last line, we have taken the limit ǫ → 0 implicitly. Consequently, we arrive at In the same way, the equilibrium spin-singlet pairing density is calculated as follows: