Macroscopic properties of triplon Bose-Einstein condensates

As it has been shown by Yukalov [28] the spontaneous gauge symmetry breaking is the necessary and sufficient condition for Bose-Einstein condensation and can be realized by the Bogoliubov shift in the field operator as

where for the uniform case the function $v({\mathbf{r}})$ is a constant defining the density of condensed particles as

Since by the definition of the average of $\psi^{\dagger}({\mathbf{r}})\psi({\mathbf{r}})$ is the total number of particles:

with the density of triplons per unit cell $\rho=N/V$, from normalization condition

one immediately obtains

Therefore the field operator $\widetilde{\psi}(\mathbf{r})$ determines the density of uncondensed particles. Note that $\widetilde{\psi}(\mathbf{r})$ and the condensate function should be considered as independent variables being orthogonal to each other in terms of the expectation values:

Consequently, passing to the momentum space

and inserting (10) into (8) one can easily see that in (10) the summation by momentum in the finite volume should not include $\mathbf{k=0}$ states. Below we indicate this rule by introducing prime sign in the momentum summation.

Now using (4) and (10) in (1) we present the Hamiltonian as the sum of five terms

labeled according to their order with respect to $a_{\mathbf{k}}$ and $a^{\dagger}_{\mathbf{k}}$. The zero-order term does not contain the field operators of uncondensed triplons

The first order term is identically zero

due to the orthogonality condition (8). The second-order term is

The third-order

as well as fourth-order

terms are rather complicated and a diagonalization procedure is needed. For this purpose in HFB approximation the following procedure is usually implemented [33]:

where $\langle a_{\mathbf{k}}^{\dagger}a_{\mathbf{p}}\rangle=\delta_{\mathbf{k},\mathbf{p}}n_{\mathbf{k}},\quad\langle a_{\mathbf{k}}a_{\mathbf{p}}\rangle=\delta_{\mathbf{k},-\mathbf{p}}\sigma_{\mathbf{k}}$ with $n_{\mathbf{k}}$ and $\sigma_{\mathbf{k}}$ being related to the normal ($\rho_{1}$) and anomalous ($\sigma$) densities as

Here we underline that the main difference between the HFP and HFB approximations concerns the anomalous density: neglecting $\sigma$ as well as $\langle a_{\mathbf{k}}a_{\mathbf{p}}\rangle$ in (14), (17), and (18) one arrives at the HFP approximation, which can also be obtained in variational perturbation theory [34]. However, the normal, $\rho_{1}$, and anomalous averages, $\sigma$, are equally important and neither of them can be neglected without making the theory not self-consistent [35, 36, 37]. Altough $\rho_{1}$ and $\sigma$ are functions of temperature and external magnetic field, we omit explicit dependences in the formulas when it does not cause a confusion.

Taking into account that $E_{3}=\langle H_{3}\rangle=0$, the expectation value of the energy can be evaluated as $E(T\leq T_{c})=E_{0}+E_{2}$, where

and $E_{2}=\langle\widetilde{H}_{2}\rangle$, with $\widetilde{H}_{2}$ being quadratic in $a_{\mathbf{k}},a_{\mathbf{k}}^{\dagger}$:

In the last equation

The next step, in both approximations, is the Bogoliubov transformation

to diagonalize $\widetilde{H}_{2}$. The operators $b_{\mathbf{k}}$ and $b_{\mathbf{k}}^{\dagger}$ can be interpreted as annihilation and creation operators of phonons with following properties:

where $\beta\equiv 1/T$. To determine the phonon dispersion $E_{k}$ we insert (24) into (21) and require that the coefficient of the term $b_{\mathbf{k}}b_{-\mathbf{k}}+b_{-\mathbf{k}}^{\dagger}b_{\mathbf{k}}^{\dagger}$ vanishes, i.e:

Now using the condition $u_{\mathbf{k}}^{2}-v_{\mathbf{k}}^{2}=1$ and presenting $u_{\mathbf{k}},v_{\mathbf{k}}$ as

yields

that is

where $\omega_{k}$ and $\Delta$ are given in Eqs. (22) and (23).

