Condensation and superfluidity of $SU(N)$ Bose gas

TL;DR

This paper investigates the properties of condensate and superfluid densities in an $SU(N)$ Bose gas, using large-$N$ expansion and Monte Carlo simulations to analyze thermal fluctuations and compare theoretical predictions with numerical results.

Contribution

It introduces a large-$N$ expansion approach to derive a general expression for superfluid density in multi-component Bose gases, incorporating thermal fluctuations.

Findings

Superfluid density formula reproduces Landau's at low temperatures.

Numerical results agree qualitatively with Monte Carlo simulations.

Thermal fluctuations significantly affect superfluid properties.

Abstract

We perform the comprehensive comparison of properties of the condensate and superfluid densities for the $N$-component three-dimensional Bose gas with the symmetric inter- and intraspecies short-range interaction between particles. In particular, based on the large-$N$ expansion approach for many-boson systems we obtain general expression for density of the superfluid component that at very low temperatures reproduce the well-know Landau's formula and non-trivially includes the thermal fluctuations in the finite-temperature region, and compare it to the condensate density calculated previously. The numerically evaluated temperature dependencies are in a qualitatively good agreement with the results of Monte Carlo simulations.

Equations

$$S=\int dx{\psi }_{\sigma }^{\ast }(x)\left\{{\partial }_{\tau }-\xi -i\phi (x)\right\}{\psi }_{\sigma }(x)$$

$$-\frac{N}{2g}\int dx{\phi }^2(x),$$

$$S_{\mathrm{eff}}=-\frac{NV{\phi }_0^2}{2Tg}+NV{n}_0\tilde{\mu }/T$$

$$-N{n}_0\int dx\tilde{\phi }(x){\left\{\tilde{\xi }+i\tilde{\phi }(x)-{\partial }_{\tau }\right\}}^{-1}\tilde{\phi }(x)$$

$$-\frac{N}{2g}\int dx{\tilde{\phi }}^2(x)-N\mathrm{Sp}\ln \left\{\tilde{\xi }+i\tilde{\phi }(x)-{\partial }_{\tau }\right\},$$

$$\frac{\Omega }{N}=\frac{V{\phi }_0^2}{2g}-V{n}_0\tilde{\mu }-T\sum _{\mathbf{k}}\ln \left(1-e^{-{\tilde{\xi }}_k/T}\right)$$

$$+\frac{T}{2N}\sum _K\ln \left[1+g\Pi (K)\right],$$

$${\Pi }_0(K)=\frac{n_0}{{\tilde{\xi }}_k-i{\omega }_k}+c.c.,$$

$${\Pi }_T(K)=\frac{T}{V}\sum _{K^{\prime }}\frac{1}{{\tilde{\xi }}_{k^{\prime }}-i{\omega }_{k^{\prime }}}\frac{1}{{\tilde{\xi }}_{|{\mathbf{k}}^{\prime }+\mathbf{k}|}-i{\omega }_{k^{\prime }+k}},$$

$$\mu =ng+\frac{T}{NV}\sum _K{g}_K\frac{{\tilde{\xi }}_k+i{\omega }_k}{E_K^2+{\omega }_k^2},$$

$$g=\frac{4\pi {\hslash }^{2}a}{m}+{\left(\frac{4\pi {\hslash }^{2}a}{m}\right)}^2\frac{1}{NV}\sum _{\mathbf{k}}\frac{1}{2{\epsilon }_k}\pm \dots$$

$$n_0=n-\left(1-\frac{1}{N}\right)\frac{1}{V}\sum _{\mathbf{k}}n({\tilde{\xi }}_k/T)-\frac{1}{NV}\sum _{\mathbf{k}}\frac{\partial n({\tilde{\xi }}_k/T)}{\partial {\tilde{\xi }}_k}$$

$$\times {\Sigma }_R^{(1)}({\tilde{\xi }}_k,k)-\frac{T}{NV}\sum _K\frac{i{\omega }_k+{\tilde{\xi }}_k+n_0{g}_K}{E_K^2+{\omega }_k^2}.$$

$${\Sigma }^{(1)}(K)=\frac{T}{V}\sum _{K^{\prime }}\frac{g_{K^{\prime }}({\tilde{\xi }}_{k^{\prime }}^2+{\omega }_{k^{\prime }}^2)}{E_{K^{\prime }}^2+{\omega }_{k^{\prime }}^2}\frac{1}{{\tilde{\xi }}_{|{\mathbf{k}}^{\prime }+\mathbf{k}|}-i{\omega }_{k^{\prime }+k}}.$$

$${\Omega }_{\mathbf{v}}=\Omega (\mu \to \mu -m{v}^2/2,{\omega }_k\to {\omega }_k+i\hslash \mathbf{vk}).$$

$${\Omega }_{\mathbf{v}\to 0}=\Omega +NV{n}_{s}m{v}^2/2+o(v^2).$$

$$n_s=n-\left(1-\frac{1}{N}\right)\frac{1}{V}\sum _{\mathbf{k}}n({\tilde{\xi }}_k/T)-\frac{1}{NV}\sum _{\mathbf{k}}\frac{\partial n({\tilde{\xi }}_k/T)}{\partial {\tilde{\xi }}_k}$$

$$\times {\Sigma }_R^{(1)}({\tilde{\xi }}_k,k)-\frac{2T}{3NV}\sum _K{\epsilon }_k\frac{E_K^2-{\omega }_k^2}{(E_K^2+{\omega }_k^2{)}^2}.$$

$$\frac{1}{V}\sum _{\mathbf{k}}n({\tilde{\xi }}_k/T)=\frac{1}{V}\sum _{\mathbf{k}}n({\epsilon }_k/T)-\frac{1}{V}\sum _{\mathbf{k}}\frac{\partial n({\epsilon }_k/T)}{\partial {\epsilon }_k}\tilde{\mu },$$

$$\frac{n_0}{n}\propto \delta t-\frac{4}{{\pi }^{2}N}\delta t\ln \delta t+\dots \simeq (\delta t{)}^{2\beta },$$

$$\frac{n_s}{n}\propto \delta t-\frac{16}{3{\pi }^{2}N}\delta t\ln \delta t+\dots \simeq (\delta t{)}^{2{\beta }_s},$$

$$\beta =\frac{1}{2}-\frac{\sin \left(\frac{D-2}{2}\pi \right)\Gamma (D-2)}{2\pi N{\Gamma }^2\left(\frac{D}{2}\right)},$$

$${\beta }_s=\frac{1}{2}-\frac{2\sin \left(\frac{D-2}{2}\pi \right)\Gamma (D-2)}{\pi DN{\Gamma }^2\left(\frac{D}{2}\right)},$$

Full text

Superfluidity, Bose-Einstein condensation, Large-$N$ expansion

I Introduction

Superfluidity is the first observable phenomenon that comes to mind when one talks about macroscopic low-temperature behavior of interacting many-boson systems. Being induced by quantum effects, discovered more than 80 years ago Kapitza ; Allen in liquid 4He and described Tisza ; Landau by means of the two-fluid concept, the superfluid behavior is a consequence of emergence of the off-diagonal (quasi)-long-range order in the Bose system Pitaevskii . Particularly, for spatial dimensionalities larger than two it always goes together with the truly microscopic characteristic of bosonic systems, namely, the Bose-Einstein condensation (BEC) phenomenon. Therefore, it is generally believed that the occurrence of superfluidity should be associated with BEC transition in these systems. This notion is supported by experimental measurements Onofrio ; Raman , results of the Monte Carlo (MC) simulations Kashurnikov ; Arnold ; Prokofev ; Nho and analytical approaches Floerchinger ; Isaule , nonetheless, an interesting alternative points of view often appear Cooper_etal ; Kruglov .

Experimental realization of almost uniform BECs Gaunt and recent measurements of the ground-state Bogoliubov condensate depletion Lopes_etal in the homogeneous three-dimensional Bose gas together with the observation Ville of sound propagation in the two-dimensional superfluids at finite temperatures allow for further experimental investigations of the interplay between quantum and thermal effects in the many-boson systems with continuous translational symmetry. From the theoretical point of view, the account of quantum and thermal fluctuations is fundamentally important for the description of these systems in the whole temperature range including both very low temperatures near absolute zero and a narrow region of critical point. The standard Bogoliubov prescription adopted for thermodynamics of a dilute Bose gas is capable to capture only the low-temperature properties of a system and in order to obtain the correct (at least on the qualitative level) behavior near the BEC transition temperature one has to use more sophisticated treatments. An extension of Bogoliubov’s theory that naturally combines quantum and thermal fluctuations in the simplest fashion is given by the large-$N$ expansion approach Andersen ; Chien ; Diehl , which was recently shown to predict a qualitatively correct shift of the critical temperature Hryhorchak and thermodynamics Hryhorchak_19 for the Bose gas with a point-like two-body repulsive potential in a whole range of the interaction parameter. The objective of present study is to explore, by means of the $1/N$-expansion techniques, the superfluid properties of this system.

