Current-density functional theory for bosonic superfluids

A finite-temperature current-density functional theory for bosonic superfluids (sf-CDFT) in the thermal equilibrium state is proposed herein. In the sf-CDFT, hydrodynamic physical quantities, such as particle number density, current density, and the order parameter of the Bose–Einstein condensation, are chosen as the basic variables. This theory enables the simultaneous reproduction of the particle number and current densities of both the superfluid and normal fluid components with incorporating effects of the interaction between these components. Specifically, these components are determined by solving two single-particle equations, i.e., the Gross–Pitaevskii–Kohn–Sham and Kohn–Sham equations. Furthermore, using the continuity equation of superfluids, we present the sum rule for the exchange-correlation energy functional of the sf-CDFT, which is useful for developing the approximate form.

Considering the theories of superfluids, the two-fluid model [15,16], which is known as the standard model, in which a fluid consists of superfluid and normal fluid components, is widely used for describing the macroscopic behavior of superfluids. Among the microscopic theories, the Gross-Pitaevskii (GP) theory [17][18][19] and Bogoliubov theory [20,21] are widely known and have been used for the investigation of superfluids. The GP theory is suitable for treating a weakly interacting boson system in the ground state. In this theory, the behavior of the order parameter of the BEC can be obtained by solving the GP equation [17][18][19], which allows the discussion of the dynamics of the particle number and current densities of the superfluid in the ground state [5]. The Bogoliubov theory was developed for dealing with the excited states of superfluids consisting of weakly interacting bosons. In this theory, the quasi-particle description is applied to the excited states of superfluids, assuming that their ground state is given by the solution of the steady-state GP equation. Thus, the Bogoliubov theory enables an approximate calculation of the particle number and current densities in the excited states. Because the set of quasiparticle excitations may be considered a normal fluid component in a superfluid, the particle number and current densities of the normal fluid component can be calculated using the Bogoliubov theory. By subtracting the particle number and current densities of the normal fluid component from the whole ones, the particle number and current densities of the superfluid component can be obtained.
In the GP and Bogoliubov theories, the interactions acting on bosons are assumed to be of the delta function type, which leaves room for improvement for systems composed of bosons with relatively strong interactions, such as liquid helium-4. Indeed, for liquid helium-4, some experimental results cannot be described by the GP or Bogoliubov theories [22,23]. Furthermore, the aforementioned conventional theories cannot simultaneously describe hydrodynamic physical quantities, such as the density and current density, of the superfluid and normal