Further requirement for $E_{k}$ follows from the Hugenholtz-Pines theorem [38]: at small momentum $k$ the spectrum should be gapless, and, therefore, the phonon dispersion is linear: $E_{k}\sim ck+O(k^{2})$, where $c$ can be considered as the sound speed²²2It can be shown that [39, 35] $\Delta$ is related to the normal ($\Sigma_{\rm n}$) and anomalous ($\Sigma_{\rm a}$) self-energies as $\Sigma_{\rm n}=\Delta+\mu$ and $\Sigma_{\rm a}=\Delta$, respectively. The last two equations give $\Sigma_{\rm n}-\Sigma_{\rm a}=\mu$ which is again in agreement with Hugenholtz-Pines theorem [40, 41]. This linearity can be achieved by setting

which together with Eq.(23) yields

With this choice one obtains

with the sound speed

i.e. $\Delta=mc^{2}$. Here $m$ has the meaning of the triplon effective mass, characterizing the dispersion in the limit of small momenta $\varepsilon_{k}\approx k^{2}/2m$.

Having fixed $u_{\mathbf{k}}$ and $v_{\mathbf{k}}$ one can find normal and anomalous densities as well as the energy by using equations (19)-(28). For example, inserting (24) into (19) results in

where we used the relation $u_{\mathbf{k}}v_{\mathbf{k}}=-{\Delta}/{2E_{k}}$ and introduced notation $W_{k}={1}/{2}+f_{B}(E_{k})$. Similarly one obtains:

with  $\omega_{k}=\varepsilon_{k}+\Delta$.

When the bare dispersion $\varepsilon_{k}$ is isotropic, the momentum summation can be done with:

As a result, at $T=0$ the quantities in Eqs.(35)-(37) can be represented as:

The divergences in these integrals can be regularized by introducing a cutting parameter $\Lambda$ ($\Lambda\rightarrow\infty$ at the end of the calculations) or equivalently by using the dimensional regularization scheme. Now, assuming for the moment that, for $T=0$ case $\varepsilon_{k}={k}^{2}/{2m}$ and using dimensional regularization one obtains: ³³3Note that the second integral in Eq. (41) can be treated with the Veltman formula, see e.g., H. Kleinert Hagen and V. Schulte-Frohlinde, Critical properties of $\phi^{4}$-theories, World Scientific Publishing Company (2001).

Summarizing this subsection we rewrite the above formulas as:

where Bose distribution of phonons $f_{B}\left(E_{{k}}\right)$ is defined in Eq.(26).

To perform the MFA calculations one starts by solving Eqs. (23) and (32) with $\rho_{1}$ and $\sigma$ given by Eqs. (45) and (46). As outlined above, in contrast to the BEC of atomic gases, in the triplon problem the chemical potential $\mu$ is the input parameter, whereas the densities are the output ones. Bearing this in mind, we rewrite the main Eqs. (23) and (32) as

The system of coupled equations, (45), (46) and (48), (49) has to be solved for given  $T$ and $\mu$ to evaluate the triplon density

which is proportional to the measured sample magnetization density $M$. Note that by formally setting in all above formulas $\sigma\equiv 0$, one arrives at the HFP approximation and particularly

Before discussing the normal $T>T_{c}$ phase we evaluate the temperature of BEC transition. It is well-known that in the MFA the system of interacting Bose condensate particles at $T\rightarrow T_{c}$ behaves like an ideal gas. In fact, assuming $\rho_{0}(T\rightarrow T_{c})=0$, $\rho_{1}(T\rightarrow T_{c})=\rho_{c}$, $\sigma(T\rightarrow T_{c})=0$, $\Delta(T\rightarrow T_{c})=0$, and $E_{k}=\varepsilon_{k}$ one concludes from Eqs.(45),(48) that

With given $\mu=\mu_{B}gH_{\rm ext}-\Delta_{\rm st}$ and coupling constant $U$, the critical temperature $T_{c}\equiv 1/\beta_{c}$ can be found as a solution of Eq. (52). Note that for the parabolic dispersion $\varepsilon_{k}=k^{2}/{2m}$, the integral (52) can be evaluated analytically giving following well known relation [42]:

where $\zeta(x)$ is the Riemann function.

As it has been underlined in Introduction, the bare dispersion of magnons, $\varepsilon_{k}$, plays the crucial role in determining $T_{c}$. Figure 2 presents the transition temperature as a function of magnetic field $H_{\rm ext}$ and corresponding triplon density $\rho$, that is the sample magnetization, for parabolic $\varepsilon_{k}={k}^{2}/{2m}$ and ”relativistic”