II Formulation

The discussed model is described by the imaginary-time action Popov , which after application of the Hubbard-Stratonovich transformation takes the following form

[TABLE]

where the point-like two-body interaction potentials (both intra- and interspecies) between Bose particles is characterized by a single coupling constant $g/N$. The summation over repeating index $\sigma=1,\ldots,N$ is assumed. Other notations are typical: $x\equiv(\tau,{\bf r})$, integration $\int dx=\int_{0}^{1/T}d\tau\int_{V}d{\bf r}$ is carried out over the $3+1$ ‘volume’ $V/T$ (here $T$ is the temperature of the system), the $1/T$-periodic in imaginary-time variable $\tau$ complex fields $\psi_{\sigma}(x)$ describe $N$ species of bosons and the second-order differential operator $\xi=-\hbar^{2}\nabla^{2}/2m-\mu$ contains the chemical potential $\mu$ that fixes the total number of particles in a system.

In the low-temperature phase, the Bose condensate occurs. It means that the spatially uniform part $\psi_{0}$ of fields $\psi_{\sigma}(x)=\psi_{0}+\tilde{\psi}_{\sigma}(x)$ describing particle degree of freedom should be singled out. We also have to explicitly separate the saddle-point value $\varphi_{0}$ of the real field $\varphi(x)=\varphi_{0}+\tilde{\varphi}(x)$ that incorporate the collective behavior of the system. After that, the shifted complex fields $\tilde{\psi}_{\sigma}(x)$ can be integrated out to yield the following effective action

[TABLE]

where $\tilde{\xi}=\xi|_{\mu\to\tilde{\mu}}$, ($\tilde{\mu}=\mu-i\varphi_{0}$), $n_{0}=|\psi_{0}|^{2}$ is the Bose condensate density of each constituent and ${\rm Sp}\ldots$ denotes the trace of the appropriate differential operator. Equations $-\left(\partial\Omega/\partial\mu\right)_{n_{0},\varphi_{0}}=NVn$, $\left(\partial\Omega/\partial\varphi_{0}\right)_{\mu,n_{0}}=0$, $\left(\partial\Omega/\partial n_{0}\right)_{\mu,\varphi_{0}}=0$ relate the Bose condensate density $n_{0}$, the chemical potential of the system and total particle number $NVn$ to non-zero expectation value of field $\varphi(x)$.

The further analysis trivializes in the large-$N$ limit. The main contribution, which is of order $N$, is now coming from the ideal-gas term, while the quadratic in fields $\tilde{\varphi}(x)$ part of action $S_{\rm eff}$, which impacts the terms of order unity in the thermodynamic potential, can be also exactly taken into account

[TABLE]

where the ‘four-vector’ $K=(\omega_{k},{\bf k})$ stands for the bosonic Matsubara frequency and wave-vector. Here, $\tilde{\xi}_{k}=\varepsilon_{k}-\tilde{\mu}$, $\varepsilon_{k}=\hbar^{2}k^{2}/2m$ is the free-particle dispersion and note that in both sums the terms with ${\bf k}=0$ are omitted. Additionally, we have introduced the polarization operator $\Pi(K)=\Pi_{0}(K)+\Pi_{T}(K)$, where

[TABLE]

denotes the condensate contribution to the dynamic structure factor of ideal Bose gas and

[TABLE]

represents the impact of density fluctuations of the thermally stimulated single-particle excitations. It should be understood that the structure of the higher-order terms in (II) is clear Vakarchuk_12 and there are no principal problems to include them. In practice, however, the numerical computations of the beyond-$1/N$-terms to the thermodynamic characteristics are substantially complicated.

Minimization of $\Omega$ with respect to parameters $\varphi_{0}$ and $\psi_{0}$ leave us with $\varphi_{0}=-ing$ (and appropriately $\tilde{\mu}=\mu-ng$) and

[TABLE]