Formulation OF sf-CDFT
Physical quantities that describe superfluid phenomena include the particle number density, current density, and the order parameter of the BEC. In this section, a density functional theory for bosonic superfluids that enables the prediction of these physical quantities is proposed. ) denotes the interaction potential between two bosons. In equation (1), v r ( ) denotes an external potential, for example caused by a trapping potential, and A r ( ) denotes a gauge field, for example caused by rotating the vessel [48]. The Hamiltonian (1) is similar to the Hamiltonian considered in the current-density functional theory for fermion systems [29]. For later convenience, we rewrite equation ( where T, W , n r ( ) and j r ( ) denote the kinetic energy, boson-boson interaction, particle number density and current-density operators, respectively, which are given by While formulating the conventional CDFT for fermion systems [29] as well as the conventional DFT [24,25], physical quantities coupled with external potentials can be chosen as the basic variables. For the aforementioned Hamiltonian, n r ( ) and j r ( ) are coupled with external potentials v r ( ) and A r ( ) so that the particle number density and current density can be chosen as basic variables. It is useful to reproduce not only the particle number density and the current density but also the order parameter of the BEC in describing the properties of the superfluid. In the present theory, the order parameter of the BEC is chosen, in addition to the particle number density and the current density, as the basic variables. For this aim, instead of applying the extended constrained-search theory [30,31] to the system, artificial fields D r ( ) and D r ( ) ⁎ are introduced as a mathematical device [39] in order to treat the order parameters of the BEC y r ( ) and y r ( ) † as basic variables. These artificial fields are coupled with y r ( ) and y r ( ) † in the Hamiltonian and are eventually set to zero. The Hamiltonian with these artificial field terms is given by The proof of equation (12) is given in appendix B.

Noninteracting reference system
Similar to DFT and its extensions [24][25][26][27][28][29][30][31][32][33][34][35], we introduce a noninteracting reference system in which the basic variables of the real system at the thermal equilibrium state are reproduced. The Hamiltonian of the reference system H Ŝ is given by are the effective potentials, which are determined later so that the thermal equilibrium values of basic variables in the reference system coincide with those in the real system. This reference system touches the reserver whose temperature and chemical potential are q and m, respectively.
The essential point about the reference system is that H Ŝ contains the effective potential term D , Ŝ which explicitly breaks the conservation of the boson number. In other words, we intentionally prepare a device that can express the symmetry-breaking state of the U(1) gauge for the reference system. This device is similar to that of the Bogoliubov theory, in which an approximation breaking the U(1) gauge symmetry is made for the Hamiltonian [20,21].
Similarly, as in the case of the real system, we can prove the HK theorem for the reference system. That is, the class of effective potentials  In equation (18), n r , Seq ( ) j r , Seq ( ) y r Seq ( ) and y* r Seq ( ) denotes the thermal equilibrium values of the basic variables in the reference system and are given by where X again represents n, j , y, or y . † The effective potentials are determined by using HK theorem for the real system (equation (12)) and that for the reference system (equation (18) The derivation of equations (20)-(23) is given in appendix C. By the analogy with the vector potential induced from the motion of charged particles, the second term of equation (21), which is hereafter denoted by A r , xc ( ) may be considered as a sort of vector potential induced by the fluid current. The vector potential induced from the motion of charged particles is attributable to the Coulomb interaction between charged particles and is coupled with the current of the charged particles, while A r xc ( ) is attributable to the interaction between bosons, as shown in equation (21), and is coupled with the current density of the fluid, as shown in equation (13). Because the effective potential D r S ( ) is coupled with the order parameter of the BEC, as shown in equation (13), y r eq ( ) is induced by D r S ( ) in the reference system. This is similar to the case where the order parameter of the superconducting state is induced by the effective pair potential in the reference system of the DFT for superconductors [40,41].

Single-particle equations and basic variables
As mentioned previously, the Hamiltonian H Ŝ does not commute N; hence, the total particle number of the reference system is not the conserved quantity. Thus, calculating traces such as equation (19) is challenging. If H Ŝ had been quadratic with respect to the field operator of bosons, it could have been diagonalized with respect to the boson quasi-particle via the Bogoliubov-Valatin transformation [21,49,50]. However, H Ŝ includes the potential term D , Ŝ which is first order with respect to the field operator. To resolve this issue, we introduce the following operator and classical field with reference to the Bogoliubov theory [20,21]: where Y r ( ) is a classical field with a complex function, and y r ( ) is the operator indicating the difference between y r ( ) and Y r .
( ) In the Bogoliubov theory, the solution of the GP equation is used for Y r ( ) so that the first-order terms with respect to y r ( ) vanish in its approximate Hamiltonian [20]. As for the present case, we determine Y r ( ) so that the first-order terms with respect to y r ( ) vanishes in H . Ŝ This procedure shall be explained below. Substituting equation (25) in equation (13) The first and second terms on the right-hand side of equation (26) are the zeroth-order terms of y r ( ) and y r .

( ) †
The third and fourth terms are first-order terms, and the last term is a second-order term. Here, we shall determine the complex functions Y r ( ) and Y r ( ) ⁎ to satisfy the following equations: are determined by equations (28) and (29), the first-order terms (third and fourth terms of equation (26)) are eliminated. Note that equations (28) and (29) are equivalent, because h Ŝ is a Hermitian operator 3 . Hereafter, equation (28) will be referred to as the Gross-Pitaevskii-Kohn-Sham (GP-KS) equation.
Considering these functions, we can rewrite equation (26) as the quadratic form of y r ( ) and y r : The single-particle Hamiltonian h Ŝ is a Hermitian operator, since both A r xc ( ) and v s (r) are given by equations (21) and (20), respectively, and are therefore real functions.
where the zeroth order terms (first and second terms on the right-hand side of equation (26)) are replaced by K 0 for convenience. Thus, the aforementioned difficulty can be avoided by determining Y r ( ) and Y r ( ) ⁎ satisfying equations (28) and (29), respectively.
Next, we expand y r ( ) and y r ( ) † using the eigenfunctions of h . Ŝ Suppose the eigenfunctions and eigenvalues of h Ŝ are denoted by f r i ( ) and e , i respectively, i.e., where a î and a î † denote the annihilation and creation operators, respectively, of a certain boson particle with the state given by f r .
Finally, we express the basic variables using the eigenvalues and eigenfunctions of the KS equation (equation (31)) and GP-KS equation (equation (28)). By solving the KS equation and using the resultant solutions, the Hamiltonian of the reference system is written in the diagonalized form given by equation (33). Using equation (33), the basic variables of the reference system can be calculated by using equations (5), (6), (19), (25), (31) and (32). We have It should be noted that the classical field Y r , ( ) which is introduced to eliminate the first-order terms of y r ( ) and y r , ( ) † is just the order parameter of the BEC. Therefore, the GP-KS equation (equation (28)) is recognized as the equation that should be fulfilled by the order parameter of the BEC. The order parameter of the BEC corresponds to the eigenfunction of the first-order reduced-density matrix whose occupied number (eigenvalue) is O(N). More specifically, Y r ( ) is equal to such an eigenfunction multiplied by the square root of its eigenvalues. Therefore, the first term on the right side of equation (34) can be recognized as the density of the superfluid component. Because n r eq ( ) denotes the total density of the system, the remaining part on the right side of equation (34), i.e., the second term, corresponds to the density of the normal fluid component. Because n r eq ( ) and Y r ( ) are reproduced correctly in the reference system, the remaining part can be reproduced correctly. Similarly, the first term on the right side of equation (35) can be recognized as the current density of the superfluid component, and the second term as the current density of the normal fluid component.
If we denote the density and current density of the superfluid component by n r sf ( ) and j r sf ( ) and those of the normal fluid component by n r nf ( ) and j r , nf ( ) equations (34) and (35) can be rewritten as = + n n n r r r, 3 8 eq sf nf Taking the additional current term caused by A r ( ) into consideration, the total current density of the real system J r eq ( ) is given by  23)). Thus, the sf-CDFT enables a simultaneous treatment of the superfluid and normal fluid components. This is one of the significant features of the present theory.
Next, we shall comment on the calculation procedure. The calculation procedure is shown in figure 1 and is as follows: (i) An input set of basic variables is generated. (iii) The GP-KS equation (equation (28)) is solved using the effective potentials presented in step (ii). (iv) The KS equation (equation (31)) is also solved using the effective potentials presented in step (ii).
(v) By using solutions of both the GP-KS and KS equations (steps (iii) and (iv)), the basic variables are calculated using equations (34)-(37).
(vi) By comparing the thus-obtained set of basic variables with the trial set of basic variables (step (i)), selfconsistency is verified. If the variables are inconsistent within some accuracy, the calculations are restarted from (i), and the input set is changed. The calculations are repeated until the aforementioned selfconsistency is attained. After the self-consistency is verified, the GP-KS and KS equations give the same particle density and current density for both the superfluid and normal fluid components as those in the previous iteration. This means that each equation is solved in conjunction.
Thus, using the self-consistent solutions of the GP-KS and KS equations, the hydrodynamic physical quantities of superfluid and normal fluid components can be obtained simultaneously.
At the end of this section, we summarize some of the distinctive features of the present scheme.   3. Discussion