bare dispersion, typical for a system with a gapped spectrum [24, 43]. Here the effective exchange parameter $J=2\sqrt{\Delta_{\rm st}/m}$ is chosen to match the parabolic and the relativistic $\varepsilon_{k}$ at small $k$. The dashed line obtained directly from Eq.(53) displays $T_{c}\sim\rho^{2/3}$ behavior. The solid one is a result of numerical solution of Eq.(52) with the dispersion in Eq.(54) and shows a crossover from $T_{c}\sim\rho^{2/3}$ at lower to $T_{c}\sim\rho^{0.5}$ at higher temperatures in agreement with the experimental data.⁴⁴4We mention that for linear dispersion $\varepsilon_{k}\sim k$, the exponent would be equal to $1/3$, and, therefore, in the experimental regime, the triplon gas is between the nonrelativistic and strongly relativistic realizations. The decrease in $\Delta_{\rm st}$ enhances the role of relativistic features in the spectrum and leads to a faster deviation from the $T_{c}\sim\rho^{2/3}$ behavior, as can be seen in Fig.2. Here and below we mainly use the set of input parameters as $m=0.0204$ K^-1, $\Delta_{\rm st}=7.1$ K, $U=313$ K, and $g=2.06$ [32] valid for the weakly anisotropic quantum antiferromagnet TlCuCl₃.

In the normal phase the symmetry in Eq.(3) is not broken and the Bogoliubov shift is not needed either. Here the anomalous density vanishes, $\sigma(T>T_{c})=0$, and hence, both approximations, HFB and HFP coincide. As a result we obtain for the triplon density

with $\omega_{k}=\varepsilon_{k}-\mu+2U\rho\equiv\varepsilon_{k}-\mu_{\rm eff}$. Similarly, by setting $\Delta=\rho_{0}=\sigma=0,\ \rho_{1}=\rho,\ E_{{k}}=\omega_{{k}}$ in Eq.(47) one obtains following equation for the energy per unit cell

The density of triplons as a function of temperature, evaluated in the HFP and HFB approximations is presented in Fig.3. It is seen that the former leads to a discontinuity in the magnetization near the critical temperature, while the latter gives a continuous behavior in accordance with the experimental data [32, 44]. However, in the condensate phase, the triplon density is higher in the HFB than in the HFP approximation. Therefore, the validity of the description of the BEC in the HFP approximation can be checked in the magnetization measurement experiments.

Triplon contribution to the constant volume specific heat $C_{v}$, can be calculated by differentiation of the energy with respect to the temperature at given chemical potential, that is at given external field [45, 46]:

where the energies for condensate and normal states are given by Eqs. (47) and (56), respectively.

Before discussing the magnetic susceptibility, we note that the macroscopic properties of the systems are related to their response to external fields. An example is given by the isothermal compressibility

If $\kappa_{T}\rightarrow\infty$ the system becomes unstable [46, 47] since an infinitesimal fluctuation of pressure $P$ will lead to its collapse or explosion. As to the triplons their density is proportional to the magnetization $M$, while the $H_{\rm ext}$ simulates the pressure, and hence the relevant parameter, which determines the stability of the system with respect to the magnetic field, is the susceptibility

The susceptibility can be calculated with the triplon density as:

We will omit $H_{\rm ext}$ from the arguments of $\chi(H_{\rm ext},T)$ and $C_{v}(H_{\rm ext},T)$ below.

We begin with the normal phase, where $\rho_{0}=\sigma=\Delta=0$. The derivative of the density with respect to the temperature can be evaluated here with Eq.(LABEL:2.50) as

Calculating the derivative of Eq.(56) with respect to $T$ one obtains:

Here we introduced dimensionless quantities

and used the relation ${\partial f_{B}(\omega)}/{\partial T}=-\beta\omega f_{B}^{\prime}(\omega)$, where $f_{B}^{\prime}\left(\omega\right)={\partial f_{B}(\omega)}/{\partial\omega}=-\beta\exp({\beta\omega})f_{B}^{2}\left(\omega\right)$ and $\omega_{k}$ is given in Eq.(22).

Similarly, taking the derivative of the self-consistency Eq.(LABEL:2.50), which we present here in the form:

