**Preface to "Symmetry and Mesoscopic Physics"**

Nature loves to play dice with us, providing this or that riddle along the way of our understanding of its diversity. It is like the faces of an n-dimensional facially regular polyhedron, the numerous faces of which from time to time indicate different problems for us. Depending on our luck, we can get the Ariadne's thread, which could lead us to the laws of nature based on beauty and harmony. Since ancient times people have realized that very often harmony is associated with symmetry. In this issue, we tried to present some of the results of the manifestation of symmetries and symmetry breaking in finite quantum systems that have peculiarities compared to macroscopic samples. Some authors find the effects of symmetry and asymptotic symmetry breaking in the analysis of Bose-Einstein condensate. Other authors have shown how symmetries, due to spin-orbit interaction in quantum dots and graphene, can explain the effects of electron transport in these systems. The problem of scaling symmetry, related to fractals, is studied as well, which allows us to understand the mechanisms of fabrication of complex materials with predefined physical properties and functionalities.

> **Rashid G. Nazmitdinov, Vyacheslav Yukalov** *Editors*

### *Article* **Particle Fluctuations in Mesoscopic Bose Systems**

#### **Vyacheslav I. Yukalov 1,2**


Received: 10 April 2019; Accepted: 25 April 2019; Published: 1 May 2019

**Abstract:** Particle fluctuations in mesoscopic Bose systems of arbitrary spatial dimensionality are considered. Both ideal Bose gases and interacting Bose systems are studied in the regions above the Bose–Einstein condensation temperature *Tc*, as well as below this temperature. The strength of particle fluctuations defines whether the system is stable or not. Stability conditions depend on the spatial dimensionality *d* and on the confining dimension *D* of the system. The consideration shows that mesoscopic systems, experiencing Bose–Einstein condensation, are stable when: (i) ideal Bose gas is confined in a rectangular box of spatial dimension *d* > 2 above *T<sup>c</sup>* and in a box of *d* > 4 below *Tc*; (ii) ideal Bose gas is confined in a power-law trap of a confining dimension *D* > 2 above *T<sup>c</sup>* and of a confining dimension *D* > 4 below *Tc*; (iii) the interacting Bose system is confined in a rectangular box of dimension *d* > 2 above *Tc*, while below *Tc*, particle interactions stabilize the Bose-condensed system, making it stable for *d* = 3; (iv) nonlocal interactions diminish the condensation temperature, as compared with the fluctuations in a system with contact interactions.

**Keywords:** Bose systems; asymptotic symmetry breaking; Bose–Einstein condensation; particle fluctuations; stability of Bose systems

#### **1. Introduction**

The theory of Bose systems has recently attracted high attention triggered by experimental studies of cold trapped atoms (see, e.g., the books and review articles [1–19]). Special attention has been payed to the study of particle fluctuations, mainly considering three-dimensional macroscopic Bose systems or harmonically-trapped atoms. The importance of this problem has been emphasized after the appearance of a number of papers claiming the occurrence of thermodynamically-anomalous particle fluctuations in the whole region below the condensation temperature *T<sup>c</sup>* even for equilibrium three-dimensional interacting systems (a list of the papers containing such claims has been summarized in [20]). The origin of the arising fictitious anomalies and the ways of avoiding them have been discussed in detail in reviews [16–18].

It would not be strange if anomalously strong fluctuations would be found at the point of a second-order phase transition. This would be natural, since at the point of a phase transition, the system is unstable and fluctuations in a system can drastically increase. It is exactly the system instability that drives the phase transition and forces the system to transfer to another state. However, as soon as the transition to the other state has happened, the real system becomes stable and has to exhibit thermodynamically normal fluctuations. It is therefore more than strange how thermodynamically-anomalous fluctuations could arise in realistic three-dimensional interacting systems.

Moreover, Bose–Einstein condensation is necessarily accompanied by the spontaneous breaking of global gauge symmetry. From the mathematical point of view, the similar breaking of continuous symmetry occurs under magnetic phase transitions [21]; hence, anomalous fluctuations of the order parameter should appear in magnets below *Tc*. However, thermodynamically-anomalous fluctuations imply the system instability [22]. Therefore, if such fluctuations would really arise in the whole range below *Tc*, then neither superfluids nor magnets would exist. Fortunately, it has been shown [23,24] that thermodynamically-anomalous fluctuations in interacting three-dimensional equilibrium systems, discussed in theoretical papers, are just calculational artifacts caused, briefly speaking, by the use of a second-order approximation for calculating fourth-order terms.

The aim of the present paper is to extend the investigation of particle fluctuations in Bose systems in several aspects: First, we consider mesoscopic systems that are finite, although containing many particles *N* 1. Taking into account a finite number of particles requires modifying the definition of the Bose function by introducing a finite cutoff responsible for the existence of a minimal wave vector prescribed by the system geometry. Second, we analyze particle fluctuations above, as well as below *T<sup>c</sup>* for the Bose systems of arbitrary dimensionality, which allows us to find the critical spatial dimension above which the system is stable. Third, we consider two types of Bose systems, confined either in a rectangular box or in a power-law trap. Fourth, the influence of nonlocal interactions on particle fluctuations is analyzed, as compared to that of local interactions.

Throughout the paper, the system of units is employed where the Planck and Boltzmann constants are set to one, ¯*h* = 1 and *k<sup>B</sup>* = 1.

#### **2. Particle Fluctuations and Stability**

Here and in what follows, we consider mesoscopic systems that are finite, containing a finite number of particles *N*, although with this number is rather large, *N* 1.

Observable quantities are given by statistical averages <sup>h</sup>*A*ˆ<sup>i</sup> of Hermitian operators *<sup>A</sup>*ˆ. Fluctuations of the observable quantities are characterized by the variance:

$$\text{var}(\hat{A}) \equiv \langle \hat{A}^2 \rangle - \langle \hat{A} \rangle^2.$$

The observable is called extensive when:

$$
\langle \hat{A} \rangle \propto N \qquad (N \gg 1) \,\tag{1}
$$

which is equivalent to the condition:

$$\frac{\langle \hat{A} \rangle}{N} \simeq \text{const} \qquad \left( N \gg 1 \right). \tag{2}$$

Fluctuations are termed thermodynamically normal if the inequalities:

$$0 \le \frac{\text{var}(\hat{A})}{|\langle \hat{A} \rangle|} < \infty \tag{3}$$

are valid for any *N*, which can also be represented as the condition:

$$\frac{|\text{var}(\hat{A})|}{|\langle \hat{A} \rangle|} \simeq \text{const} \qquad (\text{N} \gg 1) \,. \tag{4}$$

When these conditions do not hold, the fluctuations are called thermodynamically anomalous. Sometimes, instead of the terms thermodynamically normal or thermodynamically anomalous, one says, for short, that fluctuations are just normal or anomalous.

Particle fluctuations, describing the fluctuations of the number of particles, characterized by the number-of-particles operator *N*ˆ , are quantified by the relative variance:

$$\frac{\text{var}(\hat{\mathsf{N}})}{N} = \frac{1}{N} \left( \langle \hat{\mathsf{N}}^2 \rangle - \langle \hat{\mathsf{N}} \rangle^2 \right) \tag{5}$$

where *<sup>N</sup>* <sup>=</sup> <sup>h</sup>*N*<sup>ˆ</sup> <sup>i</sup>. The fluctuations are normal when:

$$0 \le \frac{\text{var}(\hat{N})}{N} < \infty \tag{6}$$

for any *N*, or in other words, when:

$$\frac{\text{var}(\hat{N})}{N} \simeq \text{const} \qquad (N \gg 1) \,. \tag{7}$$

The strength of particle fluctuations characterizes the system stability, since these fluctuations are directly connected to the isothermal compressibility:

$$\kappa\_T \equiv -\frac{1}{V} \left( \frac{\partial V}{\partial P} \right)\_{TN} = \frac{1}{\rho N} \left( \frac{\partial N}{\partial \mu} \right)\_{TV} \tag{8}$$

by the equality:

$$\kappa\_T = \frac{\text{var}(\hat{N})}{\rho T N} \qquad \left(\rho \equiv \frac{N}{V}\right) \text{ .} \tag{9}$$

with *ρ* being the average particle density. The system stability requires that:

$$0 \le \varkappa\_T < \infty \tag{10}$$

for any *N*, which yields Conditions (6) and (7). The above relations give us one of the ways for calculating the relative variance:

$$\frac{\text{var}(\hat{N})}{N} = \rho T \kappa\_T = \frac{T}{N} \left( \frac{\partial N}{\partial \mu} \right)\_{TV} \,. \tag{11}$$

#### **3. Ideal Gas in a Rectangular Box**

Bose systems in a rectangular box are not merely an interesting object allowing for detailed calculations, but it can also be realized experimentally inside box-shaped traps [25–27].

#### *3.1. Modified Bose Function*

The grand Hamiltonian of a gas in a rectangular box of volume *V* reads as:

$$H = \hat{H} - \mu \hat{N} = \int \psi^\dagger(\mathbf{r}) \left( -\frac{\nabla^2}{2m} - \mu \right) \psi(\mathbf{r}) \, d\mathbf{r} \, \tag{12}$$

where the integration is over the given volume *V*. Assuming periodic continuation of the box, the field operators can be expanded in plane waves,

$$
\psi(\mathbf{r}) = \sum\_{k} a\_{k} \varphi\_{k}(\mathbf{r}) \; , \qquad \varphi\_{k}(\mathbf{r}) = \frac{1}{\sqrt{V}} \; \varepsilon^{i\mathbf{k} \cdot \mathbf{r}} \; , \tag{13}
$$

which gives:

$$H = \sum\_{k} \omega\_{k} a\_{k}^{\dagger} a\_{k} \qquad \left(\omega\_{k} = \frac{k^{2}}{2m} - \mu\right) \,. \tag{14}$$

The total number of particles is the sum:

$$N = N\_0 + N\_1 \ \ \ \ \quad \quad N\_1 = \sum\_{k \neq k\_0} n\_k \ \ \ \ \ \tag{15}$$

where *N*<sup>0</sup> is the number of condensed particles, while *N*<sup>1</sup> is the number of uncondensed particles, with the momentum distribution:

$$n\_k \equiv \langle a\_k^\dagger a\_k \rangle = \left( e^{\oint \omega\_k} - 1 \right)^{-1}. \tag{16}$$

Here, *β* = 1/*T* is the inverse temperature. For a large number of particles, the sums over momenta can be represented as the integrals,

$$
\sum\_{k} n\_k \to V \int n\_k \, \frac{d\mathbf{k}}{(2\pi)^d} \,\,\prime
$$

where *d* is spatial dimensionality. In the case of isotropic functions under the integrals, it is possible to pass to spherical coordinates. However, it is necessary to take into account that for a finite system, the values of the wave vectors start not from zero, but from a finite minimal momentum *k*<sup>0</sup> that can be estimated as:

$$k\_0 = \frac{2\pi}{L} = \frac{2\pi}{aN^{1/d}}\,\,\,\,\,\tag{17}$$

with the box volume:

$$V = L^d \,, \qquad L = aN^{1/d} \,,$$

where *a* is the mean interparticle distance. Thus, the integration over the momenta takes the form:

$$\int \frac{d\mathbf{k}}{(2\pi)^d} \rightarrow \frac{2}{(4\pi)^{d/2}\Gamma(d/2)} \int\_{k\_0}^{\infty} k^{d-1} \, dk \, \tag{18}$$

where the lower limit is given by the cutoff prescribed by the minimal quantity *k*0. Then, the number of uncondensed particles becomes proportional to the modified Bose function:

$$g\_n(z) \equiv \frac{1}{\Gamma(n)} \int\_{u\_0}^{\infty} \frac{zu^{n-1}}{e^u - z} \, du\,,\tag{19}$$

with *z* ≡ exp(*βµ*) being the fugacity and where the lower limit is given by the cutoff:

$$
\mu\_0 = \frac{k\_0^2}{2mT} = \frac{\varepsilon\_0}{T} \tag{20}
$$

defined by the minimal energy:

$$
\varepsilon\_0 = \frac{2\pi^2}{ma^2} \,\mathrm{N}^{-2/d} \,. \tag{21}
$$

Since the minimal energy (21) tends to zero for large *N*, it is admissible to keep in mind that:

$$
u\_0 \ll 1 \qquad \left(N \gg 1\right). \tag{22}$$

In this way, the relative variance (11) can be expressed through the derivative of the modified Bose function (19). The latter differs from the standard Bose function by the existence of a nonzero lower integration limit defined by the minimal wave vector.

#### *3.2. Fluctuations above the Condensation Temperature*

At temperatures above the condensation point, there are no condensed particles, so that the total number of particles reads as:

$$N = \frac{V}{\lambda\_T^d} \text{ g}\_{d/2}(\mathbf{z}) \qquad (T \ge T\_\mathbf{c}) \text{ .} \tag{23}$$

where:

$$
\lambda\_T \equiv \sqrt{\frac{2\pi}{mT}}
$$

is the thermal wavelength. Hence, the relative variance (11) is:

$$\frac{\text{var}(\mathcal{N})}{N} = \frac{z}{\rho \lambda\_T^d} \frac{\partial g\_{d/2}(z)}{\partial z} \qquad (T > T\_c) \; \tag{24}$$

where *ρ* ≡ *N*/*V* is the particle density.

Estimating the Bose function above *Tc*, where *z* < 1, we find:

$$\log\_{\mathbb{H}}(z) = -\frac{z}{(1-z)\Gamma(1+n)} \left[ u\_0^n - \frac{nu\_0^{1+n}}{(1+n)(1-z)} \right] \qquad (n < 0 \text{ } z < 1) \,. \tag{25}$$

In particular,

$$\log\_{-1/2}(z) = -\frac{z}{\sqrt{\pi}(1-z)} \left( u\_0^{-1/2} + \frac{u\_0^{1/2}}{1-z} \right) \qquad (z < 1) \tag{26}$$

and:

$$\mathbf{g}\_0(\mathbf{z}) = -\frac{z}{1-z} \qquad (z < 1) \,. \tag{27}$$

Calculating the derivatives of the modified Bose functions requires being attentive, since some of the derivatives are different from those for the standard Bose functions. Generally, we have:

$$\frac{\partial g\_n(z)}{\partial z} = \frac{1}{z} \, g\_{n-1}(z) + \frac{u\_0^{n-1}}{\Gamma(n)(1 - z + u\_0)} \,. \tag{28}$$

Using the smallness of *u*0, we can write:

$$\frac{u\_0^{n-1}}{1 - z + \mu\_0} \asymp \frac{1}{1 - z} \left( u\_0^{n-1} - \frac{u\_0^n}{1 - z} \right) \qquad (z < 1) \; .$$

Therefore, for *n* < 1, we find:

$$\frac{\partial g\_n(z)}{\partial z} = -\frac{u\_0^n}{(1-z)^2 \Gamma(1+n)} \qquad (n < 1 \text{ , } z < 1) \text{ ;}\tag{29}$$

while for *n* > 1, keeping the main terms, we get:

$$\frac{\partial g\_n(z)}{\partial z} = \frac{1}{z} \, g\_{n-1}(z) \qquad (n > 1 \,\, \, z < 1) \,. \tag{30}$$

We shall also need the derivatives:

$$\frac{\partial g\_{1/2}(z)}{\partial z} = -\frac{2u\_0^{1/2}}{\sqrt{\pi}(1-z)^2} \qquad (d=1, \; z < 1),$$

$$\frac{\partial g\_1(z)}{\partial z} = -\frac{u\_0}{(1-z)^2} \qquad (d=2, \; z < 1),$$

$$\frac{\partial g\_{3/2}(z)}{\partial z} = \frac{1}{z} \; g\_{1/2}(z) \qquad (d=3, \; z < 1).$$

For the relative variance, depending on the space dimensionality, we obtain:

$$\frac{\text{var}(\hat{N})}{N} = -\frac{2z}{\sqrt{\pi}(1-z)^2 \rho \lambda\_T} \left(\frac{\varepsilon\_0}{T}\right)^{1/2} \qquad (d=1, \ T > T\_c),$$

$$\frac{\text{var}(\hat{N})}{N} = -\frac{z}{(1-z)^2 \rho \lambda\_T^2} \left(\frac{\varepsilon\_0}{T}\right) \qquad (d=2, \ T > T\_c),$$

*Symmetry* **2019**, *11*, 603

$$\frac{\text{var}(\hat{N})}{N} = \frac{1}{\rho \lambda\_T^3} \ g\_{1/2}(z) \qquad (d = 3 \,\,\, T > T\_\mathfrak{c}) \,. \tag{31}$$

The negative values of the variance for *d* = 1 and *d* = 2 show that these low-dimensional systems are unstable. However, the gas is stable in three dimensions.

As follows from the derivative:

$$\frac{\partial g\_{d/2}(z)}{\partial z} = \frac{1}{z} \, g\_{(d-2)/2}(z) \qquad (d > 2 \,\, , \, z < 1) \,. \tag{32}$$

the system is stable for *d* > 2. That is, the critical spatial dimension, above which the uncondensed gas in a rectangular box is stable, is *d<sup>c</sup>* = 2, so that the stability condition is:

$$d\_{\mathbb{C}} > d\_{\mathbb{C}} = 2 \qquad \left(T > T\_{\mathbb{C}}\right). \tag{33}$$

#### *3.3. Condensation Temperature of a Gas in a Rectangular Box*

At the temperature of Bose condensation, the chemical potential becomes zero, *µ* = 0, because of whuich *z* = 1. The total number of particles:

$$N = \frac{V}{\lambda\_T^d} \text{ g}\_{d/2}(1) \qquad (T = T\_c) \tag{34}$$

defines the critical temperature:

$$T\_c = \frac{2\pi}{m} \left[\frac{\rho}{g\_{d/2}(1)}\right]^{2/d}.\tag{35}$$

For different dimensionalities, we have:

$$g\_{1/2}(1) = \frac{2}{\sqrt{\pi}} \ u\_0^{-1/2} \qquad (d=1)$$

$$g\_1(1) = -\ln u\_0 \qquad (d=2)$$

$$g\_{3/2}(1) = \zeta(3/2) \qquad (d=3)$$

This gives the critical temperatures for a one-dimensional system:

$$T\_{\mathfrak{c}} = \frac{\pi \mathfrak{p}}{\sqrt{2m}} \,\, \mathfrak{e}\_0^{1/2} \qquad (d=1) \; , \tag{36}$$

and for a two-dimensional system:

$$T\_c = \frac{2\pi\rho}{m\ln(T\_c/\varepsilon\_0)}\qquad(d=2)\,. \tag{37}$$

Iterating the latter equation and taking into account that:

$$\frac{T\_c}{\varepsilon\_0} \ll \exp\left(\frac{2\pi\rho}{m\varepsilon\_0}\right), \qquad \varepsilon\_0 \propto N^{-2/d} \qquad (N \gg 1)$$

we obtain:

$$T\_{\varepsilon} = \frac{2\pi\rho}{m\ln(2\pi\rho/m\varepsilon\_0)}\qquad(d=2)\,. \tag{38}$$

For a three-dimensional box, we have the known result:

$$T\_c = \frac{2\pi}{m} \left[\frac{\rho}{\zeta(3/2)}\right]^{2/3} \qquad (d=3) \,\text{.}\tag{39}$$

The critical temperatures at a large *N* 1 scale as:

$$T\_c \propto \frac{1}{N} \qquad (d=1)$$

$$T\_c \propto \frac{1}{\ln N} \qquad (d=2)$$

$$T\_c \propto \text{const} \qquad (d=3)$$

For *d* ≤ 2, the critical temperature diminishes to zero as *N* increases. It remains finite for *d* > 2. Recall that, as is found above, the system is unstable for *d* ≤ 2. Thus, the Bose gas in a box is stable in the case where the critical temperature remains finite for large *N*.

#### *3.4. Fluctuations below Critical Temperature*

Below the critical temperature, there appears the Bose–Einstein condensate, so that the number-of-particle operator becomes the sum of the number of condensed particles *N*<sup>0</sup> and the number-of-particle operator *N*ˆ <sup>1</sup> of uncondensed particles,

$$
\hat{N} = N\_0 + \hat{N}\_1 \qquad \left( T < T\_c \right) \,. \tag{41}
$$

When the condensate function *η* is introduced by means of the Bogolubov shift [28–30] of the field operator:

$$
\psi(\mathbf{r}) \to \eta(\mathbf{r}) + \psi\_1(\mathbf{r}) \,\mathrm{s}
$$

then particle fluctuations are defined by the fluctuations of uncondensed particles (see the detailed explanations in [3,9,16–18]),

$$\text{var}(\hat{N}) = \text{var}(\hat{N}\_1) \dots$$

The average number of uncondensed particles is:

$$N\_1 = \frac{V}{\lambda\_T^d} \ g\_{d/2}(1) \qquad \left(T < T\_c\right). \tag{42}$$

Therefore, the relative particle variance reads as:

$$\frac{\text{var}(\hat{N})}{N} = \frac{1}{\rho \lambda\_T^d} \lim\_{z \to 1} \frac{\partial g\_{d/2}(z)}{\partial z} \qquad (T < T\_c) \,. \tag{43}$$

For a mesoscopic system, we have:

$$\lim\_{z \to 1} \frac{\partial g\_n(z)}{\partial z} = g\_{n-1}(1) + \frac{u\_0^{n-2}}{\Gamma(n)}.\tag{44}$$

Notice that the last term here would be absent for a macroscopic system. In particular,

$$\lim\_{z \to 1} \frac{\partial g\_{1/2}(z)}{\partial z} = g\_{-1/2}(1) + \frac{u\_0^{-3/2}}{\sqrt{\pi}} \, .$$

Since:

$$\mathcal{g}\_{-1/2}(1) = -\frac{1}{3\sqrt{\pi}} \,\mu\_0^{-3/2} \,\mu$$

we find:

$$\lim\_{z \to 1} \frac{\partial g\_{1/2}(z)}{\partial z} = \frac{2}{3\sqrt{\pi}} \,\,\mu\_0^{-3/2} \,\, .$$

*Symmetry* **2019**, *11*, 603

We also need the limits:

$$\lim\_{z \to 1} \frac{\partial \mathcal{g}\_1(z)}{\partial z} = \frac{1}{u\_0} \text{ }.$$

and:

$$\lim\_{z \to 1} \frac{\partial g\_{3/2}(z)}{\partial z} = g\_{1/2}(1) + \frac{2}{\sqrt{\pi}} \ u\_0^{-1/2} \ .$$

Using *g*1/2(1) from the previous subsection, we get:

$$\lim\_{z \to 1} \frac{\partial g\_{3/2}(z)}{\partial z} = \frac{4}{\sqrt{\pi}} \,\mu\_0^{-1/2} \,.$$

Generally, from the expression:

$$\lim\_{z \to 1} \frac{\partial g\_{d/2}(z)}{\partial z} = g\_{(d-2)/2}(1) + \frac{u\_0^{(d-4)/2}}{\Gamma(d/2)}\tag{45}$$

we see that the last term here increases with *N* by a power-law, when *d* < 4, while it increases logarithmically for *d* = 4,

$$\lim\_{z \to 1} \frac{\partial g\_2(z)}{\partial z} = 1 - \ln u\_0 \qquad (d = 4) \; .$$

For the relative variance (43), we obtain:

$$
\frac{\text{var}(\hat{N})}{N} = \frac{2}{3\sqrt{\pi}\rho\lambda\_T} \left(\frac{T}{\varepsilon\_0}\right)^{3/2} \qquad (d=1),
$$

$$
\frac{\text{var}(\hat{N})}{N} = \frac{1}{\rho\lambda\_T^2} \left(\frac{T}{\varepsilon\_0}\right) \qquad (d=2),
$$

$$
\frac{\text{var}(\hat{N})}{N} = \frac{4}{\sqrt{\pi}\rho\lambda\_T^3} \left(\frac{T}{\varepsilon\_0}\right)^{1/2} \qquad (d=3),
$$

$$
\frac{\text{var}(\hat{N})}{N} = \frac{1}{\rho\lambda\_T^4} \left(\frac{T}{\varepsilon\_0}\right) \qquad (d=4).
\tag{46}
$$

Keeping in mind that *ε*<sup>0</sup> ∝ *N*−2/*<sup>d</sup>* , the scaling of these expressions with respect to *N* is as follows:

$$\frac{\text{var}(\hat{N})}{N} \propto N^3 \qquad (d=1)$$

$$\frac{\text{var}(\hat{N})}{N} \propto N \qquad (d=2)$$

$$\frac{\text{var}(\hat{N})}{N} \propto N^{1/3} \qquad (d=3)$$

$$\frac{\text{var}(\hat{N})}{N} \propto \ln N \qquad (d=4) \,. \tag{47}$$

This shows that for all dimensions below and including four, particle fluctuations are anomalous, corresponding to an unstable systems. In that sense, the dimension four is critical, implying that the stability condition for a condensed gas in a box is:

$$d\_{\mathbb{C}} > d\_{\mathbb{C}} = 4 \qquad \left( T < T\_{\mathbb{C}} \right) \,. \tag{48}$$

#### **4. Ideal Gas in a Power-Law Trap**

Power-law traps are the most often used devices for trapping particles. Here, we study particle fluctuations and the related stability of mesoscopic clouds in such traps.

#### *4.1. Modified Semiclassical Approximation*

The general form of confining potentials, employed in power-law traps, can be represented as:

$$\mathcal{U}(\mathbf{r}) = \sum\_{a=1}^{d} \frac{\omega\_a}{2} \left| \begin{array}{c} r\_a \\ I\_a \end{array} \right|^{n\_a} \,, \tag{49}$$

where:

$$l\_a \equiv \frac{1}{\sqrt{m\omega\_a}}$$

is the effective trap radius in the *α* direction. As a whole, a trap can be characterized by the effective trap frequency *ω*<sup>0</sup> and effective length *l*<sup>0</sup> connected by the relations:

$$
\omega\_0 \equiv \left(\prod\_{a=1}^d \omega\_a\right)^{1/d} = \frac{1}{m l\_0^2} \; \; \quad l\_0 \equiv \left(\prod\_{a=1}^d l\_a\right)^{1/d} = \frac{1}{\sqrt{m \omega\_0}} \; \; \tag{50}
$$

In the limit *n<sup>α</sup>* → ∞, we return to a rectangular box.

When the effective trap frequency is much lower than temperature,

$$\frac{\omega\_0}{T} \ll 1\,\tag{51}$$

it is possible to resort to the semiclassical approximation that, however, needs to be modified for considering mesoscopic systems [18,31].

In the semiclassical approximation, one defines the density of states:

$$\rho(\varepsilon) = \frac{(2m)^{d/2}}{(4\pi)^{d/2}\Gamma(d/2)} \int\_{\mathbb{V}\_{\varepsilon}} [\varepsilon - \mathcal{U}(\mathbf{r})]^{d/2 - 1} \, d\mathbf{r} \,\,\mu$$

in which:

$$\mathbb{V}\_{\varepsilon} \equiv \{ \mathbf{r} : \ U(\mathbf{r}) \le \varepsilon \}$$

is the volume available for particle motion.

For trapped particles, an important notion is the confining dimension [18,31]:

$$D \equiv d + \sum\_{\alpha=1}^{d} \frac{2}{n\_{\alpha}} \,. \tag{52}$$

The density of states for the power-law potential (49) reduces to:

$$\rho(\varepsilon) = \frac{\varepsilon^{D/2 - 1}}{\gamma\_D \Gamma(D/2)} \; , \tag{53}$$

where we use the notation:

$$\gamma\_D \equiv \frac{\pi^{d/2}}{2^{D/2}} \prod\_{\alpha=1}^d \frac{\omega\_\alpha^{1/2 + 1/n\_\alpha}}{\Gamma(1 + 1/n\_\alpha)}$$

In the normal state above *Tc*, the number of particles is given by the formula:

$$N = \frac{T^{D/2}}{\gamma\_D} \ g\_{D/2}(z) \qquad \left(T \ge T\_c\right). \tag{54}$$

.

We again meet the Bose function that has to be modified according to the definition (19) by using the integral cutoff:

$$
\mu\_0 = \frac{\varepsilon\_0}{T} \qquad \left(\varepsilon\_0 \sim \omega\_0\right) \tag{55}
$$

with *ε*<sup>0</sup> being the lowest energy level in the trap, which is of the order of *ω*0.

#### *4.2. Condensation Temperature of a Gas in a Power-Law Trap*

At the critical temperature *Tc*, we have *µ* = 0 and *z* = 1. Then, Equation (54) yields:

$$T\_c = \left[\frac{\gamma\_D N}{g\_{D/2}(1)}\right]^{2/D} \,. \tag{56}$$

The modified Bose function, depending on the confining dimension, takes the forms:

$$g\_{D/2}(1) = \frac{2}{(2-D)\Gamma(D/2)} \left(\frac{T}{\varepsilon\_0}\right)^{1-D/2} \qquad (D < 2) \,\,\mu$$

$$g\_1(1) = \ln \,\, \frac{T}{\varepsilon\_0} \qquad (D = 2) \,\,\,\,$$

$$g\_{D/2}(1) = \zeta\left(\frac{D}{2}\right) \qquad (D > 2) \,\,.$$

If *D* < 2, the spatial dimension can only be *d* = 1, when:

$$
\gamma\_D = \frac{\sqrt{\pi}}{\Gamma(D)} \left(\frac{\omega\_0}{2}\right)^{D/2} \qquad (d=1)\dots
$$

Then, the critical temperature is:

$$T\_c = \frac{\sqrt{\pi}}{\Gamma(D)} \left(1 - \frac{D}{2}\right) \Gamma\left(\frac{D}{2}\right) \left(\frac{\omega\_0}{2\varepsilon\_0}\right)^{D/2} N\varepsilon\_0 \qquad (D < 2, \, d = 1) \,. \tag{57}$$

The confining dimension equals two, *D* = 2, when *d* = 1 and *n* = 2, so that:

$$
\gamma\_2 = \omega\_0 \qquad (D = 2 \text{ , } d = 1 \text{ , } n = 2) \text{ .}
$$

This yields the critical temperature:

$$T\_c = \frac{N\omega\_0}{\ln(T\_c/\varepsilon\_0)} \qquad \left(D = 2, \ d = 1\right). \tag{58}$$

For large *N*, one has:

$$\frac{T\_c}{\varepsilon\_0} \ll \exp\left(\frac{\omega\_0}{\varepsilon\_0} \, N\right) \,\prime$$

since for a one-dimensional harmonic oscillator, *ε*<sup>0</sup> = *ω*0/2. Because of this, the critical temperature (58) can be simplified to:

$$T\_{\mathcal{L}} = \frac{N\omega\_0}{\ln(2N)}\,. \tag{59}$$

For the confining dimension larger than two, the critical temperature is:

$$T\_c = \left[\frac{\gamma\_D N}{\zeta(D/2)}\right]^{2/D} \qquad (D > 2) \,. \tag{60}$$

In the case of harmonic traps, when *n<sup>α</sup>* = 2, hence *D* = 2*d* and *γ<sup>D</sup>* = *ω<sup>d</sup>* 0 , the critical temperature becomes:

$$T\_{\mathfrak{c}} = \left[\frac{N}{\zeta(d)}\right]^{1/d} \omega\_0 \qquad (n\_a = 2, \ d > 1) \dots$$

#### *4.3. Scaling with Respect to the Particle Number*

As is explained in Section 2, extensive observables are proportional to the number of particles *N*, when this number is large. This definition prescribes the scaling of the system characteristics with respect to *N*. As a representative of an observable quantity, we may take, e.g., internal energy:

$$
\langle \hat{H} \rangle = \langle H \rangle + \mu \mathcal{N} \,. \tag{61}
$$

This is an extensive quantity satisfying the condition:

$$\frac{\langle \hat{H} \rangle}{N} \simeq \text{const} \qquad \text{(N \gg 1)}\,. \tag{62}$$

For the considered case of a gas in a power-law trap, we have:

$$\frac{\langle \hat{H} \rangle}{N} = \frac{D \lg\_{1 + D/2}(z)}{2N \gamma\_D} \ T^{1 + D/2} \,. \tag{63}$$

The function *g*1+*D*/2(*z*) is finite for all *D* > 0 and all *z*. Hence, the condition (62) implies:

$$N\gamma\_D \simeq \text{const} \qquad (N \gg 1); \tag{64}$$

To make the consideration slightly less cumbersome, let us set the powers *n<sup>α</sup>* = *n* for the trapping potential. Then, the confining dimension is:

$$D = \left(1 + \frac{2}{n}\right)d\,. \tag{65}$$

*γ<sup>D</sup>* becomes:

$$\gamma\_D = \frac{\pi^{d/2}}{\Gamma^d (1 + 1/n)} \left(\frac{\omega\_0}{2}\right)^{D/2}$$

which tells us that:

$$
\gamma\_D \propto \omega\_0^{D/2} \qquad (N \gg 1) \; .
$$

Therefore, *ω*<sup>0</sup> scales as:

$$
\omega\_0 \propto \frac{1}{N^{2/D}} \qquad \left(N \gg 1\right). \tag{66}
$$

,

Using this scaling and the fact that *ω*<sup>0</sup> ∼ *ε*0, we see that the critical temperatures from the previous subsection behave as:

$$T\_{\mathbb{C}} \propto \frac{1}{N^{2/D - 1}} \qquad (D < 2)$$

$$T\_{\mathbb{C}} \propto \frac{1}{\ln N} \qquad (D = 2)$$

$$T\_{\mathbb{C}} \propto \text{const} \qquad (D > 2) \,. \tag{67}$$

#### *4.4. Fluctuations above the Condensation Temperature*

Particle fluctuations above the condensation temperature are described by the formula:

$$\frac{\text{var}(\hat{N})}{N} = \frac{T^{D/2}}{N\gamma\_D} z \frac{\partial \mathbf{g}\_{D/2}(z)}{\partial z} \qquad (T > T\_c) \tag{68}$$

where *z* < 1. For the modified Bose function, we have:

$$\frac{\partial g\_m(z)}{\partial z} = \frac{1}{z} \, g\_{m-1}(z) + \frac{1}{(1-z)\Gamma(m)} \left( u\_0^{m-1} - \frac{u\_0^m}{1-z} \right) \qquad (z < 1) \,\tag{69}$$

with the value:

$$g\_{m-1}(z) = -\frac{z}{(1-z)\Gamma(m)} \left[ u\_0^{m-1} - \frac{m-1}{m(1-z)} u\_0^m \right] \qquad (m < 1, \ z < 1) \tag{70}$$

for *m* < 1. Summarizing, we have the derivatives:

$$\frac{\partial g\_m(z)}{\partial z} = -\frac{u\_0^m}{(1-z)^2 \Gamma(1+m)} \qquad (m < 1 \; , z < 1)$$

$$\frac{\partial g\_1(z)}{\partial z} = -\frac{u\_0}{(1-z)^2} \qquad (m = 1 \; , z < 1)$$

$$\frac{\partial g\_m(z)}{\partial z} = \frac{1}{z} \; g\_{m-1}(z) \qquad (m > 1 \; , z < 1) \; .$$

From here, we find the relative variance:

$$\frac{\text{var}(\hat{\mathbf{N}})}{N} = -\frac{zT^{D/2}}{(1-z)^2 N \gamma\_D \Gamma(1+D/2)} \left(\frac{\varepsilon\_0}{T}\right)^{D/2} \qquad (D < 2, \ T > T\_c),$$

$$\frac{\text{var}(\hat{\mathbf{N}})}{N} = -\frac{zT}{(1-z)^2 N \gamma\_2} \left(\frac{\varepsilon\_0}{T}\right) \qquad (D = 2, \ T > T\_c),$$

$$\frac{\text{var}(\hat{\mathbf{N}})}{N} = \frac{T^{D/2}}{N \gamma\_D} \ g\_{D/2-1}(z) \qquad (D > 2, \ T > T\_c) \,, \tag{71}$$

characterizing particle fluctuations above the critical temperature. For *D* ≤ 2, the variance is negative, which means instability. The system is stable only for *D* > 2, giving the stability condition:

$$d + \sum\_{a=1}^{d} \frac{2}{n\_a} > 2 \qquad \left( T > T\_{\mathcal{E}} \right) \,. \tag{72}$$

#### *4.5. Fluctuations below the Condensation Temperature*

Below the condensation temperature, where *µ* = 0 and *z* = 1, the number of uncondensed particles reads as:

$$N\_1 = \frac{T^{D/2}}{\gamma\_D} \ g\_{D/2}(1) \qquad \left(T \le T\_c\right). \tag{73}$$

The variance of the total number of particles coincides with that of the uncondensed particles, which leads to:

$$\frac{\text{var}(\hat{N})}{N} = \frac{T^{D/2}}{N\gamma\_D} \lim\_{z \to 1} \frac{\partial \mathbf{g}\_{D/2}(z)}{\partial z} \qquad (T < T\_c) \,. \tag{74}$$

For the derivative in the right-hand side of the above formula, we have:

$$\lim\_{z \to 1} \frac{\partial g\_{D/2}(z)}{\partial z} = g\_{D/2 - 1}(1) + \frac{u\_0^{D/2 - 2}}{\Gamma(D/2)}.\tag{75}$$

Employing the values:

$$g\_{D/2-1}(1) = \frac{2}{(4-D)\Gamma(D/2-1)} \ u\_0^{D/2-2} \qquad (D < 4) \; \therefore$$

*Symmetry* **2019**, *11*, 603

$$g\_1(1) = -\ln u\_0 \qquad (D = 4) \; .$$

we get the derivatives:

$$\lim\_{z \to 1} \frac{\partial g\_{D/2}(z)}{\partial z} = \left[ \frac{2}{(4 - D)\Gamma(D/2 - 1)} + \frac{1}{\Gamma(D/2)} \right] u\_0^{D/2 - 2} \qquad (D < 4),$$

$$\lim\_{z \to 1} \frac{\partial g\_2(z)}{\partial z} = -\ln u\_0 \qquad (D = 4),$$

$$\lim\_{z \to 1} \frac{\partial g\_{D/2}(z)}{\partial z} = \mathcal{g}\_{D/2 - 1}(1) = \zeta \left( \frac{D}{2} - 1 \right) \qquad (D > 4).$$

In that way, we come to the relative variances:

$$\frac{\text{var}(\hat{N})}{N} = \frac{T^{D/2}}{N\gamma\_{D}} \left[ \frac{2}{(4-D)\Gamma(D/2-1)} + \frac{1}{\Gamma(D/2)} \right] \left( \frac{T}{\varepsilon\_{0}} \right)^{2-D/2} \qquad (D < 4) \text{ .}$$

$$\frac{\text{var}(\hat{N})}{N} = \frac{T^{2}}{N\gamma\_{4}} \ln \left( \frac{T}{\varepsilon\_{0}} \right) \qquad (D = 4) \text{ .}$$

$$\frac{\text{var}(\hat{N})}{N} = \frac{T^{D/2}}{N\gamma\_{D}} \zeta \left( \frac{D}{2} - 1 \right) \qquad (D > 4) \text{ .} \tag{76}$$

Keeping in mind that *ε*<sup>0</sup> ∝ *N*−2/*<sup>D</sup>* results in the scaling:

$$
\frac{\text{var}(\hat{N})}{N} \propto N^{(4-D)/D} \qquad (D < 4)
$$

$$
\frac{\text{var}(\hat{N})}{N} \propto \ln N \qquad (D = 4)
$$

$$
\frac{\text{var}(\hat{N})}{N} \propto \text{const} \qquad (D > 4)
\text{.}\tag{77}
$$

This tells us that the system is stable only for *D* > 4. Therefore, the stability condition is:

$$d + \sum\_{a=1}^{d} \frac{2}{n\_a} > 4 \qquad \left( T < T\_\varepsilon \right). \tag{78}$$

Notice that in the case of the often considered harmonic potential, when *n<sup>α</sup>* = 2, we have *D* = 2*d* and *γ<sup>D</sup>* = *ω<sup>d</sup>* 0 . Then, the stability condition (78) reduces to the condition *d* > 2. The relative particle variance reads as:

$$\frac{\text{var}(\hat{N})}{N} = \frac{\zeta(d-1)}{\zeta(d)} \left(\frac{T}{T\_{\varepsilon}}\right)^{d} \qquad (n\_{a} = 2, \ d > 2) \text{ .}$$

#### **5. Interacting Bose System above the Condensation Temperature**

The grand Hamiltonian for a system of interacting Bose particles is:

$$H = \int \psi^\dagger(\mathbf{r}) \left( -\frac{\nabla^2}{2m} - \mu \right) \psi(\mathbf{r}) \, d\mathbf{r} +$$

$$+ \frac{1}{2} \int \psi^\dagger(\mathbf{r}) \psi^\dagger(\mathbf{r}') \Phi(\mathbf{r} - \mathbf{r}') \psi(\mathbf{r}') \psi(\mathbf{r}) \, d\mathbf{r} d\mathbf{r}' \,. \tag{79}$$

For generality, we consider a nonlocal isotropic interaction potential Φ(**r**) = Φ(*r*), where *r* ≡ |**r**|. The integration is assumed to be over a rectangular box of volume *V* confining the system.

In the Hartree–Fock approximation, the Hamiltonian takes the form:

$$H\_{HF} = E\_{HF} + \int \psi^\dagger(\mathbf{r}) \left( -\frac{\nabla^2}{2m} - \mu \right) \psi(\mathbf{r}) \, d\mathbf{r} +$$

$$+ \int \Phi(\mathbf{r} - \mathbf{r}') \left[ \rho(\mathbf{r}') \psi^\dagger(\mathbf{r}) \psi(\mathbf{r}) + \rho(\mathbf{r}', \mathbf{r}) \psi^\dagger(\mathbf{r}') \psi(\mathbf{r}) \right] \, d\mathbf{r} d\mathbf{r}',\tag{80}$$

where:

$$E\_{HF} = -\frac{1}{2} \int \Phi(\mathbf{r} - \mathbf{r}') \left[ \rho(\mathbf{r}) \rho(\mathbf{r}') + |\rho(\mathbf{r}, \mathbf{r}')|^2 \right] \, d\mathbf{r} d\mathbf{r}' $$

and the notations are used for the single-particle density matrix:

$$
\rho(\mathbf{r}, \mathbf{r}') = \langle \psi^\dagger(\mathbf{r}')\psi(\mathbf{r})\rangle \tag{81}
$$

and the particle density:

$$
\rho(\mathbf{r}) = \rho(\mathbf{r}, \mathbf{r}) = \left< \psi^\dagger(\mathbf{r}) \psi(\mathbf{r}) \right> . \tag{82}
$$

Employing the expansion of the field operators over plane waves, as in Equation (13), we get the Hamiltonian:

$$H\_{\rm HF} = E\_{\rm HF} + \sum\_{k} \omega\_{k} a\_{k}^{\dagger} a\_{k} \, , \tag{83}$$

in which:

$$E\_{HF} = -\frac{1}{2} \rho \Phi\_0 N - \frac{1}{2V} \sum\_{kp} n\_k n\_p \Phi\_{k+p} \dots$$

Φ*k* is a Fourier transform of Φ(**r**), and:

$$
\Phi\_0 = \int \Phi(\mathbf{r}) \, d\mathbf{r} \,. \tag{84}
$$

The momentum distribution is given by the expression (16), with the spectrum:

$$
\omega\_k = \frac{k^2}{2m} + \rho \Phi\_0 + \frac{1}{V} \sum\_p n\_p \Phi\_{k+p} - \mu \,. \tag{85}
$$

The function *n<sup>p</sup>* possesses a maximum at *p* → 0, because of which it is possible to use the approximation [32,33]:

$$\sum\_{p} n\_{p} \Phi\_{k+p} \cong \Phi\_{k} \sum\_{p} n\_{p} \tag{86}$$

giving:

$$
\omega\_k = \frac{k^2}{2m} + \rho(\Phi\_0 + \Phi\_k) - \mu \,. \tag{87}
$$

Introducing the effective interaction radius by the relation:

$$r\_{eff}^2 \equiv \frac{\int \Phi(\mathbf{r}) r^2 d\mathbf{r}}{\int \Phi(\mathbf{r}) d\mathbf{r}} = \frac{4\pi}{\Phi\_0} \int\_0^\infty \Phi(\mathbf{r}) r^4 \,d\mathbf{r} \tag{88}$$

shows that the long-wave limit of Φ*<sup>k</sup>* is:

$$
\Phi\_k \cong \left(1 - \frac{1}{6} \, k^2 r\_{eff}^2\right) \Phi\_0 \,. \tag{89}
$$

Then, the spectrum (87) can be represented as:

$$
\omega\_k \simeq \frac{k^2}{2m^\*} - \mu\_{eff} \qquad (k \to 0) \,\,\,\,\tag{90}
$$

*Symmetry* **2019**, *11*, 603

with the effective mass:

$$m^\* \equiv \frac{m}{1 - \rho \Phi\_0 r\_{eff}^2 / 3} \tag{91}$$

and effective chemical potential:

$$
\mu\_{eff} \equiv \mu - 2\rho \Phi\_0 \,. \tag{92}
$$

In this approximation, the number of particles acquires the same form (23), however with the notations:

$$
\lambda\_T \equiv \sqrt{\frac{2\pi}{m^\*T}}\,,\qquad z \equiv \exp(\beta \mu\_{eff})\,.\tag{93}
$$

Using again the modified Bose function (19) and following the same analysis as in Section 3, we come to the conclusion that the system is stable for *d* > 2, when *T* > *Tc*. The difference is that now instead of mass *m*, there is the effective mass *m*∗ , and at the critical temperature, we have:

$$
\mu\_{eff} = 0 \,, \qquad \mu = 2\rho\Phi\_0 \qquad (T = T\_c) \, .
$$

Therefore, the critical temperature becomes:

$$T\_c = \frac{2\pi}{m^\*} \left[\frac{\rho}{g\_{d/2}(1)}\right]^{2/d}.\tag{94}$$

As an example, let us consider the realistic three-dimensional case. Using the Robinson representation (see the details in review [18]), we can find the behavior of the effective chemical potential at high temperatures:

$$
\mu\_{eff} = T \ln \left( \rho \lambda\_T^3 \right) \qquad \left( T \gg T\_\mathfrak{c} \right) \tag{95}
$$

and at the temperature approaching the critical point from above,

$$
\mu\_{eff} \simeq -T \frac{\zeta^2(3/2)}{4\pi} \left[ 1 - \left( \frac{T\_c}{T} \right)^{3/2} \right]^2. \tag{96}
$$

Then, the isothermal compressibility:

$$\kappa\_T = \frac{g\_{1/2}(z)}{\rho^2 T \lambda\_T^3} \tag{97}$$

at high temperatures is:

$$
\kappa\_T \simeq \frac{1}{\rho T} \qquad \left(T \gg T\_\mathbb{C}\right) \,, \tag{98}
$$

while close to the critical point, it is:

$$
\kappa\_T \simeq \frac{0.921}{\rho T} \left[ 1 - \left( \frac{T\_c}{T} \right)^{3/2} \right]^{-1}; \tag{99}
$$

respectively, particle fluctuations, described by the relative variance:

$$\frac{\text{var}(\hat{N})}{N} = \rho T \kappa\_T = \frac{g\_{1/2}(z)}{\rho \lambda\_T^3} \, , \tag{100}$$

at high temperatures behave as:

$$\frac{\text{var}(\hat{N})}{N} \simeq 1 \qquad (T \gg T\_{\mathfrak{c}}),\tag{101}$$

and close to the critical point, we get:

$$\frac{\text{var}(\hat{N})}{N} \simeq 0.921 \left[ 1 - \left( \frac{T\_c}{T} \right)^{3/2} \right]^{-1}. \tag{102}$$

Outside of the critical temperature itself, particle fluctuations are thermodynamically normal. The divergence of the compressibility at the critical point signifies a second-order phase transition. At the point of the phase transition, the system is not stable, and the fluctuations do not need to be finite.

#### **6. Interacting Bose System below the Condensation Temperature**

In Section 3, it is proven that the ideal Bose gas, confined in a box, is stable below the condensation temperature only for *d* > 4. In the present section, we show that interactions stabilize the system, making it stable already for *d* = 3.

#### *6.1. Self-Consistent Approach*

For describing a Bose system with the Bose–Einstein condensate, we employ the self-consistent approach [16–18,24,34,35], providing a gapless spectrum, correct thermodynamics, the validity of all conservation laws, and good agreement with Monte Carlo simulations and experimental data.

The energy Hamiltonian has the form:

$$
\hat{H} = \int \hat{\psi}^{\dagger}(\mathbf{r}) \left( -\frac{\nabla^2}{2m} \right) \hat{\psi}(\mathbf{r}) \, d\mathbf{r} +
$$

$$
+\, \frac{1}{2} \int \hat{\psi}^{\dagger}(\mathbf{r}) \hat{\psi}^{\dagger}(\mathbf{r}') \Phi(\mathbf{r} - \mathbf{r}') \hat{\psi}(\mathbf{r}') \hat{\psi}(\mathbf{r}) \, d\mathbf{r} d\mathbf{r}' \,. \tag{103}
$$

The genuine Bose–Einstein condensation necessarily requires global gauge symmetry breaking [6,9,17,18]. Finite systems, strictly speaking, do not exhibit this symmetry breaking. However, a system with a large number of particles *N* 1 enjoys asymptotic symmetry breaking [36] in the sense that the system properties asymptotically, with respect to *N*, are close to the system with broken symmetry. The global gauge symmetry can be broken by the Bogolubov shift [28–30]:

$$
\hat{\psi}(\mathbf{r}) = \eta(\mathbf{r}) + \psi\_1(\mathbf{r}) \,. \tag{104}
$$

in which the condensate function *η*(**r**) and the operator of uncondensed particles *ψ*1(**r**) are mutually orthogonal,

$$
\int \eta^\*(\mathbf{r}) \psi\_1(\mathbf{r}) \, d\mathbf{r} = 0 \tag{105}
$$

and the operator of uncondensed particles satisfies the condition:

$$
\langle \psi\_1(\mathbf{r}) \rangle = \mathbf{0} \,. \tag{106}
$$

The number of condensed particles is:

$$N\_0 = \int |\eta(\mathbf{r})|^2 \, d\mathbf{r} \,\tag{107}$$

while the number of uncondensed particles is given by the average:

$$N\_1 = \langle \hat{N}\_1 \rangle \,, \qquad \hat{N}\_1 = \int \psi\_1^\dagger(\mathbf{r}) \psi\_1(\mathbf{r}) \, d\mathbf{r} \, . \tag{108}$$

*Symmetry* **2019**, *11*, 603

The grand Hamiltonian reads as:

$$H = \hat{H} - \mu\_0 \mathbf{N}\_0 - \mu\_1 \hat{\mathbf{N}}\_1 - \hat{\mathbf{A}} \tag{109}$$

where:

$$
\hat{\Lambda} = \int \left[ \lambda(\mathbf{r}) \psi\_1^\dagger(\mathbf{r}) + \lambda^\*(\mathbf{r}) \psi\_1(\mathbf{r}) \right] \, d\mathbf{r}
$$

and *µ*0, *µ*1, and *λ*(**r**) are Lagrange multipliers guaranteeing the validity of the normalizations (107) and (108), as well as condition (106).

The evolution equation for the condensate function can be written as:

$$i\hbar\frac{\partial}{\partial t}\,\eta(\mathbf{r},t) = \left\langle \frac{\delta H}{\delta \eta^\*(\mathbf{r},t)} \right\rangle \tag{110}$$

and the equation for the operator of uncondensed particles as:

$$i\hbar\frac{\partial}{\partial t}\,\psi\_1(\mathbf{r},t) = \frac{\delta H}{\delta\psi\_1^\dagger(\mathbf{r},t)}\,\,. \tag{111}$$

Keeping in mind, as usual, the periodic continuation of the box, we expand the field operators in plane waves, as in (13), and assume the existence of the Fourier representation for the interaction potential:

$$\Phi\_k = \int \Phi(\mathbf{r}) e^{-i\mathbf{k}\cdot\mathbf{r}} \, d\mathbf{r} \,, \qquad \Phi(\mathbf{r}) = \frac{1}{V} \sum\_k \Phi\_k e^{i\mathbf{k}\cdot\mathbf{r}} \,. \tag{112}$$

Then, we get the normal density matrix:

$$\rho\_1(\mathbf{r}, \mathbf{r}') = \langle \psi\_1^\dagger(\mathbf{r}')\psi\_1(\mathbf{r})\rangle = \frac{1}{V} \sum\_{k \neq 0} n\_k e^{i\mathbf{k} \cdot (\mathbf{r} - \mathbf{r}')} \tag{113}$$

and the anomalous matrix:

$$
\sigma\_1(\mathbf{r}, \mathbf{r}') = \langle \psi\_1(\mathbf{r}')\psi\_1(\mathbf{r})\rangle = \frac{1}{V} \sum\_{k \neq 0} \sigma\_k e^{i\mathbf{k} \cdot (\mathbf{r} - \mathbf{r}')} \,\tag{114}
$$

in which:

$$m\_k \equiv \langle a\_k^\dagger a\_k \rangle \; , \qquad \sigma\_k \equiv \langle a\_k a\_{-k} \rangle \; . \tag{115}$$

The condensate function *η*(**r**) = *η* defines the condensate density:

$$\rho\_0 \equiv \frac{N\_0}{V} = |\eta|^2. \tag{116}$$

The density of uncondensed particles is:

$$\rho\_1 \equiv \frac{N\_1}{V} = \rho\_1(\mathbf{r}, \mathbf{r}) = \frac{1}{V} \sum\_k n\_k \,. \tag{117}$$

The diagonal anomalous matrix gives the anomalous average:

$$
\sigma\_1 \equiv \sigma\_1(\mathbf{r}, \mathbf{r}) = \frac{1}{V} \sum\_k \sigma\_k \,. \tag{118}
$$

The average density of particles is the sum:

$$
\rho \equiv \frac{N}{V} = \rho\_0 + \rho\_1 \,. \tag{119}
$$

*Symmetry* **2019**, *11*, 603

Then, we use the Hartree–Fock–Bogolubov approximation and accomplish the Bogolubov canonical transformation:

$$a\_k = u\_k b\_k + v\_{-k}^\* b\_{-k}^\dagger \qquad \quad b\_k = u\_k^\* a\_k - v\_k^\* a\_{-k}^\dagger \,\_2$$

where *u<sup>k</sup>* and *v<sup>k</sup>* are chosen so as to diagonalize the Hamiltonian. In that way, we obtain the diagonalized Hamiltonian:

$$H\_B = E\_B + \sum\_k \varepsilon\_k b\_k^\dagger b\_k \tag{120}$$

in which:

$$E\_B = -\frac{1}{2}N\rho\Phi\_0 - \rho\_0\sum\_p (n\_p + \sigma\_p)\Phi\_p - \frac{1}{2V}\sum\_{kp} (n\_k n\_p + \sigma\_k \sigma\_p)\Phi\_{k+p} + \frac{1}{2}\sum\_k (\varepsilon\_k - \omega\_k)\Phi\_k$$

the particle spectrum is:

$$
\varepsilon\_k = \sqrt{\omega\_k^2 - \Delta\_k^2} \tag{121}
$$

and where:

$$
\omega\_k = \frac{k^2}{2m} + \Delta + \rho\_0(\Phi\_k - \Phi\_0) + \frac{1}{V} \sum\_p n\_p (\Phi\_{k+p} - \Phi\_p) \tag{12}
$$

$$
\Delta\_k = \rho\_0 \Phi\_k + \frac{1}{V} \sum\_p \sigma\_p \Phi\_{k+p} \, \prime \qquad \Delta \equiv \lim\_{k \to 0} \Delta\_k = \rho\_0 \Phi\_0 + \frac{1}{V} \sum\_p \sigma\_p \Phi\_p \, \, \, \tag{122}
$$

For the expressions in (115), we find:

$$
\sigma\_k = \frac{\omega\_k}{2\varepsilon\_k} \coth\left(\frac{\varepsilon\_k}{2T}\right) - \frac{1}{2}\,'\,, \qquad \sigma\_k = -\frac{\Delta\_k}{2\varepsilon\_k} \coth\left(\frac{\varepsilon\_k}{2T}\right)\,. \tag{123}
$$

The chemical potentials are:

$$
\mu\_0 = \rho \Phi\_0 + \frac{1}{V} \sum\_k (n\_k + \sigma\_k) \Phi\_k \ , \qquad \mu\_1 = \rho \Phi\_0 + \frac{1}{V} \sum\_k (n\_k - \sigma\_k) \Phi\_k \ . \tag{124}
$$

In the long-wave limit, we can use the expansion:

$$
\Phi\_{k+p} \simeq \Phi\_p + \frac{k^2}{2} \Phi\_p'' \qquad (k \to 0) \ \lambda
$$

where:

$$
\Phi\_p'' \equiv \frac{\partial^2 \Phi\_p}{\partial p^2} \dots
$$

Then, the spectrum (121) becomes of the phonon type:

$$
\varepsilon\_k \simeq ck \qquad (k \to 0) \tag{125}
$$

with the sound velocity:

$$c = \sqrt{\frac{\Delta}{m\_{eff}}}\tag{126}$$

and with the notation for the effective mass:

$$m\_{eff} \equiv \frac{m}{1 + \frac{m}{V} \sum\_{p} (n\_p - \sigma\_p) \Phi\_p^{\prime\prime}} \,. \tag{127}$$

Actually, Expression (126), which can be written as:

$$
\boldsymbol{m}\_{eff}\boldsymbol{c}^2 = \boldsymbol{\Delta}\_{\boldsymbol{\nu}}
$$

is the equation:

$$\frac{mc^2}{1 + \frac{m}{V} \sum\_{p} (n\_p - \sigma\_p) \Phi\_p''} = \rho\_0 \Phi\_0 + \frac{1}{V} \sum\_{p} \sigma\_p \Phi\_p \tag{128}$$

defining the sound velocity *c*.

To simplify the consideration, we can resort to the approximation (86), similarly, to which we can write:

$$\sum\_{p} \sigma\_{p} \Phi\_{k+p} \cong \Phi\_{k} \sum\_{p} \sigma\_{p} \,. \tag{129}$$

This gives:

$$\frac{1}{V} \sum\_{p} (n\_p - \sigma\_p) \Phi\_p^{\prime\prime} = (\rho\_1 - \sigma\_1) \Phi\_0^{\prime\prime} \ .$$

where:

$$\Phi\_0^{\prime\prime} = \lim\_{p \to 0} \Phi\_p^{\prime\prime} = -\frac{4\pi}{3} \int\_0^\infty \Phi(\mathbf{r}) r^4 \,d\mathbf{r} \,\dots$$

In view of the notation for the effective interaction radius (88), we get:

$$
\Phi\_0'' = -\frac{1}{3} \left. \Phi\_0 r\_{eff}^2 \right| \,\mathrm{s}
$$

Then, the effective mass (127) acquires the form:

$$m\_{eff} = \frac{m}{1 + (\sigma\_1 - \rho\_1)\Phi\_0 m r\_{eff}^2 / 3} \,. \tag{130}$$

In the approximations (86) and (129), the chemical potentials (124) become:

$$
\mu\_0 = \rho \Phi\_0 + (\rho\_1 + \sigma\_1) \Phi\_0 \ , \qquad \mu\_1 = \rho \Phi\_0 + (\rho\_1 - \sigma\_1) \Phi\_0 \ . \tag{131}
$$

Furthermore, we have:

$$
\omega\_k = \frac{k^2}{2m} + \Delta + \rho(\Phi\_k - \Phi\_0) \qquad \Delta\_k = (\rho\_0 + \sigma\_1)\Phi\_k \, \, \, \, \, \, \Delta = (\rho\_0 + \sigma\_1)\Phi\_0 \, \,. \tag{132}
$$

The spectrum (121) can be written as:

$$\varepsilon\_{k}^{2} = \left[\frac{k^{2}}{2m} + (\rho\_{1} - \sigma\_{1})(\Phi\_{k} - \Phi\_{0})\right] \left[\frac{k^{2}}{2m} + \rho(\Phi\_{k} - \Phi\_{0}) + (\rho\_{0} + \sigma\_{1})(\Phi\_{k} + \Phi\_{0})\right].\tag{133}$$

The density of uncondensed particles is:

$$\rho\_1 = \int \left[\frac{\omega\_k}{2\varepsilon\_k}\coth\left(\frac{\varepsilon\_k}{2T}\right) - \frac{1}{2}\right] \frac{d\mathbf{k}}{(2\pi)^3} \,. \tag{134}$$

The anomalous average (118) can be represented in the form:

$$
\sigma\_1 = -\int \frac{\Delta\_k}{2\varepsilon\_k} \frac{d\mathbf{k}}{(2\pi)^3} - \int \frac{\Delta\_k}{2\varepsilon\_k} \left[ \coth\left(\frac{\varepsilon\_k}{2T}\right) - 1 \right] \frac{d\mathbf{k}}{(2\pi)^3} \,. \tag{135}
$$

When the first term here diverges, which happens for the local interaction, we can use dimensional regularization [18].

#### *6.2. Particle Fluctuations*

The number-of-particles variance can be found by involving the formula:

$$\frac{\text{var}(\hat{N})}{N} = 1 + \rho \int [\mathbf{g}(\mathbf{r}) - 1] \, d\mathbf{r} \, \tag{136}$$

in which:

$$\log(\mathbf{r}\_{12}) = \frac{1}{\mathcal{S}^2} \left< \psi^\dagger(\mathbf{r}\_1) \psi^\dagger(\mathbf{r}\_2) \psi(\mathbf{r}\_2) \psi(\mathbf{r}\_1) \right> \tag{137}$$

is the pair correlation function, with **r**<sup>12</sup> ≡ **r**<sup>1</sup> − **r**2.

Accomplishing the Bogolubov shift (104), we use the Hartree–Fock–Bogolubov (HFB) decoupling for the expressions containing the operators *ψ*1. Since, mathematically, the HFB approximation is of second order with respect to the products of the operators *ψ*1, it is necessary to leave in the pair correlation function only the terms of second order with respect to these operators [3,16–18,23,24]. As a result, we obtain:

$$\int [g(\mathbf{r}) - 1] \, d\mathbf{r} = \frac{2}{\rho} \lim\_{k \to 0} \left( n\_k + \sigma\_k \right) \,. \tag{138}$$

In this way, for the relative variance, we find:

$$\frac{\text{var}(\hat{N})}{N} = 1 + 2 \lim\_{k \to 0} \left( n\_k + \sigma\_k \right) \,. \tag{139}$$

,

For small *k*, when *ε<sup>k</sup>* tends to zero, we have:

$$n\_k \simeq \frac{T\Delta\_k}{\varepsilon\_k^2} + \frac{\Delta\_k}{12T} + \frac{T}{2\Delta\_k} - \frac{1}{2} + \left(\frac{\Delta\_k}{3T} - \frac{T}{\Delta\_k} - \frac{\Delta\_k^3}{90T^3}\right) \frac{\varepsilon\_k^2}{8\Delta\_k^2},$$

$$\sigma\_k \simeq -\frac{T\Delta\_k}{\varepsilon\_k^2} - \frac{\Delta\_k}{12T} + \frac{\Delta\_k \varepsilon\_k^2}{720T^3} \qquad (\varepsilon\_k \to 0) \,. \tag{140}$$

Therefore:

$$\lim\_{k \to 0} \left( n\_k + \sigma\_k \right) = \frac{1}{2} \left( \frac{T}{\Delta} - 1 \right)$$

with:

$$
\Delta = m\_{eff}c^2 = (\rho\_0 + \sigma\_1)\Phi\_0 \dots
$$

Thus, we come to the expression:

$$\frac{\text{var}(\hat{N})}{N} = \frac{T}{m\_{eff}c^2} \; ; \tag{141}$$

respectively, the compressibility is:

$$\kappa\_T = \frac{\text{var}(\hat{N})}{N\rho T} = \frac{1}{\rho m\_{eff} c^2} \,. \tag{142}$$

Taking into account Formula (126) leads to the variance:

$$\frac{\text{var}(\hat{N})}{N} = \frac{T}{(\rho\_0 + \sigma\_1)\Phi\_0} \,. \tag{143}$$

Note that Expression (143) is valid at zero temperature, as well. This is easy to check considering the quantities (123) at zero temperature,

$$n\_k = \frac{\sqrt{\varepsilon\_k^2 + \Delta\_k^2}}{2\varepsilon\_k} - \frac{1}{2} \quad \quad \sigma\_k = -\frac{\Delta\_k}{2\varepsilon\_k} \qquad (T = 0) \; .$$

From here, in the long-wave limit, we have:

$$m\_k \simeq \frac{\Delta\_k}{2\varepsilon\_k} + \frac{\varepsilon\_k}{4\Delta\_k} - \frac{1}{2} \qquad (\varepsilon\_k \to 0 \text{ , } T = 0) \text{ .}$$

Hence:

$$\lim\_{k \to 0} \left( n\_k + \sigma\_k \right) = -\frac{1}{2} \qquad \left( T = 0 \right).$$

and:

$$\frac{\text{var}(\hat{\mathcal{N}})}{N} = 0 \qquad (T = 0) \text{ .}$$

The above result for the relative variance (143) can be generalized for nonuniform systems [37] by involving the local-density approximation, which yields:

$$\frac{\text{var}(\hat{N})}{N} = \frac{T}{N} \int \frac{\rho(\mathbf{r})}{\Delta(\mathbf{r})} \, d\mathbf{r} \,\tag{144}$$

where:

$$
\Delta(\mathbf{r}) = [\rho\_0(\mathbf{r}) + \sigma\_1(\mathbf{r})] \Phi\_0 \,. \tag{145}
$$

Particle fluctuations in a three-dimensional Bose-condensed system of interacting particles are thermodynamically normal in both cases, when particles are in a box or in a nonuniform external potential.

#### **7. Conclusions**

Particle fluctuations in Bose systems were studied. Investigating the behavior of these fluctuations is important because they are directly connected with isothermal compressibility and define the system stability with respect to pressure variations. Thermodynamically-anomalous fluctuations signify system instability; while thermodynamically-normal fluctuations mean that the equilibrium system is stable. The obtained results are as follows.

The ideal Bose gas confined in a rectangular box is stable, depending on the temperature, in spatial dimensions:

$$d > 2 \qquad \left(T > T\_c\right).$$

$$d > 4 \qquad \left(T < T\_c\right).$$

The stability of the ideal Bose gas in a power-law trap depends on the confining dimension:

$$D \equiv d + \sum\_{\alpha=1}^{d} \frac{2}{n\_{\alpha}} \dots$$

This gas is stable for the confining dimensions:

*D* > 2 (*T* > *Tc*) , *D* > 4 (*T* < *Tc*) .

Interactions stabilize Bose-condensed systems, so that an interacting system with Bose–Einstein condensate becomes stable at *d* = 3 for either a system in a box or in an external potential.

Nonlocal interactions with a stronger strength or with a larger interaction radius increase the effective mass, hence diminishing the condensation temperature.

**Funding:** This research received no external funding.

**Acknowledgments:** The author is grateful to E.P. Yukalova for useful discussions.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


© 2019 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