Condition for the exchange-correlation energy functional
We first derive the continuity equation related to the GP-KS equation. Taking the complex conjugate of the equation, which is the multiplication of equation (28) and Y r , ( ) and subtracting its original one from the conjugate equation, we get By following the same procedure for the KS equation, the following equation can be derived: Adding both sides of equation (47) and those of equation (48) and using equations (38) and (39) using equations (21) and (22). Equation (50) is considered the relation that should be satisfied by y y j F n, , , .
The approximate form of the xc energy functional in the original DFT or its extensions has been successfully developed using relations to be fulfilled by the xc energy functional as the restrictive conditions [51][52][53][54][55]. For instance, the generalized gradient approximation [51][52][53] for the xc energy functional of the DFT and the vorticity expansion approximation for the xc energy functional of the CDFT [54,55] have been developed so far. Therefore, equation (50) may be used to devise an approximate form of y y j F n, , , . In contrast, the sf-CDFT can reproduce the particle number density, current density and order parameter of the BEC. As mentioned in section 2, the density and current density can be expressed using two terms. One is written using the order parameter given as a solution of the GP-KS equation, and the other is given by the solutions of the KS equation. The former represents the contribution of the superfluid component, and the latter represents that of the normal fluid component. Therefore, the sf-CDFT enables the calculation of the hydrodynamic physical quantities, such as the particle number density and current density, corresponding to both the fluid components simultaneously.
We can also show that the solution of the GP-KS equation Y r ( ) includes some additional terms other than those in the solution of the GP equation y r . GP ( ) To compare Y r ( ) with y r , GP ( ) we can rewrite equation (28) to obtain where p is the momentum operator, and v r H ( ) and v r xc ( ) are the second and third terms of the right-hand side of equation (20), respectively, i.e.,   (28)) is reduced to the GP equation for the system under the gauge field A r .
( ) Therefore, equation (55) is reduced to Y y = r r.
GP ( ) ( ) In contrast, the second term on the right-hand side of equation (55) is considered an additional term attributable to the effects of the xc energy functional and classical interaction between the bosons in the superfluid and normal states. Thus, while the order parameter in the GP equation is determined by incorporating the effect of the interaction only between bosons in the superfluid state, the order parameter reproduced in the sf-CDFT is determined by incorporating the interaction terms of the bosons in the superfluid and normal states via v r H nf , ( ) in addition to the interaction only between bosons in the superfluid state. Furthermore, the order parameter reproduced in the sf-CDFT contains the xc effects via v r , xc ( ) A r xc ( ) and D r s ( ) that are given as the functional derivative of the xc energy functional. Also in the KS equation, effects of the interaction between bosons in the superfluid and normal fluid components are included. Therefore, the hydrodynamic physical quantities of the superfluid and normal fluid components are interrelated, and they are determined using the GP-KS and KS equations of the sf-CDFT. This is a noteworthy feature of the sf-CDFT.