with respect to $\mu$ and solving the equation for $\partial\rho/\partial\mu$ one finds $\partial\rho/\partial\mu=S_{2}/U(2S_{2}-1)$. As a result, we obtain with Eq.(60) the normal phase susceptibility:

Now we proceed with the condensate phase, where the dependence of $\Delta$ and the corresponding normal and anomalous densities on temperature and magnetic field should be taken into account.

Here $C_{v}(T)$ is obtained by taking the derivative of $E\left(T<T_{c}\right)$ given by Eq. (47). The latter can be rewritten as:

where we used the relation $\rho_{0}=\mu/U-2\rho_{1}+\sigma$ and present ${E}_{2}$ defined in Eq.(37) as:

As a result:

We begin with calculation of ${\partial\Delta}/{\partial T}$, which is the key ingredient in the specific heat. It is obtained by differentiating both sides of the main equation (48) with respect to the temperature, and for known ${\partial\Delta}/{\partial T}$ the other two derivatives ${\partial\rho_{1}}/{\partial T}$ and ${\partial\sigma}/{\partial T}$ can be evaluated directly from Eqs.(45), (46). As a result one obtains the derivatives necessary to evaluate the specific heat:

where we introduced notation for the frequently used expression:

The expression for

completes the set of derivatives necessary for evaluation of the specific heat.

The derivative ${\partial\Delta}/{\partial\mu}$ can be found by differentiating both sides of the main equation (48) with respect to $\mu$. As a result,

As to ${\partial\rho}/{\partial\mu}$, needed for evaluation of the magnetic susceptibility in Eq.(60), it is obtained by differentiating the equation $\rho=\left(\Delta+\mu\right)/2U,$ with respect to $\mu$. This yields

Below we apply these equations to the cusps in the specific heat and susceptibility and to the qualitative effects such as the instability in strong magnetic fields. We will show that at a given temperature, there exists a critical field $H^{\rm cr}_{\rm ext}(T)$, such that when $H_{\rm ext}$ approaches $\geq H^{\rm cr}_{\rm ext}(T)$, the magnetic susceptibility, $\chi(H_{\rm ext},T)$ diverges.

The cusp in the specific heat defined as

is an interesting quantity in the theory of phase transitions. In accordance with the Ehrenfest classification, a phase transition with the discontinuity in $C_{v}$ near the transition point, is the second order one. Particularly, it is well-known that [45, 46] for the ideal gas $\Delta C_{v}$ near the transition into BEC is zero i.e. the specific heat is continuous. We shall illustrate this fact for completeness. The specific heat per particle of the ideal Bose gas is given by [45, 46]

where $g_{p}(z)$ is defined as:

related to the Riemann function as $\zeta(p)=g_{p}(z=1)$. When $T\rightarrow T_{c}+0$, the fugacity $z=\exp({\beta\mu})<1$ tends to unity, i. e. $\displaystyle\lim_{T\rightarrow T_{c}+0}z=1$ and the second term in (78) vanishes, because of the divergence in $g_{1/2}$, i.e. $g_{1/2}(z)\sim(1-z)^{-1/2}$ in this region, while the first term exactly coincides with (LABEL:3.18). As a result,

and hence the specific heat of the ideal gas is continuous although being plotted as a function of temperature $C_{v}^{U=0}(T)$ behaves similarly to the $\lambda$ curve.

We proceed with calculation of $\Delta C_{v}$. Bearing in mind that for $T>T_{c}$ HFB and HFP approximations coincide, from (62) we obtain

For $T<T_{c},$ $C_{v}$ is given by Eqs. (68)-(72). Assuming in the last equations $E_{{k}}=\varepsilon_{k}$ and $\Delta=0,$ one finds

which is finite. Using (82) in (70)-(72) and setting in (68) $\rho=\mu/2U$, $\rho_{0}=0,\rho_{1}=\rho,\sigma=0$, one obtains for the maximum value of the specific heat $\mbox{max}\{C_{v}^{\rm HFB}(T)\}\equiv\lim_{T\rightarrow T_{c}-0}C_{v}(T)$:

Subtracting (81) from (83) we finally obtain

where ${\partial\Delta}/{\partial T}$ is given by Eq. (82).

In a quite similar way one can calculate the cusp in the susceptibility:

For the parabolic dispersion, calculations can be done analytically, as presented below. The specific heat is expressed as: