Stability of normal quantum-fluid mixtures

V.I. Yukalov1,2

1Bogolubov Laboratory of Theoretical Physics,
Joint Institute for Nuclear Research, Dubna 141980, Russia

2Instituto de Fisica de São Carlos, Universidade de São Paulo,
CP 369, São Carlos 13560-970, São Paulo, Brazil

E-mail<EMAIL_ADDRESS>

###### Abstract

Mixtures of quantum fluids, that is gases or liquids, are considered with the
emphasis on the conditions characterizing the stability of the mixtures. The
mixtures, that can be formed by cold atoms or molecules, are assumed to be
quantum requiring the description using quantum techniques, but not so cold
that to exhibit superfluidity or superconductivity. Reviewing the stability
conditions of such normal quantum systems is important for the comparison of
these conditions with the stability conditions of, e.g., Bose-condensed
mixtures. The behavior of observable quantities under the stratification of
quantum mixtures is discussed.

Keywords: quantum mixtures, normal quantum systems, stability conditions,
response functions, dynamic form-factors, structure factors, isothermic
compressibility, sum rules, role of temperature

##  1 Introduction

Quantum systems composed of atoms or molecules are nowadays intensively
studied due to the developed techniques of cooling and trapping [1]. The major
interest has been paid to the study of Bose systems in the condensed state,
which several books [2, 3, 4] and review articles [5, 6, 7, 8, 9, 10, 11, 12,
13, 14, 15, 16] are devoted to. Of particular interest are the mixtures of
quantum superfluids, Bose-Bose mixtures [17], Bose-Fermi mixtures [18], and
Fermi-Fermi superfluid mixtures [19].

One of the main points in considering quantum mixtures is the problem of their
stability with respect to their stratification into separate components
occupying different spatial locations. This type of stability, for brevity, is
called mixture stability. The mixture stability of superfluid Bose-Bose
mixtures is usually analyzed in the coherent approximation [20, 21] or
Bogolubov approximation [22]. The Hugenholtz-Pines relation [23] for
superfluid systems can be generalized for superfluid mixtures [24, 25]. The
stability of superfluid mixture in the Hartree-Fock-Bogolubov approximation
has been recently analyzed [26].

The mixture stability for normal (not superfluid) quantum fluids has been
considered to a much lesser extent. The aim of the present paper is to
compensate this deficiency and to survey what is known on the stability
conditions of normal quantum mixtures.

The principal difference in the study of mixture stability for superfluid and
normal fluids is in the following. The single-particle spectrum of a
superfluid, where the gauge symmetry is broken, coincides with the spectrum of
collective excitations [27, 28]. Therefore, to investigate the dynamic
stability of a superfluid spectrum it is sufficient to consider the poles of
the single-particle Green function. Contrary to this, for a normal quantum
system, where the gauge symmetry is preserved, it is necessary to study the
poles of the two-particle Green function or of the response function (dynamic
susceptibility). For a normal system, the single-particle Green function does
not define the spectrum of collective excitations.

It is interesting and important to compare the conditions of mixture stability
for superfluid and normal quantum systems. Also, it is useful to compare
dynamic and thermodynamic stability [29, 30]. Dynamic stability requires that
the spectrum of collective excitations be positive for finite momenta. And
thermodynamic stability studies conditions where the thermodynamic potential
of a mixture is minimal.

We consider the systems composed of neutral atoms or molecules, whose density
fluctuations enjoy sound modes. The particles can be bosons or fermions,
provided that sound modes are well defined for them. The instability of a
mixture with respect to the spatial separation of the components can be caused
by the interspecies interactions. The instability can arise when the repulsive
interspecies interactions become sufficiently strong. The derivation of
conditions characterizing the mixture stability of normal quantum-fluid
mixtures is the main topic of the present survey.

Throughout the paper, the system of units is employed where the Planck and
Boltzmann constants are set to one.

##  2 Quantum-fluid mixture

Let us start with the analysis of dynamic stability of a quantum-fluid
mixture. We consider particles without internal degrees of freedom, such as
spin. In the case of particles with spin, this could imply that their spins
are frozen in one direction and thus are switched off the particle dynamics.

The grand Hamiltonian of the mixture of several components, enumerated by the
index $i$, has the form

$H=\sum_{i}\int\psi_{i}^{\dagger}({\bf
r})\;\left(-\;\frac{\nabla^{2}}{2m}+U_{i}-\mu_{i}\right)\;\psi_{i}({\bf
r})\;d{\bf r}\;+$ $+\;\frac{1}{2}\sum_{ij}\psi_{i}^{\dagger}({\bf
r})\;\psi_{j}^{\dagger}({\bf r}^{\prime})\;\Phi_{ij}({\bf r}-{\bf
r}^{\prime})\;\psi_{j}({\bf r}^{\prime})\;\psi_{i}({\bf r})\;d{\bf r}d{\bf
r}^{\prime}\;,$ (1)

where $\psi_{i}({\bf r})=\psi_{i}({\bf r},t)$ are field operators (either Bose
or Fermi), $U_{i}=U_{i}({\bf r},t)$ is an external potential, $\mu_{i}$ are
chemical potentials, and $\Phi_{ij}({\bf r})$ is an interaction potential.

The response functions are the variational derivatives

$\chi_{ij}({\bf r},t,{\bf r}^{\prime},t^{\prime})=\frac{\delta\rho_{i}({\bf
r},t)}{\delta U_{j}({\bf r}^{\prime},t^{\prime})}$ (2)

of the density

$\rho_{i}({\bf r},t)=\langle\;\psi_{i}^{\dagger}({\bf r},t)\;\psi_{i}({\bf
r},t)\;\rangle=\pm iG_{i}({\bf r},t,{\bf r},t)\;,$ (3)

where we assume the causal Green functions

$G_{i}({\bf r},t,{\bf
r}^{\prime},t^{\prime})=-i\;\langle\;\hat{T}\psi_{i}({\bf
r},t)\;\psi_{i}^{\dagger}({\bf r}^{\prime},t^{\prime})\;\rangle\;.$ (4)

The Green functions at coinciding times are defined as

$G_{i}({\bf r},t,{\bf r}^{\prime},t)=\lim_{t^{\prime}\rightarrow
t+0}\;G_{i}({\bf r},t,{\bf r}^{\prime},t^{\prime})\;.$ (5)

Some interaction potentials are known to be non-integrable and not allowing
for their Fourier transformation. This obstacle can be overpassed by accepting
an effective interaction potential smoothed by the pair correlation function
[31]. It has been shown [7, 32, 33] that, starting from such a correlated
approximation, it is possible to formulate a self-consistent iterative theory
of all orders. Keeping this in mind, we assume that the interaction potential
is integrable and allows for the Fourier transformation