Conclusion
We developed the sf-CDFT, which applies to superfluids consisting of interacting bosons, such as a dilute Bose gas and liquid helium-4. In this theory, hydrodynamic physical quantities of the superfluid, i.e., the particle number density, current density, and order parameter of the BEC, can be reproduced by solving the GP-KS and KS equations in conjunction and self-consistently. The noteworthy feature of this theory is that the particle number densities and current densities of the superfluid and normal fluid components can be calculated simultaneously. In addition, the proposed theory can describe the interaction between the superfluid and normal fluid components. This means that the sf-CDFT may provide the microscopic basis of the two-fluid model.
It should be worthwhile to mention that the sf-CDFT applies to superfluid systems consisting of not only weakly interacting bosons but also strongly interacting bosons such as liquid helium-4.
In the proposed theory, a normal fluid component is described by the solutions of the KS equation. In contrast, the superfluid component is described by the solutions of the GP-KS equation. These equations are interrelated via the classical interaction and xc effective potentials given by the functional derivative of the xc energy functional. In other words, the normal fluid and superfluid components are in motion while mutually influencing each other. Thus, in this theory, the mutual friction between the normal fluid and superfluid components is incorporated in the form of classical interaction and xc effective potentials. Furthermore, the backflow effect [56][57][58] may be incorporated by devising the xc energy functional. Backflow occurs when the constituent particles of a fluid push other particles out of the way as they move. In other words, the backflow is caused by the interaction between bosons. Therefore, this effect can be incorporated into the xc energy functional, which is a topic for future studies.
Finally, we briefly comment on the actual calculation using the sf-CDFT. To perform actual calculations based on the sf-CDFT, the approximate form of the xc energy functional must be devised. As explained in section 3.1, we have derived a restrictive condition for the approximate form of the xc energy functional. In addition to developing the approximate form of the xc energy functional, we must develop a numerical calculation scheme for the differential equation that includes the vector potential term. We have already developed such a scheme for the case where a uniform magnetic field is applied [59,60]. With reference to this scheme, we may develop a numerical calculation scheme for the GP-KS and KS equations of the sf-CDFT.
where N is the total particle number operator. Subsequently, we can calculate y y j n r r r r , , , eq eq eq eq where X represents n, j , y, or y . † Therefore, the class of variables y y j n r r r r , , , eq eq eq eq is determined if the class of potentials is given.
Next, we prove that the class of potentials is uniquely determined when r eq is given. Assume that the density matrix is obtained from two different classes of potentials where Ĥ and ¢ Ĥ denote the Hamiltonian under two different classes of potentials respectively. Thus, equation (A3) can be rewritten as Because the right-hand side of equation (A4) is a c-number, -¢ H Ĥˆin the left-hand side should also be a c-number. This means that v r , respectively, by more than a constant. This contradicts the assumption. Thus, the class of potentials is determined uniquely when r eq is given. Consequently, the class of potentials and r eq have a one-to-one correspondence.
Next, we show that the class of variables and r eq are in one-to-one correspondence using Gibbs's variational principle [26]. According to Gibbs's variational principle, the functional of r¢, Equations (A6) and (A7) result in an inconsistent inequality: eq v D D eq eq eq eq eq eq eq eq eq eq eq eq eq eq eq eq eq eq eq eq eq eq , , , , where the universal functional is defined by equation (11)

Appendix C: Derivation of effective potentials
The universal functional defined in equation (11) can be rewritten as y y y y y y q y y = + - [ ] is the xc energy functional of the sf-CDFT for superfluids and is defined as the functional that contains the difference between the functionals of the real and reference systems.
Using equation (C4) and HK theorem II in the real system (equation (12)), we have Similarly, HK theorem II in the reference system (equation (18)) results in      S eq S eq S eq S eq { ( ) ( ) ( ) ( )}of the reference system. When equations (C6)-(C9) coincide with equations (C10)-(C13), respectively, the effective potentials included in equations (C10)-(C13) are considered potentials that can reproduce y y* j n r r r r , , , eq eq eq eq { ( ) ( ) ( ) ( )}in the reference system. Specifically, we can derive the following formulae for effective potentials, which reproduce the same basic variables as those of the real system: