Variational Approach to Many-Body Problems Incorporating Many-Body Effects at Finite Temperature

TL;DR

This paper introduces a finite-temperature variational method that self-consistently includes many-body correlations, revealing enhanced dynamical processes and potential qualitative changes in quasiparticle properties in weakly interacting Bose-Einstein condensates.

Contribution

A novel finite-temperature variational formalism that incorporates many-body effects through Green's functions and self-consistent parameters.

Findings

Many-body effects lower free energy compared to mean-field approximations.

Enhanced pair creation/annihilation processes at long wavelengths.

3/2-body correlations significantly affect quasiparticle self-energy.

Abstract

We develop a variational approach at finite temperature that incorporates many-body correlation self-consistently. The grand potential is constructed in terms of Green's function expressed by the variational parameters. We apply this formalism to weakly interacting Bose-Einstein condensates to incorporate the dynamical 3/2-body processes, which are considered important in the dynamical properties. The processes lower the free energy below the mean-field Hartree--Fock--Bogoliubov's value in the same way as a previous zero-temperature formalism. From our numerical results, the pair creation or annihilation processes neglected in the Popov--Shohno approximation are enhanced, particularly in the long wavelength region, owing to the many-body effects. Because the 3/2-body correlations give a finite contribution to the self-energy of quasiparticles, they may change the microscopic properties qualitatively, even in the weak-coupling region.

Equations

$${\Omega }_{\mathrm{v}}=⟨\hat{H}+{\beta }^{-1}\ln {\hat{\rho }}_{\mathrm{v}}⟩\ge \Omega ;\beta =\frac{1}{k_{B}T}$$

$$\hat{H}=\sum _{\mathbf{k}}({\epsilon }_k-\mu ){\hat{c}}_{\mathbf{k}}^{†}{\hat{c}}_{\mathbf{k}}+\frac{1}{2\mathcal{V}}\sum _{\mathbf{k},{\mathbf{k}}^{\prime },{\mathbf{q}}^{\prime }}{U}_q{\hat{c}}_{\mathbf{k}+\mathbf{q}}^{†}{\hat{c}}_{{\mathbf{k}}^{\prime }-\mathbf{q}}^{†}{\hat{c}}_{{\mathbf{k}}^{\prime }}{\hat{c}}_{\mathbf{k}},$$

$$\hat{H}={\hat{H}}_0+{\hat{H}}_1+{\hat{H}}_{3/2}+{\hat{H}}_2.$$

$${\hat{H}}_0\equiv$$

$${\hat{H}}_1\equiv$$

$$+\frac{1}{2\mathcal{V}}\text{sideset}{}^{\prime }{\sum }_{\mathbf{k}}{U}_k\left({\hat{c}}_{\mathbf{k}}^{†}{\hat{c}}_{-\mathbf{k}}^{†}{\hat{c}}_0{\hat{c}}_0+{\hat{c}}_0^{†}{\hat{c}}_0^{†}{\hat{c}}_{\mathbf{k}}{\hat{c}}_{-\mathbf{k}}\right),$$

$${\hat{H}}_{3/2}\equiv$$

$${\hat{H}}_2\equiv$$

$${\hat{\rho }}_{\mathrm{v}}=\exp \left[\beta \left({\Omega }_{\mathrm{v},\mathrm{LW}}-{\hat{H}}_{\mathrm{v}}\right)\right]$$

$${\Omega }_{\mathrm{v},\mathrm{LW}}=-\frac{1}{\beta }\ln \mathrm{Tr}{e}^{-\beta {\hat{H}}_{\mathrm{v}}}.$$

$${\Omega }_{\mathrm{v}}=\mathrm{Tr}{\hat{\rho }}_{\mathrm{v}}(\hat{H}-{\hat{H}}_{\mathrm{v}})+{\Omega }_{\mathrm{v},\mathrm{LW}}.$$

$${\hat{H}}_{\mathrm{v}}\equiv {\hat{H}}_{\mathrm{v},1}+{\hat{H}}_{\mathrm{v},3/2}.$$

$${\hat{H}}_{\mathrm{v},1}\equiv$$

$${\hat{H}}_{\mathrm{v},3/2}\equiv$$

$${\hat{\gamma }}_{\mathbf{k}}=u_{\mathbf{k}}{\hat{c}}_{\mathbf{k}}-v_{\mathbf{k}}{\hat{c}}_{-\mathbf{k}}^{†},$$

$$\left[\begin{array}[]{c}u_{\bm{k}}\\ v_{\bm{k}}\end{array}\right]\equiv\frac{1}{\sqrt{1-|\phi_{\bm{k}}|^{2}}}\left[\begin{array}[]{c}1\\ \phi_{\bm{k}}\end{array}\right],\ \ \ \phi_{\bm{k}}^{*}=\phi_{\bm{k}}=\phi_{-\bm{k}},$$

$${\hat{c}}_{\mathbf{k}}=u_{\mathbf{k}}{\hat{\gamma }}_{\mathbf{k}}+v_{\mathbf{k}}{\hat{\gamma }}_{-\mathbf{k}}^{†}.$$

$${\rho }_{\mathbf{k}}\equiv$$

$$=$$

$$F_{\mathbf{k}}\equiv$$

$$=$$

$$W_{{\mathbf{k}}_1{\mathbf{k}}_2;{\mathbf{k}}_3}\equiv$$

$$=$$

$${\mathcal{N}}_0\equiv \mathrm{Tr}{\hat{\rho }}_{\mathrm{v}}{\hat{c}}_0^{†}{\hat{c}}_0.$$

$$\mathrm{Tr}{\hat{\rho }}_{\mathrm{v}}{\hat{c}}_0^{†}{\hat{c}}_0^{†}{\hat{c}}_0{\hat{c}}_0\approx$$

$$\mathrm{Tr}{\hat{\rho }}_{\mathrm{v}}{\hat{c}}_0^{†}{\hat{c}}_0{\hat{c}}_{\mathbf{k}}^{†}{\hat{c}}_{\mathbf{k}}\approx$$

$$\mathrm{Tr}{\hat{\rho }}_{\mathrm{v}}{\hat{c}}_0^{†}{\hat{c}}_0^{†}{\hat{c}}_{\mathbf{k}}{\hat{c}}_{-\mathbf{k}}\approx$$

$$\mathrm{Tr}{\hat{\rho }}_{\mathrm{v}}{\hat{c}}_0^{†}{\hat{c}}_{-{\mathbf{k}}_3}^{†}{\hat{c}}_{{\mathbf{k}}_2}{\hat{c}}_{{\mathbf{k}}_1}\approx$$

$$\mathrm{Tr}{\hat{\rho }}_{\mathrm{v}}{\hat{c}}_{{\mathbf{k}}_1^{\prime }}^{†}{\hat{c}}_{{\mathbf{k}}_2^{\prime }}^{†}{\hat{c}}_{{\mathbf{k}}_2}{\hat{c}}_{{\mathbf{k}}_1}\approx$$

$$+{\delta }_{{\mathbf{k}}_1^{\prime },-{\mathbf{k}}_2^{\prime }}{\delta }_{{\mathbf{k}}_1,-{\mathbf{k}}_2}{F}_{{\mathbf{k}}_1^{\prime }}{F}_{{\mathbf{k}}_1}.$$

$${\Omega }_{\mathrm{v}}=$$

$$+\text{sideset}{}^{\prime }{\sum }_{\mathbf{k}}\left[{\epsilon }_k-\mu +{\bar{n}}_0(U_0+U_k)\right]{\rho }_{\mathbf{k}}+{\bar{n}}_0\text{sideset}{}^{\prime }{\sum }_{\mathbf{k}}{U}_k{F}_{\mathbf{k}}$$

$$+\frac{\sqrt{{\mathcal{N}}_0}}{\mathcal{V}}\text{sideset}{}^{\prime }{\sum }_{{\mathbf{k}}_1{\mathbf{k}}_2{\mathbf{k}}_3}{\delta }_{{\mathbf{k}}_1+{\mathbf{k}}_2+{\mathbf{k}}_3,0}{U}_{k_1}\left(W_{{\mathbf{k}}_1{\mathbf{k}}_2;{\mathbf{k}}_3}+W_{{\mathbf{k}}_1{\mathbf{k}}_2;{\mathbf{k}}_3}^{\ast }\right)$$

$$+\frac{1}{2\mathcal{V}}\text{sideset}{}^{\prime }{\sum }_{\mathbf{k}{\mathbf{k}}^{\prime }}\left[(U_0+U_{|\mathbf{k}-{\mathbf{k}}^{\prime }|}){\rho }_{\mathbf{k}}{\rho }_{{\mathbf{k}}^{\prime }}+U_{|\mathbf{k}-{\mathbf{k}}^{\prime }|}{F}_{\mathbf{k}}{F}_{{\mathbf{k}}^{\prime }}\right]$$

$$\displaystyle-\sideset{}{{}^{\prime}}{\sum}_{\bm{k}}E_{\bm{k}}\braket{\hat{\gamma}_{\bm{k}}^{\dagger}\hat{\gamma}_{\bm{k}}}-\frac{1}{3!}\sideset{}{{}^{\prime}}{\sum}_{\bm{k}_{1}\bm{k}_{2}\bm{k}_{3}}b_{\bm{k}_{1}\bm{k}_{2}\bm{k}_{3}}\biggl{(}\braket{\hat{\gamma}_{\bm{k}_{1}}^{\dagger}\hat{\gamma}_{\bm{k}_{2}}^{\dagger}\hat{\gamma}_{\bm{k}_{3}}^{\dagger}}$$

$$\displaystyle+\braket{\hat{\gamma}_{\bm{k}_{1}}\hat{\gamma}_{\bm{k}_{2}}\hat{\gamma}_{\bm{k}_{3}}}\biggr{)}+\Omega_{\rm{v,LW}},$$

$${\bar{n}}_0\equiv {\mathcal{N}}_0/\mathcal{V}.$$

$$⟨{\hat{\gamma }}_{\mathbf{k}}^{†}{\hat{\gamma }}_{\mathbf{k}}⟩=$$

$$⟨{\hat{\gamma }}_{{\mathbf{k}}_1}^{†}{\hat{\gamma }}_{{\mathbf{k}}_2}^{†}{\hat{\gamma }}_{{\mathbf{k}}_3}^{†}⟩=$$

$${\mathcal{G}}_{\mathbf{k}}({\tau }_1-{\tau }_2)\equiv -⟨{\hat{T}}_{\tau }{\hat{\gamma }}_{\mathbf{k}}({\tau }_1){\hat{\gamma }}_{\mathbf{k}}^{†}({\tau }_2)⟩,$$

$${\mathcal{G}}_{\mathbf{k}}(\tau )=\frac{1}{\beta }\sum _{n=-\infty }^{\infty }{e}^{-i{\epsilon }_n\tau }{\mathcal{G}}_{\mathbf{k}}(i{\epsilon }_n),{\epsilon }_n\equiv 2n\pi /\beta .$$

Full text

1 Introduction

The use of a variational approach is one of the most powerful methods for studying many-particle systems. In particular, some of the variational wave functions for the ground states have been constructed to consider correlations between particles or spins. [1, 2, 3, 4, 5] One of the most outstanding examples is the Gutzwiller wave function [1], which contains on-site correlations and has been used for the Hubbard model of electrons [6, 7] and bosons[8, 9, 10]. At finite temperatures, the variational principle is given in terms of the grand potential as follows:[11, 12]

[TABLE]

where $\Omega$ is the exact grand potential, and $\braket{\cdots}$ denotes a statistical average using the variational density matrix $\hat{\rho}_{\rm{v}}$. There also exist some variational approaches at finite temperatures for incorporating correlations, particularly in strongly correlated electron systems [13] and quantum spin systems [14] because it is difficult to describe these systems using the standard mean-field approaches, e.g., Hartree-Fock theory. On the other hand, most of the approaches adopt mean-field $\hat{\rho}_{\rm{v}}$.

The importance of the many-body effects is not limited to the strongly correlated system. Although they may contribute to the thermodynamic and dynamic properties of the system, they are often neglected because it is difficult to incorporate them. For example, in Bose-Einstein condensates (BECs), the dynamical 3/2-body processes, where one of two colliding non-condensed particles drops into the condensate and vice versa, are sources of many-body correlations. The importance of these processes is emphasized in the dynamics of BECs at finite temperatures or near the critical temperature.[15, 16] In recent years, however, we constructed a variational wave function for the ground state of weakly interacting bosons that self-consistently incorporates the dynamical 3/2-body processes.[17, 18, 19] Using it, we found that these processes contribute to the ground state energy with the same order as the mean-field energies. Thus, the 3/2-body correlations should be incorporated self-consistently, even in the collisionless regime.

In this study, we develop a variational approach at finite temperatures to investigate the many-body effects in equilibrium states, which does not depend on the specific form of Hamiltonians. We then apply it to the weak-coupling Bose-Einstein condensed phase and include the 3/2-body correlation self-consistently. Unlike the zero-temperature formalism[17, 18, 19], it is difficult to evaluate the expectation values of Bogoliubov’s quasiparticle operators.[20, 21] We evaluate them in terms of Green’s function of these operators expressed by the variational parameters.

The remainder of this paper is organized as follows. Section 2 constructs a variational density matrix with the 3/2-body correlation, and obtains an expression for the grand potential. Section 3 presents numerical results for the free energy, condensate fraction, and pair correlation. Section 4 summarizes this paper and describes some applications to other systems. Appendix discusses the connection with the zero-temperature formalism.[17] Hereafter, we adopt the units $\hbar=k_{B}=1$.

2 Formulation

2.1 Hamiltonian

We consider a system of $\mathcal{N}$ identical bosons with mass $m$ and spin 0 in a box of volume $\mathcal{V}$ described by the grand-canonical Hamiltonian

[TABLE]

where $\varepsilon_{k}\equiv k^{2}/2m$ is the kinetic energy, $\mu$ denotes the chemical potential, $(\hat{c}_{\bm{k}}^{\dagger},\hat{c}_{\bm{k}})$ are the field operators satisfying the Bose commutation relations, and $U_{q}$ is the interaction potential. We aim to describe the equilibrium states of Eq. (2) with Bose-Einstein condensation on the $\bm{k}=\bm{0}$ state. Hence, it is convenient to classify $\hat{H}$ according to the number of non-condensed operators involved as

[TABLE]

Each contribution on the right-hand side is given in terms of the primed sum $\displaystyle\sideset{}{{}^{\prime}}{\sum}_{\bm{k}}\equiv\sum_{\bm{k}}\left(1-\delta_{\bm{k},\bm{0}}\right)$ as

[TABLE]

2.2 Variational density matrix

We select $\hat{\rho}_{\rm{v}}$ in terms of a variational Hamiltonian $\hat{H}_{\rm{v}}$ as

[TABLE]

with

[TABLE]

Substituting Eq. (5) into Eq. (1), we obtain an expression for the variational grand potential as

[TABLE]

An appropriate choice of $\hat{H}_{v}$ is crucial for our variational study. Here, we express $\hat{H}_{\rm{v}}$ as the sum of two distinct contributions:

[TABLE]

The terms on the right side are defined as

[TABLE]

where $\hat{\gamma}_{\bm{k}}$ denotes the quasiparticle operator

[TABLE]

with

[TABLE]

so that $[\hat{\gamma}_{\bm{k}},\hat{\gamma}_{\bm{k}^{\prime}}^{\dagger}]=\delta_{\bm{k},\bm{k}^{\prime}}$ is satisfied. The inverse of Eq. (10) is given by

[TABLE]

By definition, the coefficient $b_{\bm{k}_{1}\bm{k}_{2}\bm{k}_{3}}$ is real and symmetric with respect to every permutation among $(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$. Using Eq. (12), one can evaluate various expectations of the operators in Eq. (2.3) in terms of Eq. (5) with Eq. (2.8). It is convenient for this purpose to introduce the following quantities for $\bm{k}\neq\bm{0}$:

[TABLE]

We define the number of condensed particles as

[TABLE]

Using this, we can express the expectations of the operator products in Eq. (2.3) as

[TABLE]

Because we consider the thermodynamic limit, we have approximated $(\hat{c}_{\bm{0}}^{\dagger})^{n}\hat{c}_{\bm{0}}^{m}$ as $\mathcal{N}_{0}^{(n+m)/2}$. In deriving Eqs. (15d) and (15e), we have also used the fact that the numbers of differences between the $\hat{\gamma}_{\bm{k}}^{\dagger}$ and $\hat{\gamma}_{\bm{k}}$ in any of their products should be multiples of 3 to give finite contributions in terms of Eq. (5) with Eq. (2.8). Using Eq. (2.14), we obtain an expression for Eq. (7) as

[TABLE]

with

[TABLE]

It follows that the two basic expectations in Eq. (16) are derivable from Eq. (6) as

[TABLE]

2.3 Expression of $\Omega_{\rm{v,LW}}$ and basic expectations

To obtain them, we introduce the quasiparticle Green’s function $\mathcal{G}_{\bm{k}}(\tau)$ as

[TABLE]

where $\hat{\gamma}_{\bm{k}}(\tau)\equiv e^{\tau\hat{H}_{v}}\hat{\gamma}_{\bm{k}}e^{-\tau\hat{H}_{v}},\ \hat{\gamma}_{\bm{k}}^{\dagger}(\tau)\equiv e^{\tau\hat{H}_{\rm{v}}}\hat{\gamma}_{\bm{k}}^{\dagger}e^{-\tau\hat{H}_{\rm{v}}}$, and $\hat{T}_{\tau}$ is a time-ordered operator for imaginary time $\tau$. It has the symmetry $\mathcal{G}_{\bm{k}}(\tau)=\mathcal{G}_{\bm{k}}(\tau+\beta)$, so we can expand it as

[TABLE]

Equation (6) can be regarded as the grand potential for the system of quasiparticles described by the effective Hamiltonian $\hat{H}_{\rm{v}}$. As shown by Luttinger and Ward,[22] it is expressible in terms of $\mathcal{G}_{\bm{k}}(i\varepsilon_{n})$ as

[TABLE]

where $p=(\bm{k},i\varepsilon_{n})$. $\Phi_{\rm{v,LW}}[\mathcal{G}]$ consists of all the skeleton diagrams in the simple perturbation expansion with respect to $\hat{H}_{{\rm{v}},3/2}$ for $\Omega_{\rm{v,LW}}$ with replacement of unperturbed Green’s functions $\mathcal{G}^{(0)}$ with $\mathcal{G}$, and the self-energy $\mathcal{S}_{\bm{k}}(i\varepsilon_{n})$ is derivable from it by

[TABLE]

Factor $e^{i\varepsilon_{n}0_{+}}$ can be omitted for the self-energy corresponding to the interaction of Eq. (9b). Eq. (21) has an important property in that it is stationary:

[TABLE]

with respect to the variation in $\mathcal{G}_{\bm{k}}(i\varepsilon_{n})$ that obeys Dyson’s equation:

[TABLE]

The key quantity in the functional of Eq. (21) is $\Phi_{\rm{v,LW}}$. We adopt the approximation of the lowest order for it as

[TABLE]

Substituting in Eq. (20), we can transform Eq. (25) into

[TABLE]

The corresponding self-energy is calculated by Eq. (22) to be

[TABLE]

From the stationary condition of Eq. (23), we can perform the differentiations of Eq. (2.17) with respect to the explicit dependences in Eq. (21) with Eq. (26) as

[TABLE]

2.4 Stationary conditions

Now, we derive the stationary conditions of Eq. (16) with respect to $\mathcal{N}_{0},\phi_{\bm{k}},E_{\bm{k}},$ and $b_{\bm{k}_{1}\bm{k}_{2}\bm{k}_{3}}$ along with the equation for the average particle number. Below, we set $b_{\bm{k}_{1}\bm{k}_{2}\bm{k}_{3}}=b_{-\bm{k}_{1}-\bm{k}_{2}-\bm{k}_{3}}$.

First, the equation $\mathcal{N}=-\partial\Omega_{\rm{v}}/\partial\mu$ for the average particle number gives an expression for the particle density $\bar{n}\equiv\mathcal{N}/\mathcal{V}$ as

[TABLE]

Next, the stationary condition $\delta\Omega_{\rm{v}}/\delta\mathcal{N}_{0}=0$ yields an expression for the chemical potential as

[TABLE]

To calculate stationary conditions with respect to $(\phi_{\bm{k}},E_{\bm{k}},b_{\bm{k}_{1}\bm{k}_{2}\bm{k}_{3}})$, it is convenient to introduce the following quantities:

[TABLE]

and

[TABLE]

With these preliminaries, the three conditions $\delta\Omega_{\rm{v}}/\delta\phi_{\bm{k}}=0,\ \delta\Omega_{\rm{v}}/\delta E_{\bm{k}}=0$, and $\delta\Omega_{\rm{v}}/\delta b_{\bm{k}_{1}\bm{k}_{2}\bm{k}_{3}}=0$ yield

[TABLE]

We note that Eq. (23) is equivalent to $\delta\Omega_{\rm{v}}/\delta\mathcal{G}_{\bm{k}}(i\varepsilon_{n})=0$ in our formalism. This means that $\Omega_{\rm{v}}$ is also the functional of $\mathcal{G}_{\bm{k}}(i\varepsilon_{n})$, and the dynamics of quasiparticles are optimized by the variational principle. While $\mathcal{G}_{\bm{k}}$ and $\mathcal{S}_{\bm{k}}$ are represented by $\phi_{\bm{k}},E_{\bm{k}}$, and $b_{\bm{k}_{1}\bm{k}_{2}\bm{k}_{3}}$, these parameters depend on Green’s functions through the basic expectations. Thus, we must determine $\phi_{\bm{k}},E_{\bm{k}},b_{\bm{k}_{1}\bm{k}_{2}\bm{k}_{3}}$, and $\mathcal{G}_{\bm{k}}$ self-consistently by solving Eqs. (33a)-(33c) with Eqs. (24), (27), (28a), (28b), and (2.28)-(2.31).

3 Numerical Results

3.1 Models and numerical procedures

Numerical calculations were performed for the contact interaction potential $U_{k}=U$.[23] For convenience, we alternatively express $U$ as

[TABLE]

A cutoff wavenumber $k_{\rm{c}}$ is introduced into every summation over $\bm{k}$ as

[TABLE]

to remove the ultraviolet divergence inherent in the potential. The $s$-wave scattering length $a$ of the potential is obtained by

[TABLE]

which yields

[TABLE]

We choose $k_{\rm{c}}$ such that $k_{\rm{c}}a_{U}\ll 1$ is satisfied, i.e., $a\approx a_{U}$.

The units of energy and wavenumber of this system are defined by

[TABLE]

respectively. They are used to transform Eq. (16) into the dimensionless form $\Omega_{\rm{v}}/\mathcal{N}\varepsilon_{U}$ for numerical calculations. Sums over $\bm{k}$ are transformed into integrals as follows [17]:

[TABLE]

where $\tilde{k}\equiv k/k_{U}$. $a_{U}^{3}\bar{n}$ gives the density within its interaction range. Here, we consider this coupling constant so that $a_{U}^{3}\bar{n}\ll 1$ is satisfied, corresponding to dilute ultracold atomic gases. The free energy $\mathcal{F}=\Omega+\mu\mathcal{N}$ is expressible as

[TABLE]

$c_{1}$ represents the contribution of $\braket{\hat{\gamma}_{\bm{k}}^{\dagger}\hat{\gamma}_{\bm{k}}}$ from $\hat{H}_{1}$, while $c_{2}$ is the contribution of the 3/2-body and 2-body correlations. The last term expresses an entropy defined by $S\equiv-\braket{\ln{\hat{\rho}_{v}}}=\beta(\braket{\hat{H}_{\rm{v}}}-\Omega_{\rm{v,LW}})$.

Numerical procedures are as follows. We use the result of the Bogoliubov approximation[20]

[TABLE]

and $b_{\bm{k}_{1}\bm{k}_{2}\bm{k}_{3}}=0$ as initial solutions for the self-consistent equations at finite temperatures.

We now make some comments on the summations over Matsubara frequency $\varepsilon_{n}$. First, some of them can be easily done computed using an imaginary time representation. For example,

[TABLE]

We use the fast Fourier transform (FFT)[24] to evaluate the Fourier transform of $\mathcal{G}(i\varepsilon_{n})$ into $\mathcal{G}(\tau)$. This can also be applied to the transform $\mathcal{S}(\tau)\to\mathcal{S}(i\varepsilon_{n})$ by discretizing $\tau$. To perform these calculations, we must choose a cutoff frequency $\varepsilon_{n_{\rm{c}}}$. We should note that the FFT of $\mathcal{G}_{\bm{k}}(i\varepsilon_{n})$ for $k\gg k_{U}$ may suffer from substantial numerical errors. However, it can be expected that the quasiparticles behave as the ideal bosons $E_{\bm{k}}\approx\varepsilon_{k}$ in the region $k\gg k_{U}$, so we can approximate $\mathcal{G}_{\bm{k}}(i\varepsilon_{n})\approx\mathcal{G}_{\bm{k}}^{(0)}(i\varepsilon_{n})\approx(i\varepsilon_{n}-\varepsilon_{k})^{-1}$ for $k\lesssim k_{\rm{c}}$. Because $|\mathrm{Re}\mathcal{G}^{(0)}_{\bm{k}}(i\varepsilon_{n})|\approx\varepsilon_{k}(\varepsilon_{n}^{2}+\varepsilon_{k}^{2})^{-1}$ and $|\mathrm{Im}\mathcal{G}_{\bm{k}}^{(0)}(i\varepsilon_{n})|\approx\varepsilon_{n}(\varepsilon_{n}^{2}+\varepsilon_{k}^{2})^{-1}$, it is sufficient to choose $\varepsilon_{n_{\rm{c}}}$ to satisfy $\varepsilon_{n_{\rm{c}}}\gg\varepsilon_{k_{\rm{c}}}$. We fixed the number of integration points for Matsubara frequency and imaginary time at $n_{\rm{c}}=2^{14}\sim 10^{4}$. Secondly, we must obtain the value of the Matsubara Green’s function at the discontinuity point $\tau=0_{\pm}$ from the FFT. It is known that the Fourier series $\tilde{f}_{\rm{FT}}(x)$ converges to $\frac{1}{2}\left[f(a+0_{+})+f(a+0_{-})\right]$ at the discontinuity point $x=a$. Because, $\mathcal{G}_{\bm{k}}(\tau=0_{+})-\mathcal{G}_{\bm{k}}(\tau=0_{-})=-\braket{\hat{\gamma}_{\bm{k}}\hat{\gamma}_{\bm{k}}^{\dagger}}+\braket{\hat{\gamma}_{\bm{k}}^{\dagger}\hat{\gamma}_{\bm{k}}}=-1$, $\braket{\hat{\gamma}_{\bm{k}}^{\dagger}\hat{\gamma}_{\bm{k}}}$ is expressed by

[TABLE]

Finally, because it is difficult to determine the imaginary time dependence, we evaluate the first term of Eq.(21) numerically by substituting $0_{+}=10^{-30}/\varepsilon_{n_{\rm{c}}}$. The convergence of the iteration can be checked by monitoring the grand potential energy $\Omega_{\rm{v}}$. We stopped the iteration when the magnitude of the relative difference between the old and new grand potentials decreased to below $10^{-10}$.

3.2 Results

In this section, we compare the results obtained by the effective Hamiltonian $\hat{H}_{\rm{v,1}}+\hat{H}_{\rm{v,3/2}}$ with those of $\hat{H}_{\rm{v}}=\hat{H}_{\rm{v,1}}$, i.e., the Hartree–Fock–Bogoliubov (HFB) approximation.[25] Because at the low temperatures, the 3/2-body correlations are considered to be the main sources of the many-body effects compared to the collisions between non-condensed particles, i.e., the dynamical 2-body correlations, we solved the self-consistent equations in the region $0<T<T_{\rm{c0}}/2$, where $T_{\rm{c0}}$ is the critical temperature for an ideal system.

First, we estimate the free energy per particle. We set the energy difference from the HFB theory as $\Delta\mathcal{F}\equiv\mathcal{F}-\mathcal{F}^{\rm{HFB}}$, where $\mathcal{F}^{\rm{HFB}}$ represents the free energy of the HFB approximation. Figure 1 shows the $T$ dependence of the deviations of the free energy $\Delta\mathcal{F}$ and the coefficient $\Delta c_{2}\equiv c_{2}-c_{2}^{\rm{HFB}}$. We observe that the free-energy is lower than $\mathcal{F}^{\rm{HFB}}$ ($\Delta\mathcal{F}<0$), and it decreases by an order of $a_{U}^{3}\bar{n}$. However, the deviations of $c_{1}$ and the entropy depend on the cutoff in our calculation. The inset in Fig. 1 shows that our calculation is numerically connected to the result of zero temperature formalism[17]. Figure 2 plots $\Delta\mathcal{F}$ as functions of $\log_{10}(a_{U}^{3}\bar{n})$. These results show that the free energy is reduced by enhancement of the many-body effect with density and temperature. In summary, the mean-field theory, which is the standard extension of Bogoliubov theory[20] to finite temperatures, is not effective for describing the weak-coupling BEC as well as the zero-temperature system, so the 3/2-body correlations should be incorporated self-consistently to investigate the thermodynamic properties, even in a collisionless regime.

To see the contributions of the many-body effects to the thermodynamic properties, we calculate the condensate fractions from Eq. (29). Figure 3 shows that the number of condensed particle decreases due to the 3/2-body correlations: $\Delta\mathcal{N}_{0}\equiv\mathcal{N}_{0}-\mathcal{N}_{0}^{\rm{HFB}}<0$. Because the 3/2-body correlations give the opportunity for particle exchange between condensate and non-condensate, they tend to enhance other correlations. Indeed, the pair correlations $F_{\bm{k}}=\mathrm{Tr}\hat{\rho}_{\rm{v}}\hat{c}_{\bm{k}}\hat{c}_{-\bm{k}}$ are enhanced due to the 3/2-body correlations, particularly in the long wavelength region, as shown in Fig. 4. This indicates that $F_{\bm{k}}$ should be incorporated self-consistently in BEC phase. Therefore, the Popov–Shohno approximation[27, 26], which is obtained by $\hat{H}_{\rm{v}}=\hat{H}_{\rm{v,1}}$ with $F_{\bm{k}}$ neglected to give a gapless excitation, is invalid even in the collisionless regime. This can also be confirmed from the fact that the free-energy of the Popov–Shohno approximation is higher than that of the HFB approximation; therefore, the pair creation or annihilation processes also contribute to the stabilization of the BEC system.

4 Summary

We have constructed a variational density matrix incorporating the many-body effects self-consistently. In general, if the trial density matrix includes the interaction terms, it is difficult to calculate the right hand side of Eq. (6) directly. The advantage of our formalism is that by considering $\Omega_{\rm{v,LW}}$ as the Luttinger–Ward functional and using Green’s function, it is possible to evaluate the basic expectations and derive Dyson’s equation that determines quasiparticle dynamics from the variational principle. We adopted this formalism in a BEC system, which incorporates the dynamical 3/2-body processes to give a lower free energy than the mean-field density matrix. Therefore, it turns out that the particle correlations contribute to stabilizing the system. Because the 3/2-body correlations provide the collisions between the condensate and non-condensate, they also contribute to increasing entropy. However, we could not confirm this through our numerical results.

A qualitative difference from the mean-field theory cannot be found in the above thermodynamic properties. On the other hand, the self-energy of the quasiparticles $\mathcal{S}_{\bm{k}}(i\varepsilon_{n})$ becomes finite due to the 3/2-body correlations. Since this quantity contributes to the spectrum of excitations in the form of an energy shift or width of the spectrum, qualitative changing of the single-particle excitation could occur, as discussed in the zero-temperature formalism.[17] Thus, the many-body effects may bring about qualitative changing of the microscopic properties from the mean-field theory, even in a collisionless regime. We intend to investigate this practically in the near future.

In this paper, we only considered the case of a single-component BEC phase at low temperatures. Our formalism can be used at higher temperatures, where the dynamical 2-body correlations are more important. In this case, only $\hat{H}_{\rm{v}}$ must be constructed so as to characterize these processes. It can also be adopted in other systems such as the superconductor and Bose-Fermi mixture systems as follows:

1. (i)

We construct $\hat{\rho}_{\rm{v}}$ that includes the terms considered to be the sources of the many-body effects in the systems. 2. (ii)

$\Omega_{\rm{v,LW}}$ is evaluated as a functional of Green’s function, and $\Phi_{\rm{v,LW}}$ determines the many-body effects. 3. (iii)

Because $\Omega_{\rm{v}}$ is expressible using the variational parameters, we determine them from Eq. (1).

It seems to be simple to extend to an inhomogeneous system such as the trapped system by using a general quantum number $q$ instead of the wavenumber $\bm{k}$ and Green’s function in coordinate space. Therefore, our variational approach can be used to study the many-body effects on the thermodynamic properties in other systems.

\acknowledgment\acknowledgment

W. K. is a JSPS Research Fellow, and this work was supported in part by KAKENHI Grand Number 18J13241. The numerical calculations were performed using the facilities of the Supercomputer Center, Institute for Solid State Physics, University of Tokyo. We would like to thank Editage (www.editage.jp) for English language editing.

Appendix A Connection to the Zero-Temperature Formalism

We adopt the Lehman representation for Green’s function:

[TABLE]

where $\mathcal{A}_{\bm{k}}(\varepsilon)$ satisfies the sum rule

[TABLE]

Using Eq. (43), we obtain an alternative expression for $\Phi_{\rm{v,LW}}[\mathcal{G}]$ and $\mathcal{S}_{\bm{k}}(i\varepsilon_{n})$ as

[TABLE]

[TABLE]

Here, $f(\varepsilon)$ denotes the Bose distribution function. It is convenient to express Eq. (46) in the Lehman representation too as

[TABLE]

$\varGamma_{\bm{k}}(\varepsilon)\equiv-2\rm{Im}\mathcal{S}_{\bm{k}}(\varepsilon+i0_{+})$ is given explicitly by

[TABLE]

The retarded self-energy is obtained from Eq. (47) by $\mathcal{S}_{\bm{k}}^{\rm{R}}(\varepsilon)=\mathcal{S}_{\bm{k}}(\varepsilon+i0_{+})$ as

[TABLE]

where $\rm{P}$ denotes the principal value. Using Eqs. (24) and (43), we obtain the spectral function $\mathcal{A}_{\bm{k}}(\varepsilon)=-2\rm{Im}\mathcal{G}_{\bm{k}}(\varepsilon+i0_{+})$ as

[TABLE]

The quantities in Eqs. (28a) and (28b) are expressible as

[TABLE]

Here, we show that our formalism connects with the zero-temperature formalism.[17] First, we substitute the zeroth-order of the spectral function

[TABLE]

into Eq. (48), where we approximate $E_{\bm{k}}\approx E_{\bm{k}}^{(0)}$with $E_{\bm{k}}^{(0)}$ given by Eq. (32a). Then, the first-order spectral function for the self-energy is expressible as

[TABLE]

The corresponding self-energy is given by Eq. (47) as

[TABLE]

Using Eq. (53) and neglecting $\rm{Re}\mathcal{S}^{\rm{R}}(\varepsilon)$ in Eq. (50), we obtain the first-order correction to the spectral function as

[TABLE]

where we have neglected $[\varGamma_{\bm{k}}^{(1)}(\varepsilon)/2]^{2}$ in the denominator. Substituting this expression into Eq. (51a), we obtain

[TABLE]

On the other hand, the zeroth-order expression of Eq. (52) may suffice to evaluate Eq. (51b) in the weak-coupling region, thus obtaining

[TABLE]

We clarify how the results of the zero-temperature formalism are reproduced.[17] Taking the limit of $T\rightarrow 0$ where $f(E_{\bm{k}})\rightarrow 0$ and substituting Eq. (57) with Eq.(33c) into Eq. (31c), we obtain the expression of $\chi_{\bm{k}}\equiv\lambda_{\bm{k}}u_{\bm{k}}^{2}$ that coincides with that of the zero-temperature formalism with the correspondence

[TABLE]

to the leading order. Eq. (56) in the limit of $T\rightarrow 0$ is identical to $\braket{\Phi}{\hat{\tilde{\gamma}}_{\bm{k}}^{\dagger}\hat{\tilde{\gamma}}_{\bm{k}}}{\Phi}$ in the zero-temperature formalism.

If we introduce the Lehmann representation for Green’s function in coordinate space, we can check the connection to the zero-temperature formalism in an inhomogeneous system[19] in the same manner. That is, it is considered that our formalism is a natural extension of the zero-temperature formalism[17, 18, 19]. On the other hand, we note that Eq. (56) cannot be used to solve the self-consistent equations because $\mathcal{A}_{\bm{k}}^{(1)}(\varepsilon)$ expressed by Eq. (55) does not satisfy the sum rule.