$\Phi_{ij}({\bf r})=\int\overline{\Phi}_{ij}({\bf k})\;e^{i{\bf k}\cdot{\bf
r}}\;\frac{d{\bf k}}{(2\pi)^{3}}\;,\qquad\overline{\Phi}_{ij}({\bf
k})=\int\Phi_{ij}({\bf r})\;e^{-i{\bf k}\cdot{\bf r}}\;d{\bf r}\;,$ (6)

with the Fourier transform (6) existing for all $k$, including $k=0$.

Considering a uniform system in volume $V$, we have the constant densities of
the components

$\rho_{j}\equiv\frac{N_{j}}{V}\qquad\left(N=\sum_{j}N_{j}\right)\;,$ (7)

hence the constant total density

$\rho\equiv\frac{N}{V}=\sum_{j}\rho_{j}\;.$ (8)

Assuming that the single-particle spectrum is real, for the Fourier transform
of the Green function (4) we have

$G_{i}({\bf k},\omega)=\frac{1\pm n_{i}({\bf k})}{\omega-\omega_{i}({\bf
k})+i0}\;\mp\;\frac{n_{i}({\bf k})}{\omega-\omega_{i}({\bf k})-i0}\;,$ (9)

with the single-particle spectrum given by the equation

$\omega_{i}({\bf k})=\frac{k^{2}}{2m}+\Sigma_{i}({\bf
k},\omega_{i})-\mu_{i}\;,$ (10)

where $\Sigma_{i}({\bf k},\omega)$ is a self-energy and the momentum
distribution is

$n_{i}({\bf k})=\frac{1}{\exp\\{\beta\omega_{i}({\bf k})\\}\mp 1}\qquad(\beta
T=1)\;,$ (11)

with $T$ being temperature. Here the upper sign is for Bose-Einstein
statistics, while the lower sign, for Fermi-Dirac statistics. The chemical
potentials are defined by the normalization to the density

$\rho_{i}=\int n_{i}({\bf k})\;\frac{d{\bf k}}{(2\pi)^{3}}\;.$ (12)

##  3 Random-phase approximation

The response functions are usually calculated in the random-phase
approximation [34] with the self-energy taken in the Hartree form

$\Sigma_{i}({\bf k},\omega)=\sum_{j}\Phi_{ij}\;\rho_{j}\;,$ (13)

where

$\Phi_{ij}\equiv\int\Phi_{ij}({\bf r})\;d{\bf r}=\overline{\Phi}_{ij}(0)\;.$
(14)

Then the Fourier transforms of the mixture response functions satisfy the
equations

$\chi_{ij}({\bf k},\omega)=\Pi_{ij}({\bf k},\omega)+\Pi_{ii}({\bf
k},\omega)\sum_{n}\overline{\Phi}_{in}({\bf k})\;\chi_{nj}({\bf k},\omega)\;,$
(15)

with the polarization functions

$\Pi_{ij}({\bf k},\omega)=\pm i\delta_{ij}\int G_{i}({\bf k}+{\bf
k}^{\prime},\omega+\omega^{\prime})\;G_{i}({\bf
k}^{\prime},\omega^{\prime})\;\frac{d{\bf
k}^{\prime}d\omega^{\prime}}{(2\pi)^{4}}\;.$ (16)

For the Green function (9), we get

$\Pi_{ii}({\bf k},\omega)=\frac{k^{2}}{m_{i}}\int\frac{n_{i}({\bf
k}^{\prime})}{(\omega-{\bf k}\cdot{\bf
k}^{\prime}/m_{i})^{2}-(k^{2}/2m_{i})^{2}}\;\frac{d{\bf
k}^{\prime}}{(2\pi)^{3}}\;.$ (17)

In the case of two components, the solutions to these equations are

$\chi_{11}=\frac{\Pi_{11}(1-\Pi_{22}\overline{\Phi}_{22})}{(1-\Pi_{11}\overline{\Phi}_{11})(1-\Pi_{22}\overline{\Phi}_{22})-\Pi_{11}\;\Pi_{22}\overline{\Phi}_{12}^{2}}\;,$
$\chi_{12}=\frac{\Pi_{11}\Pi_{22}\overline{\Phi}_{12}}{(1-\Pi_{11}\overline{\Phi}_{11})(1-\Pi_{22}\overline{\Phi}_{22})-\Pi_{11}\;\Pi_{22}\overline{\Phi}_{12}^{2}}\;.$
(18)

Here, for compactness, the variables ${\bf k}$ and $\omega$ are omitted. The
response functions $\chi_{21}$ and $\chi_{22}$ can be obtained by exchanging
the indices $1$ and $2$ in the above expressions.

We consider here an arbitrary potential $\Phi({\bf r})$. In particular, for
dilute systems of trapped atoms, one has

$\Phi_{ii}=4\pi\;\frac{a_{i}}{m_{i}}\;,\qquad\Phi_{ij}=2\pi\;\frac{a_{ij}}{m_{ij}}\;\qquad\left(m_{ij}\equiv\frac{m_{i}m_{j}}{m_{i}+m_{j}}\right)\;,$

where $a_{i}$ and $a_{ij}$ are scattering lengths.

##  4 Excitation spectrum

The spectrum of collective excitations is given by the poles of the response
functions, that is by the equation

$\chi_{ij}^{-1}({\bf k},\omega)=0\;.$ (19)

It is convenient to introduce the function

$f_{i}({\bf
k},\omega)\equiv\frac{1}{\rho_{i}}\int\frac{[\omega^{2}-(k^{2}/2m_{i})^{2}]n_{i}({\bf
k}^{\prime})}{(\omega-{\bf k}\cdot{\bf
k}^{\prime}/m_{i})^{2}-(k^{2}/2m_{i})^{2}}\;\frac{d{\bf
k}^{\prime}}{(2\pi)^{3}}$ (20)

satisfying the limits

$\lim_{k\rightarrow 0}f_{i}({\bf
k},\omega)=\lim_{k\rightarrow\infty}f_{i}({\bf k},\omega)=1\;.$ (21)

For the mixture of two components, the spectrum of collective excitations can
be represented [35] as

$\varepsilon^{2}_{\pm}=\frac{1}{2}\;\left[\;\varepsilon_{1}^{2}+\varepsilon_{2}^{2}\;\pm\;\sqrt{(\varepsilon_{1}^{2}-\varepsilon_{2}^{2})^{2}+4\varepsilon_{12}^{4}}\;\right]\;,$
(22)

where

$\varepsilon_{i}^{2}({\bf k})=s_{i}^{2}({\bf
k})\;k^{2}+\left(\frac{k^{2}}{2m_{i}}\right)^{2}\;,$ (23)