where we have introduced the effective two-body potential $g_{K}\equiv g/[1+g\Pi_{T}(K)]$ induced by the density fluctuations of non-condensed particles and the frequency-dependent quantity $E^{2}_{K}=\tilde{\xi}^{2}_{k}+2n_{0}g_{K}\tilde{\xi}_{k}$ which coincides with the Bogoliubov spectrum at absolute zero. The Matsubara frequency sum in the above formula should be calculated with factor $e^{i\omega_{k}\tau}$ ($\tau\to+0$) and even after that the integral over the wave-vector is divergent. This divergence is caused by the replacement of the original two-body short-range potential by the $\delta$-function. Therefore, we additionally have to rewrite everywhere in Eq. (2.6) the ‘bare’ coupling constant $g$ via the $s$-wave scattering length $a$ in the adopted approximation

[TABLE]

The subsequent substitution in Eq. (2.6) guaranties finiteness of the chemical potential Hryhorchak_19 . Finally, by using the well-know identity which relates the derivative of the thermodynamic potential with respect to the chemical potential to the particle number, we obtain the condensate density in the adopted approximation

[TABLE]

Here, $n(x)=1/(e^{x}-1)$ is the Bose distribution; the last sum should be carried out with factor $e^{+i\omega_{k}0}$ and $\Sigma^{(1)}_{R}(\omega,k)$ denotes the real part (after analytical continuation in the upper complex half-plane $i\omega_{k}\to\omega+i0$) of the leading-order self-energy $\Sigma(K)=\Sigma^{(1)}(K)/N+\ldots$ of normal Green’s function Hryhorchak_19

[TABLE]

It should be noted that Eq. (II) for the condensate density can be directly obtained by means of Green’s function technique Hryhorchak_19 , which proves the consistency of the large-$N$ expansion approach. Now the interaction-renormalized temperature of BEC transition $T_{c}$ can be found from the condition $n_{0}=0$ and is determined by the following relation: r.h.s. of Eq. (II)$=0$.

The standard procedure to obtain the density of superfluid component is well-described in literature and lies in the calculation of response of the thermodynamic potential to a slow motion of the whole system with velocity ${\bf v}$. In general, this shift of $\Omega$ can be accounted by the gauge transformation of the initial bosonic fields $\psi_{\sigma}(x)\to e^{im{\bf v}{\bf r}/\hbar}\psi_{\sigma}(x)$, $\psi^{*}_{\sigma}(x)\to e^{-im{\bf v}{\bf r}/\hbar}\psi^{*}_{\sigma}(x)$ in action (II) with the subsequent repetition of the calculation scheme (II), (II). From the practical point of view the problem reduces to the replacement of original chemical potential and Matsubara frequencies in the thermodynamic potential

[TABLE]

The effectiveness of the above-described prescription for the superfluid density calculation and for the evaluation of the exact identities relating macroscopic observables to parameters of the low-lying excitations was previously proven for one Pastukhov_q2D ; Pastukhov_InfraredStr and two-component Pastukhov_twocomp ; Konietin Bose systems. In general, the superfluidity in the $N$-component Bose mixture, due to the Andreev-Bashkin effect, is characterized by the $N\times N$ symmetric matrix $n^{\sigma\sigma^{\prime}}_{s}$ of superfluid densities. In our case, however, by supposing the same velocity for all $N$ species we smear out the information about the occurrence of this drag effect providing the diagonal structure of matrix $n^{\sigma\sigma^{\prime}}_{s}=\delta_{\sigma\sigma^{\prime}}n_{s}$ with $N$ times degenerated eigenvalue.

The further consideration should not cause any principal complications and by the brute force expansion of $\Omega_{\bf v}$ at small ${\bf v}$ we obtain

[TABLE]

where $\Omega$ is given by Eq. (II) and the calculations of superfluid density $n_{s}$ per each sort of bosons yield

[TABLE]

The latter formula together with the condensate density (II) represent the main result of this article that will be used below for the numerical computations. But before we proceed to numerics few remarks should be made. The first important thing one should keep in mind is that the condensate and superfluid densities given by formulae (II), (II) both have additional irrelevant terms of higher orders in powers of parameter $1/N$. Indeed, these two equations contain the ideal-gas dispersion $\tilde{\xi}_{k}$ with the $1/N$-corrected effective chemical potential $\tilde{\mu}$ [see Eq. (2.6)]. Therefore, terms of order unity in Eqs. (II), (II)

[TABLE]

should be expanded in the Taylor series up to the first order in $\tilde{\mu}$. By direct comparison of Eq. (2.6) and Eq. (2.9) it is easy to verify that $\tilde{\mu}=\Sigma^{(1)}(0)/N$, which not only ensures the fulfilment of the Hugengoltz-Pines theorem but also provides the infrared convergence of integrals for $n_{0}$ and $n_{s}$. Particularly, it means that at temperatures not too close to $T_{c}$ the difference $\Sigma^{(1)}_{R}(\varepsilon_{k},k)/N-\tilde{\mu}$ (which is under the wave-vector integral in the second terms of (II) and (II) (after substitution of $\tilde{\mu}$) is proportional to $k^{2}$, while $k$ goes to zero. In all other terms in Eqs. (II), (II), which are already of order $1/N$ one is free to replace $\tilde{\xi}_{k}$ by $\varepsilon_{k}$.

Secondly, it is clearly seen from the above expressions that the superfluidity in the leading order of the large-$N$ expansion vanishes with vanishing of a condensate, i.e., the temperature of superfluid transition coincides with the critical temperature of the BEC. This commonly accepted idea about the nature of superfluid transition in three-dimensional Bose systems is in contradiction with some effective field theory studies Cooper_etal . To see the coincidence of the BEC and superfluid transition temperatures in our approach, it is enough to realise that $E_{K}\to\tilde{\xi}_{k}\to\varepsilon_{k}$ (when $n_{0}\to 0$) and that last terms in (II) and (II) are equal to each other in the vicinity of critical point.

The third feature of the whole previous analysis is that the thermodynamic characteristics of $N$-component bosons with the fully symmetric interaction are expressed (in the large-$N$ limit) as a systematic series expansion in integer powers of $1/N$.

III Numerical results and discussion

Taking into account the aforementioned discussion we are now in position to present the graphical results of the numerical integration and Matsubara frequency summation. To improve the convergence of these computations we have used the peculiarities of the $1/N$-expansion. In particular, it is worth mention that in the leading order, this theory applied to Bose systems incorporates two main ingredients, namely, the Bogoluibov theory at very low temperatures and inclusion of the finite-temperature density fluctuation of the non-condensed particles in the simplest approximation. Close to the Bose-Einstein transition temperature these fluctuations are highly developed providing the non-trivial critical behavior of the system. In what follows for the numerical purposes, one has to single out the finite-temperature Bogoliubov condensate depletion in (II) and the Landau expression (with Bogoliubov’s spectrum) for the superfluid density in (II) and analytically perform the summations over frequencies. Now, the summands in remaining terms are well-behaviored at ultraviolet region and can be easily computed by the numerical methods.

Additionally, it is important to consider the zero-frequency terms explicitly because they are actually responsible for the non-analytic behavior of all thermodynamic functions at the critical point and therefore have to be discussed in more detail. From equations determining $n_{0}$ and $n_{s}$ we find out that first three terms in r.h.s. of Eqs. (II), (II) vanish linearly in the vicinity of $T_{c}$, while the last ones demonstrate logarithmic nonanalicity. It is easy to obtain the coefficient standing in front of this log-linear term and hinting the universal power-law behavior

[TABLE]

[TABLE]

where $\delta t=\frac{T_{c}-T}{T_{c}}$ is the dimensionless deviation from the critical temperature. Moreover, these estimations can be easily extended to a more general case of arbitrary dimension $D>2$, where the Bose condensation phenomenon occurs. All we need is to evaluate the small-$k$ behavior of the ‘thermal’ part of the polarization operator at zero Matsubara frequency $\Pi_{T}(K)|_{\omega_{k}=0}\propto k^{D-4}$. Then the substitution in the last term of equations for $n_{0}$ and $n_{s}$ (here prefactor $2/3$ should be replaced by $2/D$) leaves us with the result for critical exponents calculated in the accepted approximation

[TABLE]

[TABLE]

where $\Gamma(x)$ is the gamma function. These findings particularly confirm the exact scaling relation for the superfluid density, $n_{s}/n\propto(\delta t)^{2\beta-\nu\eta}$ Lifshitz , where $\nu$ and $\eta$ are correlation length and Fisher exponents, respectively. Indeed, taking into account values of $\nu$ and $\eta$ calculated to leading order in $1/N$ Abe , we obtain our estimation for $\beta_{s}$.

Having found out the slope of the condensate and superfluid density curves close to $T_{c}$ we may proceed by presenting the results of full numerical computations. Putting $N=1$ we have calculated, for three values of the gas parameter $a^{3}n=10^{-6}$, $10^{-4}$ and $10^{-2}$, the temperature dependence of both $n_{s}$ and $n_{0}$ which are depicted in Figs. 1-3. We also plotted (symbols) results of quantum Monte Carlo simulations taken from Pilati_thesis .Note that temperature in MC data points as well as in our curves is measured in units of $T_{c}$, where $T_{c}$ in our case is the $1/N$-shifted Hryhorchak critical temperature of the system, while positions of MC points are rescaled to the BEC transition temperature obtained Pilati in simulations. Of course, there is some discrepancy (maximum $\sim 5\%$) between these transition temperatures, which however, cannot substantially affect the character of the calculated temperature dependencies.

The presented results of numerical calculations clearly demonstrate the natural tendency: in general, an increase of the interaction strength makes the coincidence with the MC data worse. Particularly, in Fig. 1 we observe a quite good matching in whole temperature region. Regardless the fact that the interaction between particles is very weak the difference of calculated $n_{0}$ with the ideal Bose gas condensate density (dotted line), due to exceptional impact of the density fluctuation of the non-condensed particles, is visible. The curves in Fig. 1 should be compared to results Capogrosso-Sansone obtained by means of the Beliaev technique extended on the finite-temperature region and combined with classical MC computations Prokofev . In principle, for such weakly-interacting systems all approximate approaches Kita ; Yukalovs_1 ; Watabe_13 ; Watabe_14 ; Yukalovs_2 work well and the only exception is the narrow neighborhood of the critical point, where results can vary. Even the standard Bogoliubov theory, which predicts the same critical temperature as in the ideal Bose gas, adequately describes the condensate depletion behavior of interacting particles at non-zero $T$ not too close to $T_{c}$ (at $a^{3}n=10^{-6}$ the Bogoliubov curve lies slightly $\sim 0.01$ above the calculated one).

The quantitative agreement with MC results, as it is demonstrated in Figs. 2 and 3, is found to be worse with the growth of strength of the two-body repulsive potential. It is not surprising because at $a^{3}n=10^{-2}$ even the quantum depletion of the Bose condensate predicted by the large-$N$ expansion approach (which at $T=0$ is equivalent to the Bogoliubov theory) substantially differs from findings of the essentially exact MC simulations. The calculated temperature dependencies of the Bose condensate, however, are qualitatively correct. In order to show this we built in the inset of Fig. 3 the condensate density temperature behavior, rescaled to the zero-temperature value of $n_{0}$ obtained in MC Giorgini simulations. A comparison of regimes of various coupling strengths for $n_{s}$ and $n_{0}$ leads us to conclusion that except for the ground-state depletion the temperature behavior of the Bose condensate is less sensitive to the interaction effects. The latter suggests that even for intermediate couplings the $1/N$ expansion in its simplest approximation is capable to capture the finite-temperature behavior of the Bose condensate density, while mechanisms responsible for the formation of the observed temperature dependence of the superfluid component in Bose systems with nonweak interparticle repulsion require further investigations.

IV Conclusions

In conclusion, we have shown that the combination of two simple ingredients, namely, the Bogoliubov theory at absolute zero and the inclusion of density fluctuations of the noncondensed particles at finite temperatures, that is incorporated in the simplest non-trivial approximation of the $1/N$-expansion method, is a promising tool for the quantitative study of the Bose-Einstein condensation phenomenon and emergence of superfluidity in Bose systems with short-range repulsion. Particularly, using the large-$N$ expansion to leading order we have calculated the temperature dependence of the superfluid and condensate densities for Bose system with point-like repulsive interaction between particles and demonstrated the efficiency of the presented approach by comparing our results to the Monte Carlo simulations data.

Acknowledgements

We are grateful to Prof. Sebastiano Pilati for providing us with the results of his Monte Carlo simulations. The author also thank Prof. Andrij Rovenchak for permanent help during the completion of this research. This work was partly supported by Project FF-30F (No. 0116U001539) from the Ministry of Education and Science of Ukraine.