with the effective sound velocity

$s_{i}^{2}({\bf k})=\frac{\rho_{i}}{m_{i}}\;\overline{\Phi}_{ii}({\bf
k})f_{i}({\bf k},\varepsilon_{\pm})\;,$ (24)

and where

$\varepsilon_{12}^{2}({\bf k})=s_{12}^{2}({\bf k})\;k^{2}\;,$ (25)

with

$s_{12}^{2}({\bf k})=\sqrt{\frac{\rho_{1}\rho_{2}}{m_{1}m_{2}}\;f_{1}({\bf
k},\varepsilon_{\pm})\;f_{2}({\bf
k},\varepsilon_{\pm})}\;\overline{\Phi}_{12}({\bf k})\;.$ (26)

In the long-wave limit, the spectra of a binary mixture reduce to the sound-
type form

$\varepsilon_{\pm}({\bf k})\simeq s_{\pm}\;k\qquad(k\rightarrow 0)\;,$ (27)

with the sound velocities given by the equation

$s_{\pm}^{2}=\frac{1}{2}\;\left[\;s_{1}^{2}+s_{2}^{2}\;\pm\;\sqrt{(s_{1}^{2}-s_{2}^{2})^{2}+4s_{12}^{4}}\;\right]\;,$
(28)

where

$s_{i}^{2}\equiv s_{i}^{2}(0)=\frac{\rho_{i}}{m_{i}}\;\Phi_{ii}\;,\qquad
s_{12}^{2}\equiv
s_{12}^{2}(0)=\sqrt{\frac{\rho_{1}\rho_{2}}{m_{1}m_{2}}}\;\Phi_{12}\;.$ (29)

From here, it is seen that for the separate stability of each of the
components, it should be

$\Phi_{ii}>0\;.$ (30)

However, the interaction of different components, $\Phi_{12}$ can be of both
signs.

In the short-wave limit, one has

$\varepsilon_{\pm}({\bf
k})\simeq\frac{k^{2}}{2m_{i}}\qquad(k\rightarrow\infty)\;.$ (31)

As is clear, the spectrum branch $\varepsilon_{+}({\bf k})$, describing the
total density oscillations, is non-negative, and the branch
$\varepsilon_{-}({\bf k})$, characterizing the oscillations of the components
with respect to each other, is non-negative under the condition

$\varepsilon_{12}^{2}({\bf k})<\varepsilon_{1}({\bf k})\;\varepsilon_{2}({\bf
k})\;,$ (32)

which is the condition of mixture stability. In the long-wave limit, we have

$s_{12}^{2}<s_{1}\;s_{2}\;,$ (33)

which gives the stability condition

$\Phi_{12}^{2}<\Phi_{11}\;\Phi_{22}\;.$ (34)

It is interesting that this is the same condition as for a coherent,
completely Bose-condensed mixture [20]. In general, depending on the
interaction potentials, there can appear a finite momentum where the stability
condition (32) becomes broken, as under the roton instability in dipolar
systems [36].

##  5 Moving fluids

The stability conditions also depend on the motion of fluids [20, 30, 35]. Let
us consider a binary mixture, such that the components move with velocities
${\bf v}_{i}$. Then, acting similarly to the above analysis, we come to the
equation for the spectrum of collective excitations

$\left[\;(\omega-\alpha_{1})^{2}-\varepsilon_{1}^{2}\;\right]\left[\;(\omega-\alpha_{2})^{2}-\varepsilon_{2}^{2}\;\right]=\varepsilon_{12}^{4}\;,$
(35)

in which

$\alpha_{i}({\bf k})\equiv{\bf v}_{i}\cdot{\bf
k}=v_{i}\;k\;\cos(\vartheta_{i})\;.$ (36)

By the Descartes polynomial theorem [37] this equation enjoys two real roots
provided that

$\varepsilon_{12}^{4}({\bf k})<\left[\;\varepsilon_{1}^{2}({\bf
k})-\alpha_{1}^{2}({\bf k})\;\right]\left[\;\varepsilon_{2}^{2}({\bf
k})-\alpha_{2}^{2}({\bf k})\;\right]$ (37)

for all angles $\vartheta_{i}$. Minimizing the right-hand side of the above
inequality yields

$\varepsilon_{12}^{4}({\bf k})<\left[\;\varepsilon_{1}^{2}({\bf
k})-v_{1}^{2}\;k^{2}\;\right]\left[\;\varepsilon_{2}^{2}({\bf
k})-v_{2}^{2}\;k^{2}\;\right]$ (38)

for the angles $\vartheta_{i}=0,\pi$. For the broken condition (38), under the
angles $\vartheta_{1}=0$ and $\vartheta_{2}=\pi$ or $\vartheta_{1}=\pi$ and
$\vartheta_{2}=0$, this is called the counter-flow instability, and when
$\vartheta_{1}=\vartheta_{2}$, this is the instability of motion.

In the long-wave limit, inequality (38) reduces to the form

$s_{12}^{4}<(s_{1}^{2}-v_{1}^{2})(s_{2}^{2}-v_{2}^{2})\;.$ (39)

In terms of the interaction potentials, this is equivalent to the condition

$\Phi_{12}^{2}<\left(\Phi_{11}\;-\;\frac{m_{1}v_{1}^{2}}{\rho_{1}}\right)\left(\Phi_{22}\;-\;\frac{m_{2}v_{2}^{2}}{\rho_{2}}\right)\;.$
(40)

##  6 Dynamic form-factor

The measured quantity under neutron or x-ray scattering is the double-
differential cross section (see [5, 38])

$\frac{1}{N}\;\frac{d^{2}\sigma}{d\Omega_{k}\;d\omega}=b_{sc}^{2}\;\frac{k_{sc}}{k_{in}}\;S_{tot}({\bf
k},\omega)\;,$ (41)

in which the transferred momentum and energy are

${\bf k}\equiv{\bf k}_{in}-{\bf k}_{sc}\;,\qquad\omega\equiv E_{in}-E_{sc}\;,$

$\Omega_{k}$ is a solid angle around ${\bf k}$, and $b_{sc}$ is the scattering
length of a neutron on an atom. The dynamic form-factor at equilibrium is the
Fourier transform

$S_{tot}({\bf k},\omega)=\frac{1}{N}\sum_{ij}\int D_{ij}({\bf r},t,{\bf
r}^{\prime},0)\;e^{-i{\bf k}\cdot({\bf r}-{\bf r}^{\prime})+i\omega t}\;d{\bf
r}d{\bf r}^{\prime}dt$ (42)

of the density-density correlation function

$D_{ij}({\bf r},t,{\bf r}^{\prime},t^{\prime})\equiv\langle\;\hat{n}_{i}({\bf
r},t)\;\hat{n}_{j}({\bf r}^{\prime},t^{\prime})\;\rangle\;,$ (43)

where the density operator is

$\hat{n}_{i}({\bf r},t)\equiv\psi_{i}^{\dagger}({\bf r},t)\;\psi_{i}({\bf
r},t)\;.$ (44)

The total dynamic form-factor is the sum of two terms,

$S_{tot}({\bf k},\omega)=S_{el}({\bf k},\omega)+S({\bf k},\omega)\;.$ (45)

The first term

$S_{el}({\bf k},\omega)=\frac{1}{N}\sum_{ij}\int\rho_{i}({\bf
r},t)\;\rho_{j}({\bf r}^{\prime},0)\;e^{-i{\bf k}\cdot({\bf r}-{\bf
r}^{\prime})+i\omega t}\;d{\bf r}d{\bf r}^{\prime}dt\;,$ (46)

where

$\rho_{i}({\bf r},t)\equiv\langle\;\hat{n}_{i}({\bf r},t)\;\rangle\;,$ (47)

describes elastic scattering. The second term

$S({\bf k},\omega)=\frac{1}{N}\sum_{ij}\int\left[\;D_{ij}({\bf r},t,{\bf
r}^{\prime},0)-\rho_{i}({\bf r},t)\;\rho_{j}({\bf
r}^{\prime},0)\;\right]\;e^{-i{\bf k}\cdot({\bf r}-{\bf r}^{\prime})+i\omega
t}\;d{\bf r}d{\bf r}^{\prime}dt\;,$ (48)

also called the van Hove function, describes inelastic scattering.

For a uniform system, where

$\rho_{i}({\bf r},t)=\rho_{i}=\frac{N_{i}}{V}\;,$

the elastic scattering term is

$S_{el}({\bf
k},\omega)=(2\pi)^{4}\sum_{ij}\frac{\rho_{i}\rho_{j}}{\rho}\;\delta({\bf
k})\;\delta(\omega)\;.$ (49)

The inelastic term takes the form

$S({\bf k},\omega)=\sum_{ij}S_{ij}({\bf k},\omega)\;,$ (50)

where

$S_{ij}({\bf k},\omega)=\frac{1}{N}\int\left[\;D_{ij}({\bf r},t,{\bf
r}^{\prime},0)-\rho_{i}({\bf r},t)\;\rho_{j}({\bf
r}^{\prime},0)\;\right]\;e^{-i{\bf k}\cdot({\bf r}-{\bf r}^{\prime})+i\omega
t}\;d{\bf r}d{\bf r}^{\prime}dt\;.$ (51)

Defining the Fourier transform of the density-density correlation function

$D_{ij}({\bf k},\omega)=\int D_{ij}({\bf r},t,0,0)\;e^{-i{\bf k}\cdot{\bf
r}+i\omega t}\;d{\bf r}dt$ (52)

gives for a uniform system

$S_{ij}({\bf k},\omega)=\frac{1}{\rho}\;\left[\;D_{ij}({\bf
k},\omega)-(2\pi)^{4}\;\rho_{i}\;\rho_{j}\;\delta({\bf
k})\;\delta(\omega)\;\right]\;.$ (53)

The density-density correlation function (43) is connected with the density-
density Green function

$-i\;\langle\;\hat{T}\;\hat{n}({\bf r},t)\;\hat{n}({\bf
r},t^{\prime})\;\rangle\;.$

Using the properties of this Green function [39] results in the Fourier
transform

$D_{ij}({\bf k},\omega)=-\;\frac{2{\rm Im}\chi_{ij}({\bf
k},\omega)}{1+e^{-\beta\omega}}+(2\pi)^{4}\;\rho_{i}\;\rho_{j}\;\delta({\bf
k})\;\delta(\omega)\;,$ (54)

where $\chi_{ij}({\bf k},\omega)$ is the Fourier transform of the response
function (2). Respectively, for the dynamic form-factor (51), we get

$S_{ij}({\bf k},\omega)=-\;\frac{2{\rm Im}\chi_{ij}({\bf
k},\omega)}{\rho(1+e^{-\beta\omega})}\;.$ (55)

The response function enjoys the spectral representation

$\chi_{ij}({\bf k},\omega)=\int_{-\infty}^{\infty}\frac{\gamma_{ij}({\bf
k},\omega^{\prime})}{\omega-\omega^{\prime}}\;\frac{d\omega^{\prime}}{2\pi}\;,$
(56)

with the spectral function

$\gamma_{ij}({\bf k},\omega)=i\;[\;\chi_{ij}({\bf k},\omega+i0)-\chi_{ij}({\bf
k},\omega-i0)\;]\;,$ (57)

and the dispersion relation

${\rm Im}\chi_{ij}({\bf k},\omega)=-\;\frac{1}{2}\;\gamma_{ij}({\bf
k},\omega)\;\coth\left(\frac{\beta\omega}{2}\right)\;.$ (58)

This gives

$S_{ij}({\bf k},\omega)=\frac{\gamma_{ij}({\bf
k},\omega)}{\rho(1-e^{-\beta\omega})}\;.$ (59)

The dynamic form-factor possesses the properties

$S_{ij}(-{\bf k},\omega)=S_{ij}({\bf k},\omega)\;,\qquad S_{ij}({\bf
k},-\omega)=S_{ij}({\bf k},\omega)\;e^{-\beta\omega}\;.$ (60)

It can be written as the sum

$S_{ij}({\bf k},\omega)=S_{ij}^{+}({\bf k},\omega)+S_{ij}^{-}({\bf
k},-\omega)$ (61)

of the symmetric term

$S_{ij}^{+}({\bf k},\omega)\equiv\frac{1}{2}\;[\;S_{ij}({\bf
k},\omega)+S_{ij}({\bf k},-\omega)\;]=\frac{1}{2\rho}\;\gamma_{ij}({\bf
k},\omega)\;\coth\left(\frac{\beta\omega}{2}\right)=S_{ij}^{+}({\bf
k},-\omega)$ (62)

and the antisymmetric term

$S_{ij}^{-}({\bf k},\omega)\equiv\frac{1}{2}\;[\;S_{ij}({\bf
k},\omega)-S_{ij}({\bf k},-\omega)\;]=\frac{1}{2\rho}\;\gamma_{ij}({\bf
k},\omega)=-S_{ij}^{-}({\bf k},-\omega)$ (63)

with respect to the change of the sign of $\omega$.

##  7 Response functions

The response functions (18) can be written in the form

$\chi_{11}({\bf
k},\omega)=\left(\frac{\rho_{1}}{m_{1}}\right)\;\frac{[\omega^{2}-\varepsilon_{2}^{2}({\bf
k})]k^{2}}{[\omega^{2}-\varepsilon_{+}^{2}({\bf
k})][\omega^{2}-\varepsilon_{-}^{2}({\bf k})]}\;,$ $\chi_{12}({\bf
k},\omega)=\sqrt{\frac{\rho_{1}\rho_{2}}{m_{1}m_{2}}}\;\frac{\varepsilon_{12}^{2}({\bf
k})k^{2}}{[\omega^{2}-\varepsilon_{+}^{2}({\bf
k})][\omega^{2}-\varepsilon_{-}^{2}({\bf k})]}\;.$ (64)

Then the spectral functions (57) become

$\gamma_{11}({\bf
k},\omega)=\frac{\rho_{1}}{m_{1}}\;[\;\omega^{2}-\varepsilon_{2}^{2}({\bf
k})\;]\;k^{2}\;\gamma({\bf k},\omega)\;,$ $\gamma_{12}({\bf
k},\omega)=\sqrt{\frac{\rho_{1}\rho_{2}}{m_{1}m_{2}}}\;\varepsilon_{12}^{2}({\bf
k})\;k^{2}\;\gamma({\bf k},\omega)\;,$ (65)

where

$\gamma({\bf
k},\omega)=2\pi\delta\left(\;[\;\omega^{2}-\varepsilon_{+}^{2}({\bf
k})\;]\;[\;\omega^{2}-\varepsilon_{-}^{2}({\bf k})\;]\;\right)=$
$=\frac{\pi}{\varepsilon_{+}^{2}({\bf k})-\varepsilon_{-}^{2}({\bf
k})}\;\left\\{\frac{1}{\varepsilon_{+}({\bf
k})}\;[\;\delta(\omega-\varepsilon_{+}({\bf
k}))-\delta(\omega+\varepsilon_{+}({\bf k}))\;]\;-\right.$
$-\;\left.\frac{1}{\varepsilon_{-}({\bf
k})}\;[\;\delta(\omega-\varepsilon_{-}({\bf
k}))-\delta(\omega+\varepsilon_{-}({\bf k}))\;]\;\right\\}\;.$ (66)

##  8 Structure factor

The system structure factor

$S({\bf k})=\sum_{ij}S_{ij}({\bf k})$ (67)

is the sum of the partial structure factors

$S_{ij}({\bf k})=\int_{-\infty}^{\infty}S_{ij}({\bf
k},\omega)\;\frac{d\omega}{2\pi}\;.$ (68)

For the latter, we have

$S_{ij}({\bf k})=\frac{1}{\rho}\int_{0}^{\infty}\gamma_{ij}({\bf
k},\omega)\coth\left(\frac{\beta\omega}{2}\right)\;d\omega\;.$ (69)

Explicitly, for a binary mixture, this reads as

$S_{11}({\bf
k})=\left(\frac{\rho_{1}}{m_{1}}\right)\;\frac{k^{2}}{2\rho[\;\varepsilon_{+}^{2}({\bf
k})-\varepsilon_{-}^{2}({\bf k})\;]}\;\times$
$\times\left[\;\frac{\varepsilon_{+}^{2}({\bf k})-\varepsilon_{2}^{2}({\bf
k})}{\varepsilon_{+}({\bf k})}\;\coth\left(\frac{\beta\varepsilon_{+}({\bf
k})}{2}\right)\;-\;\frac{\varepsilon_{-}^{2}({\bf k})-\varepsilon_{2}^{2}({\bf
k})}{\varepsilon_{-}({\bf k})}\;\coth\left(\frac{\beta\varepsilon_{-}({\bf
k})}{2}\right)\;\right]$ (70)

and

$S_{12}({\bf
k})=\sqrt{\frac{\rho_{1}\rho_{2}}{m_{1}m_{2}}}\;\frac{\varepsilon_{12}^{2}({\bf
k})k^{2}}{2\rho[\varepsilon_{+}^{2}({\bf k})-\varepsilon_{-}^{2}({\bf
k})]}\;\times$ $\times\left[\;\frac{1}{\varepsilon_{+}({\bf
k})}\;\coth\left(\frac{\beta\varepsilon_{+}({\bf
k})}{2}\right)\;-\;\frac{1}{\varepsilon_{-}({\bf
k})}\;\coth\left(\frac{\beta\varepsilon_{-}({\bf k})}{2}\right)\;\right]\;.$
(71)

It is possible to check that in the short-wave limit

$\lim_{k\rightarrow\infty}\;\sum_{ij}S_{ij}({\bf k})=1\;.$ (72)

The limiting behavior of the partial structure factors depends on the relation
between the sound modes and temperature. If the sound-mode energy is much
smaller than temperature, then the structure factors tend to the limits

$S_{11}(0)=\left(\frac{\rho_{1}}{m_{1}}\right)\;\frac{s_{2}^{2}(0)T}{\rho
s_{+}^{2}(0)s_{-}^{2}(0)}\;,$
$S_{12}(0)=\sqrt{\frac{\rho_{1}\rho_{2}}{m_{1}m_{2}}}\;\frac{s_{12}^{2}(0)T}{\rho
s_{+}^{2}(0)s_{-}^{2}(0)}\qquad\left(\frac{s_{\pm}k}{T}\ll 1\right)\;.$ (73)

In the opposite limit of temperature much smaller than the sound-mode energy,
one has

$S_{11}(k)\simeq\left(\frac{\rho_{1}}{m_{1}}\right)\;\frac{[\varepsilon_{+}({\bf
k})\varepsilon_{-}({\bf k})+\varepsilon_{2}^{2}({\bf
k})]k^{2}}{2\rho\varepsilon_{+}({\bf k})\varepsilon_{-}({\bf
k})[\varepsilon_{+}({\bf k})+\varepsilon_{-}({\bf k})]}\;,$
$S_{12}(k)\simeq-\;\sqrt{\frac{\rho_{1}\rho_{2}}{m_{1}m_{2}}}\;\frac{\varepsilon_{12}^{2}({\bf
k})k^{2}}{2\rho\varepsilon_{+}({\bf k})\varepsilon_{-}({\bf
k})[\varepsilon_{+}({\bf k})-\varepsilon_{-}({\bf k})]}\;,$ (74)

provided that

$\frac{T}{s_{\pm}k}\ll 1\;.$

If both temperature and the sound-mode energy tend to zero so that
$T/s_{\pm}k$ tends to zero, then

$S_{11}({\bf
k})\simeq\left(\frac{\rho_{1}}{m_{1}}\right)\;\frac{[s_{+}(0)s_{-}(0)+s_{2}^{2}(0)]k}{2\rho
s_{+}(0)s_{-}(0)[s_{+}(0)+s_{-}(0)]}\;,$ $S_{12}({\bf
k})\simeq-\;\sqrt{\frac{\rho_{1}\rho_{2}}{m_{1}m_{2}}}\;\frac{s_{12}^{2}(0)k}{2\rho
s_{+}(0)s_{-}(0)[s_{+}(0)+s_{-}(0)]}\;,$ (75)

under the condition

$T\rightarrow 0\;,\qquad k\rightarrow 0\;,\qquad\frac{T}{s_{\pm}k}\rightarrow
0\;.$

In any case, when the mixture looses stability and starts separating, with
$s_{-}(0)$ tending to zero, the structure factors diverge.

##  9 Sum rules

One considers the integral

$K_{ij}({\bf k})=\int_{-\infty}^{\infty}\omega S_{ij}({\bf
k},\omega)\;\frac{d\omega}{2\pi}$ (76)

that, using the properties of the dynamic form-factor, can be written as

$K_{ij}({\bf k})=\frac{1}{\rho}\int_{0}^{\infty}\omega\gamma_{ij}({\bf
k},\omega)\;\frac{d\omega}{2\pi}\;.$ (77)

With the spectral functions (65), this yields

$K_{ij}({\bf
k})=\delta_{ij}\;\frac{\rho_{i}}{\rho}\;\left(\frac{k^{2}}{2m_{i}}\right)\;.$
(78)

The other important integral is

$Q_{ij}({\bf k})=\int_{-\infty}^{\infty}\frac{1}{\omega}\;S_{ij}({\bf
k},\omega)\;\frac{d\omega}{2\pi}\;,$ (79)

which can be rewritten as

$Q_{ij}({\bf
k})=\frac{1}{\rho}\int_{0}^{\infty}\frac{1}{\omega}\;\gamma_{ij}({\bf
k},\omega)\;\frac{d\omega}{2\pi}.$ (80)

Taking into account the spectral representation (56) gives

$Q_{ij}({\bf k})=-\;\frac{1}{\rho}\;{\rm Re}\;\chi_{ij}({\bf k},0)\;.$ (81)

Since the response function has the spectral representation

$\chi_{ij}({\bf k},\omega)=\int_{-\infty}^{\infty}\gamma_{ij}({\bf
k},\omega^{\prime})\;\left[\;\frac{1+n(\omega^{\prime})}{\omega-\omega^{\prime}+i0}\;-\;\frac{n(\omega^{\prime})}{\omega-\omega^{\prime}-i0}\;\right]\;\frac{d\omega^{\prime}}{2\pi}\;,$
(82)

where $n(\omega)\equiv(e^{\beta\omega}-1)^{-1}$, then using the property

$\frac{1}{\omega\pm i0}=P\;\frac{1}{\omega}\;\mp\;i\pi\delta(\omega)\;,$

where $P$ implies the principal value, we obtain the real part

${\rm Re}\;\chi_{ij}({\bf
k},\omega)=P\int_{-\infty}^{\infty}\frac{\gamma_{ij}({\bf
k},\omega^{\prime})}{\omega-\omega^{\prime}}\;\frac{d\omega^{\prime}}{2\pi}$
(83)

and the imaginary part

${\rm Im}\;\chi_{ij}({\bf k},\omega)=-\;\frac{1}{2}\;\gamma_{ij}{\bf
k},\omega)\;[\;1+2n(\omega)\;]\;.$ (84)

Here

$1+2n(\omega)=\frac{e^{\beta\omega}+1}{e^{\beta\omega}-1}=\coth\left(\frac{\beta\omega}{2}\right)\;.$

##  10 Isothermal compressibility

The response functions define the compressibility functions

$\varkappa_{i}({\bf r},t,{\bf
r}^{\prime},t^{\prime})\equiv-\;\frac{\chi_{ii}({\bf r},t,{\bf
r}^{\prime},t^{\prime})}{\rho_{i}({\bf r},t)\rho_{i}({\bf
r}^{\prime},t^{\prime})}\;,$ (85)

whose Fourier transforms are

$\varkappa_{i}({\bf k},\omega)=\int\chi_{i}({\bf r},t,0,0)\;e^{-i{\bf
k}\cdot{\bf r}+i\omega t}\;d{\bf r}dt\;.$ (86)

For an equilibrium uniform system, we have

$\varkappa_{i}({\bf k},\omega)=-\;\frac{\chi_{ii}({\bf
k},\omega)}{\rho_{i}^{2}}\;.$ (87)

The coefficient of isothermal compressibility is defined by the expression

$\varkappa_{i}\equiv{\rm Re}\;\varkappa_{i}(0,0)=-\;\frac{{\rm
Re}\chi_{ii}(0,0)}{\rho_{i}^{2}}\;.$ (88)

Hence it is given by the relation

$\varkappa_{i}=\frac{\rho}{\rho_{i}^{2}}\;Q_{ii}(0)\;.$ (89)

For a binary mixture, we get

${\rm
Re}\;\chi_{11}(0,0)=-\;\frac{\rho_{1}s_{2}^{2}(0)}{m_{1}s_{+}^{2}(0)s_{-}^{2}(0)}\;,$
${\rm
Re}\;\chi_{12}(0,0)=\sqrt{\frac{\rho_{1}\rho_{2}}{m_{1}m_{2}}}\;\frac{s_{12}^{2}(0)}{s_{+}^{2}(0)s_{-}^{2}(0)}\;.$
(90)

Therefore

$\varkappa_{1}=\frac{s_{2}^{2}(0)}{\rho_{1}m_{1}s_{+}^{2}(0)s_{-}^{2}(0)}\;,\qquad\varkappa_{2}=\frac{s_{1}^{2}(0)}{\rho_{2}m_{2}s_{+}^{2}(0)s_{-}^{2}(0)}\;.$
(91)

Thus we find another relation for the isothermal compressibilities

$\varkappa_{i}=\frac{\rho S_{ii}(0)}{\rho_{i}^{2}T}\;.$ (92)

When the mixture looses its stability and $s_{-}(0)$ tends to zero, the
compressibility coefficients diverge.

##  11 Thermodynamic stability

Thermodynamic stability conditions are based on the comparison of the free
energies of the mixed, $F_{mix}$, and separated, $F_{sep}$, systems. The
mixture is stable when

$F_{mix}<F_{sep}\;.$ (93)

In what follows, we consider a uniform system with the component densities
$\rho_{i}({\bf r})=N_{i}/V$ and use the Hartree approximation.

The free energy of a mixed system writes as

$F_{mix}=E_{mix}-TS_{mix}\;.$ (94)

The energy

$E_{mix}=\langle\;H\;\rangle+\sum_{i}\mu_{i}N_{i}=K_{mix}+E_{mix}^{ext}+E_{mix}^{int}$
(95)

includes the kinetic energy term $K_{mix}$, the energy due to external
potentials,

$E_{mix}^{ext}=\sum_{i}\int_{V}U_{i}({\bf r})\;\rho_{i}({\bf r})\;d{\bf
r}=\sum_{i}U_{i}\;N_{i}\;,$ (96)

where

$U_{i}\equiv\frac{1}{V}\int_{V}U_{i}({\bf r})\;d{\bf r}\;,$

and the interaction energy

$E_{mix}^{int}=\frac{1}{2}\sum_{ij}\int_{V}\rho_{i}({\bf r})\;\Phi_{ij}({\bf
r}-{\bf r}^{\prime})\;\rho_{j}({\bf r}^{\prime})\;d{\bf r}d{\bf
r}^{\prime}\;.$ (97)

For a uniform system, the latter is

$E_{mix}^{int}=\frac{1}{2}\sum_{ij}\Phi_{ij}\;\frac{N_{i}N_{j}}{V}\;,$ (98)

where

$\Phi_{ij}\equiv\int_{V}\Phi_{ij}({\bf r})\;d{\bf r}\;.$ (99)

The Hamiltonian of a separated system is

$H_{sep}=\sum_{i}\int\psi_{i}^{\dagger}({\bf
r})\;\left(-\;\frac{\nabla^{2}}{2m}+U_{i}-\mu_{i}\right)\psi_{i}({\bf
r})\;d{\bf r}\;+$ $+\;\frac{1}{2}\sum_{i}\int\psi_{i}^{\dagger}({\bf
r})\;\psi_{i}^{\dagger}({\bf r}^{\prime})\;\Phi_{ii}({\bf r}-{\bf
r}^{\prime})\;\psi_{i}({\bf r}^{\prime})\;\psi_{i}({\bf r})\;d{\bf r}d{\bf
r}^{\prime}\;.$ (100)

The free energy writes as

$F_{sep}=E_{sep}-TS_{sep}\;.$ (101)

The energy

$E_{sep}=\langle\;H_{sep}\;\rangle+\sum_{i}\mu_{i}N_{i}=K_{sep}+E_{sep}^{ext}+E_{sep}^{int}$
(102)

contains the kinetic energy $K_{sep}$, the energy caused by external
potentials $E_{sep}^{ext}$, and the interaction energy $E_{sep}^{int}$.

The assumption of the system uniformity presupposes that the system volumes
are sufficiently large, so that the equalities are valid:

$\frac{1}{V}\int_{V}U_{i}({\bf r})\;d{\bf
r}=\frac{1}{V_{i}}\int_{V_{i}}U_{i}({\bf r})\;d{\bf
r}\;,\qquad\int_{V}\Phi_{ii}({\bf r})\;d{\bf r}=\int_{V_{i}}\Phi_{ii}({\bf
r})\;d{\bf r}\;.$

Also, we assume that the difference in the kinetic energy, under the same
given conditions, can be neglected as compared to the difference in the
potential energy, thus setting $K_{mix}-K_{sep}\approx 0$. The component
densities of the uniform separated system are

$\rho_{i}^{loc}({\bf r})=\langle\;\psi_{i}^{\dagger}({\bf r})\;\psi_{i}({\bf
r})\;\rangle=\frac{N_{i}}{V_{i}}\;.$ (103)

Then the energy due to external potentials is

$E_{sep}^{ext}=\sum_{i}\int_{V_{i}}U_{i}({\bf r})\;\rho_{i}^{loc}({\bf
r})\;d{\bf r}=\sum_{i}U_{i}\;N_{i}\;,$ (104)

while the interaction energy of the separated system becomes

$E_{sep}^{int}=\frac{1}{2}\sum_{i}\int_{V_{i}}\rho_{i}^{loc}({\bf
r})\;\Phi_{ii}({\bf r}-{\bf r}^{\prime})\;\rho_{i}^{loc}({\bf
r}^{\prime})\;d{\bf r}d{\bf r}^{\prime}\;,$ (105)

which gives

$E_{sep}^{int}=\frac{1}{2}\sum_{i}\Phi_{ii}\;\frac{N_{i}^{2}}{V_{i}}\;.$ (106)

The condition of thermodynamic stability (93) takes the form

$E_{mix}^{int}-E_{sep}^{int}<N\;T\;\Delta S_{mix}\;,$ (107)

in which

$\Delta S_{mix}=\frac{1}{N}\;(S_{mix}-S_{sep})$ (108)

is the entropy of mixing per particle.

Mechanical equilibrium implies the equality of the pressures

$P_{i}=-\;\frac{\partial F_{sep}}{\partial V_{i}}=P_{j}=-\;\frac{\partial
F_{sep}}{\partial V_{j}}.$ (109)

From the latter condition, we have

$\frac{V_{j}}{V_{i}}=\frac{N_{j}}{N_{i}}\;\sqrt{\frac{\Phi_{jj}}{\Phi_{ii}}}\;.$

This allows us to rewrite the stability condition (107) as the inequality

$\sum_{ij}\frac{N_{i}N_{j}}{2V}\;\left(\Phi_{ij}\;-\;\sqrt{\Phi_{ii}\Phi_{jj}}\;-\;\frac{2T}{\rho}\;\Delta
S_{mix}\right)<0\;.$ (110)

Assuming that this inequality is valid for any $N_{i}$ and $N_{j}$, we come to
the stability condition [30]

$\Phi_{ij}\;-\;\sqrt{\Phi_{ii}\;\Phi_{jj}}\;<\;\frac{2T}{\rho}\;\Delta
S_{mix}\;.$ (111)

The entropy of mixing per particle can be represented as

$\Delta S_{mix}=-\sum_{i}n_{i}\;\ln
n_{i}\qquad\left(n_{i}\equiv\frac{N_{i}}{N}\right)\;.$ (112)

The mixture stability condition (111) shows that finite temperatures
facilitate the component mixing. The components could be immiscible at zero
temperature but become miscible at finite temperature. This is also in
agreement with the miscibility criterion for Bose-condensed mixtures, taking
into account strong interactions and temperature [26]. The miscibility
condition of type (34) can be invalid, but the components, anyway, can be
miscible at finite temperatures and interactions.

##  12 Conclusion

Mixtures of quantum fluids, liquids or gases, are considered. The studied
systems are treated as normal quantum systems, that is requiring quantum
techniques for their description, but exhibiting no superfluidity or
superconductivity. The main point of interest in the present review is the
conditions of quantum mixture stability with respect to the spatial separation
into different components. The mixture stability can be studied from the point
of view of dynamic or thermodynamic stability.

The study of dynamic stability is based on the investigation of spectra of
collective excitations. The principal difference between Bose-condensed and
normal fluids is in the following. In a Bose-condensed system, the global
gauge symmetry is broken. This symmetry breaking is conveniently characterized
by the Bogolubov shift [40, 41, 42]. Then the vacuum state is described by a
coherent state formed by completely condensed particles, while the uncondensed
particles play the role of excitations above the vacuum [43]. Because of this,
the single-particle Green functions of uncondensed particles, actually,
describe excitations above the coherent-state vacuum. This is why to study the
collective spectrum of a Bose-condensed system it is sufficient to consider
solely the single-particle Green functions of uncondensed particles, whose
poles give the required collective spectrum.

On the other hand, the vacuum for a normal quantum system is an empty state.
The poles of single-particle Green functions describe the single-particle
spectra, but the spectrum of collective excitations is characterized by the
poles of two-particle Green functions or by the poles of response functions.

For the description of response functions, one usually involves the random-
phase approximation. This approximation describes well the low-temperature
states, but the influence of temperature on the spectrum of collective
excitations is not as important.

The role of temperature can be investigated by resorting to the study of
thermodynamic stability, when one compares the free energies of a mixed and
separated systems. The found stability criterion demonstrates that temperature
facilitates the mixing of different quantum components. At low temperature,
the system could be immiscible, while becomes miscible when rising
temperature.

This effect also is in agreement with the study of stability of Bose-condensed
systems with strong interactions and finite temperature. For the correct
description of such systems, it is necessary to take into account anomalous
averages and the existence of two chemical potentials for each of the
components, which is caused by the global gauge symmetry breaking [44, 45,
46]. Including into the description the corrections above the coherent
approximation implies the inclusion into consideration of quantum and thermal
fluctuations. These fluctuations make the mixing easier. The Bose-condensed
mixture stability, accurately taking account of fluctuations, has been
recently analyzed [26], confirming that fluctuations facilitate the mixing of
the components.

## References

  * [1] Letokhov V 2007 Laser Control of Atoms and Molecules (New York: Oxford University)
  * [2] Lieb E H, Seiringer R, Solovej J P and Yngvason J 2005 The Mathematics of the Bose Gas and its Condensation (Basel: Birkhauser)
  * [3] Pethick C J and Smith H 2008 Bose-Einstein Condensation in Dilute Gases (Cambridge: Cambridge University)
  * [4] Ueda M 2010 Fundamentals and New Frontiers of Bose-Einstein Condensation (Singapore: World Scientific)
  * [5] Courteille P W, Bagnato V S and Yukalov V I 2001 Laser Phys. 11 659
  * [6] Andersen J O 2004 Rev. Mod. Phys. 76 599
  * [7] Yukalov V I 2004 Laser Phys. Lett. 1 435
  * [8] Bongs K and Sengstock K 2004 Rep. Prog. Phys. 67 907
  * [9] Yukalov V I and Girardeau M D 2005 Laser Phys. Lett. 2 375
  * [10] Posazhennikova A 2006 Rev. Mod. Phys. 78 1111
  * [11] Morsch O and Oberthaler M 2006 Rev. Mod. Phys. 78 179
  * [12] Yukalov V I 2007 Laser Phys. Lett. 4 632
  * [13] Moseley C, Fialko O and Ziegler K 2008 Ann. Phys. (Berlin) 17 561
  * [14] Bloch I, Dalibard J and Zwerger W 2008 Rev. Mod. Phys. 80 885
  * [15] Proukakis N P and Jackson B 2008 J. Phys. B 41 203002
  * [16] Yurovsky V A, Olshanii M and Weiss D S 2008 Adv. At. Mol. Opt. Phys. 55 61
  * [17] Castilho P C M, E Pedrozo-Penafiel E, Gutierrez E M, Mazo P L, Roati G, Farias K M and Bagnato V S 2019 Laser Phys. Lett. 16 035501
  * [18] Saarela M and Taipaleenmaki T 2003 Int. J. Mod. Phys. B 17 5227
  * [19] Adhikari S K 2007 Phys. Rev. A 76 053609
  * [20] Yukalov V I and Yukalova E P 2004 Laser Phys. Lett. 1 50
  * [21] Abad M, Recati1 A, Stringari1 S and Chevy F 2015 Eur. Phys. J. D 69 126
  * [22] Nepomnyashchii Y A 1974 Theor. Math. Phys. 20 399
  * [23] Hugenholtz N M and Pines D 1959 Phys. Rev. 116 489
  * [24] Nepomnyashchii Y A 1976 J. Exp. Theor. Phys. 43 559
  * [25] Watabe S 2021 Phys. Rev. A 103 053307
  * [26] Rakhimov A, Abdurakhmonov T, Narzikulov Z and Yukalov V I 2022 Phys. Rev. A 106 033301
  * [27] Gavoret T and Nozières P 1964 Ann. Phys. (N.Y.) 28 349
  * [28] Watabe S 2020 New J. Phys. 22 103010
  * [29] Kürten K E and Campbell C E 1982 Phys. Rev. B 26 124
  * [30] Yukalov V I 2016 Laser Phys. 26 062001
  * [31] Kirkwood J G 1965 Quantum Statistics and Cooperative Phenomena (New York: Gordon and Breach)
  * [32] Yukalov V I 1990 Phys. Rev. A 42 3324
  * [33] Yukalov V I 2016 Phys. Rev. E 94 102106
  * [34] Kadanoff L P and Baym G 1962 Quantum Statistical Mechanics (New York: Benjamin)
  * [35] Yukalov V I 1980 Acta Phys. Pol. A 57 295
  * [36] Yukalov V I 2018 Laser Phys. 28 053001
  * [37] Henrici P 1988 Applied and Computational Complex Analysis Vol. 1 (New York: Wiley)
  * [38] Yukalov V I 2007 J. Phys. Stud. 11 55
  * [39] Yukalov V I 1998 Statistical Green’s Functions (Kingston: Queen’s University)
  * [40] Bogolubov N N 1967 Lectures on Quantum Statistics (New York: Gordon and Breach) Vol. 1
  * [41] Bogolubov N N 1970 Lectures on Quantum Statistics (New York: Gordon and Breach) Vol. 2
  * [42] Bogolubov N N 2015 Quantum Statistical Mechanics (Singapore: World Scientific)
  * [43] Yukalov V I 2006 Laser Phys. 16 511
  * [44] Yukalov V I 2005 Phys. Rev. E 72 066119
  * [45] Yukalov V I 2006 Phys. Lett. A 359 712
  * [46] Yukalov V I 2008 Ann. Phys. (N.Y.) 323 461