Dimensional Crossover in the Bose–Einstein Condensation Confined to Anisotropic Three-Dimensional Lattices

Abstract

The Bose–Einstein condensation (BEC) in three-dimensional (3D) anisotropic lattices is studied. We present theoretical results for the critical temperature for BEC, chemical potential, condensate fraction and relevant thermodynamic quantities like: internal energy, entropy, specific heat and compressibility as a function of anisotropy parameter being the ratio of the nearest-neighbor in-plane ($t_\parallel$) and out-of-plane ($t_\perp$) hopping amplitudes. In particular, considered scenarios include weakly coupled two-dimensional (2D) planes ($t_\perp /t_\parallel \ll 1$, relevant for layered structures) as well as a rod-like geometry of interacting one-dimensional (1D) chains ($t_\parallel /t_\perp \ll 1$). The impact of the dimensional crossover as the system is tuned away from a set of disconnected 2D layers, or traverses from a set of separate 1D chains to a regime where a fully isotropic 3D structure emerges is elucidated. Both numerical and analytic approaches are employed, (the latter in a form of series expansions involving $t_\parallel ,t_\perp$ amplitudes) for internal energy, entropy, specific heat and isothermal compressibility. The theoretical outcome of the present study may be of interest to a number of scenarios in solid-state physics, where the relevant quasi-particles are bosonic-like, as well as might be applicable to the physics of cold bosons loaded in artificially engineered 3D optical lattices.

1 Introduction

The low energy properties of many condensed matter systems are naturally described in terms of bosonic particles [1]. This is not only the case of cold bosonic atoms loaded in optical lattices [2], but effective bosonic models have also been provided an adequate description of several systems such as superconducting Josephson junction arrays [3], or quantum magnets in a field [4]. In this context, the study of the phenomenon of the Bose–Einstein condensation (BEC) is an active area of research both in atomic [2, 5, 6] and condensed matter physics [4, 7,8,9,10,11]. The fascinating aspect of Bose–Einstein condensates is that a macroscopic number of particles act collectively forming a coherent matter wave. For finite temperature condensates, only part of the bosonic entities is participating in the macroscopic occupation of the ground state—a phenomenon which is the subject of current theoretical and experimental research activity in the context of quantum phase transitions. However, the nature of the latter strongly depends on space dimension. One can rigorously show using the Mermin–Wagner theorem [12,13,14] that there is no condensation with long-range phase coherence in the one- (1D) or two-dimensional (2D) bosonic systems at finite temperature. Owing to thermal fluctuations, in 1D or 2D systems bosons cannot undergo a conventional phase transition associated with the breaking of a continuous U(1) symmetry. Nevertheless, planar 2D systems may exhibit a kind of unconventional phase transition to a state with quasi-long-range order and this kind of phase transition was explored in 2D condensates [15].

Since many bulk three-dimensional physical systems exhibiting condensation of bosonic modes actually consist of weakly coupled (rather than isolated) planes (or chaines), one can wonder to what extent finite anisotropy may give observable signatures regarding BEC in a layered three-dimensional (3D) arrangement [2, 16]. This issue has raised a lot of interest in connection with the high-temperature superconductivity (HTSC). While the search for the mechanism of HTSC is an ongoing effort which has yet to reach a consensus, many properties of high-temperature superconductors may be understood in terms of Bose–Einstein condensation of electron (hole) pairs. It is evident that the HTSC materials are highly anisotropic, with much of the physics being governed by the coupling of two-dimensional structures [16]. There is significant evidence that in cuprates, the layers are Josephson coupled and there is a significant anisotropy. The ratio of the $\lambda _\mathrm{c}$ penetration depth (measured along the $c-$axis, perpendicular to the cooper-oxide $ab-$plane) and the $ab-$plane $\lambda _{ab}$ penetration depth: $\lambda _\mathrm{c}/\lambda _{ab}$ ranges from $\lambda _\mathrm{c}/\lambda _{ab}=5$ for YBa$_2$Cu$_3$O$_7$ through $\lambda _\mathrm{c}/\lambda _{ab}=25$ for Tl$_2$Ba$_2$Cu$_3$O$_6$ till $\lambda _\mathrm{c}/\lambda _{ab}=250$ for Bi$_2$Sr$_2$CaCu$_2$O$_8$ at the other extreme [17]. Interestingly, it is observed that the superconducting critical temperature strongly depends on the inter-layer structure [18], which means that three-dimensional coupling of planes must play an important role in the onset of superconductivity [19,20,21,22].

Interestingly, the principle for realizing superfluidity have also be extended to layered materials the so-called van der Waals heterostructures, where coupled quantum wells from atomically thin materials are stacked on top of each other [23]. In such semiconductor structures, electrons would accumulate in the opposite monolayers and form bosonic bound states—the called indirect excitons. Because of the indirect exciton mass, the indirect exciton BEC should occur at temperatures typically many orders of magnitude higher than for atoms in optical lattice. Condensation of intra-layer excitons in such heterostructures open a way to realization of superfluidity with far-reaching implications for science and technology [24, 25].

Another promising extension of BEC is to magnons spin-wave quanta that behave as bosonic quasi-particles (magnons) in magnetic nanoparticles [26, 27]. This system has unique characteristics differentiating it from atomic BEC and creating the potential for a variety of interesting behaviors and applications that include high-temperature Bose–Einstein condensation and novel nano-magnetic devices [28].

A particularly well-suited setup for the modeling of spatially periodic systems is provided by neutral ultra-cold quantum gases that are confined in optical lattices. These artificial structures [2], generated by the interference pattern of laser beams, enable the probing of many-body physics in periodic potentials in a controllable manner over a wide range of parameters. With arbitrarily polarized light beams, it is possible to engineer tunneling matrix elements (and thus the dispersion relation) for bosonic quantum gases in optical lattices in a large variety of complex configurations—that are not limited to 2D—as a truly three-dimensional lattice [29]. Besides the obvious relevance to solid-state systems, optical lattices can be used to emulate the crystalline structure of solids with ultra-cold atoms, providing a way to explore the strong anisotropy phenomena typical of condensed matter systems.

This brief account of variety of contexts, where the anisotropy plays a role in the mechanism of Bose–Einstein condensation, shows that the understanding of BEC phenomenon is of great importance in many branches of physics. For example, the finite-size scaling properties have been studied in a superfluid with respect to the film thickness and resulting crossover from 3D to 2D behavior in a film geometry [30, 31]. Dimensional crossover in the Bose–Einstein condensation in quantum gases confined within slab geometries was also studied, where the system experiences a quasi-2D Berezinskii–Kosterlitz–Thouless (BKT) transition [32]. Very recent analysis of a weakly interacting Bose gas shows the crossover from three to low dimensions in the case of a box potential with periodic boundary conditions [33, 34]. For a non-lattice system of trapped cold atoms, the crossover from 3D to 1D or 2D was explored by reducing the number of atoms in condensates which were trapped in highly elongated magnetic traps (1D) and disk-shaped optical traps (2D) [35]. The BEC in a highly anisotropic 3D harmonic trap with tight radial and weak axial confinement was also investigated [36].

Therefore, in a more general context covering the condensed matter physics (and solid-state physics in particular) it will be desirable to explore the role of the lattice anisotropy on the superfluid behavior while the system is tuned from a set of disconnected 2D layers in one scenario or crosses over from a set of 1D chains to a regime where a fully isotropic three-dimensional structure is restored. To this end, in the present work, we study the BEC in a periodic potential of the anisotropic cubic lattice structure, where the degree of the anisotropy can be smoothly varied realizing 1D$\rightarrow$3D$\rightarrow$2D dimensional crossover. Using the grand canonical ensemble, the thermodynamic properties are analyzed in terms of extensive and intensive variables allowing for a detailed account of the problem in terms of equations of state and various thermodynamic quantities calculated via analytic series expansions as well as numerical calculations (with specialized approach to a lattice systems).

The outline of the paper is as follows: In Sect. 2, we introduce the model in which bosons live on the anisotropic three-dimensional lattice. In Sect. 3, we present the basic features of the superfluid properties of Bose gas confined in anisotropic lattice environment, as the crossover of critical temperature, equation of state and the behavior of the chemical potential. In Sect. 4, we present a numerical and analytic study of the model thermodynamics in various lattice settings. The implication of our results is discussed in the conclusions section where possible applications of our findings are also given. Finally, in the “Appendix” section we discuss the technical details regarding the intricate handling of lattice integrals arising in quantum statistical mechanics of bosons in 3D anisotropic structures.

2 Bosons in the Anisotropic Lattice

2.1 Model

We consider bosons moving in an anisotropic cubic lattice. The energy eigenfunctions of a single boson moving in a periodic lattice potential are Bloch waves and energy eigenvalues form bands. The Hamiltonian of non-interacting bosons in an energy band is given by

$$\begin{aligned} \hat{H}=\hat{H}_t+\hat{H}_{\mu }, \end{aligned}$$

(1)

where the terms $H_t$ and $H_\mu$ in second quantized form are

$$\begin{aligned} \hat{H}_t= & {} - \sum _{<i,j>} t_{ij}\hat{ a}^\dagger _i \hat{a}_j, \nonumber \\ \hat{H}_{\mu }= & {} -{\mu }\sum _{i=1}^N \hat{n}_i \equiv -\mu \hat{\mathcal{N}}, \end{aligned}$$

(2)

where $\hat{a}^\dagger _i$ and $\hat{a}_i$ are the bosonic creation and annihilation operators, $\hat{\mathcal{N}}$ is the total particle number operator and $\hat{n}_i =\hat{a}^\dagger _i \hat{a}_i$ is the number operator on the lattice site i, where $i=1,\dots ,N$ with N being the number of sites. Sum over pairs of lattice sites i, j includes only pairs of nearest neighbors and $t_{ij}$ is the hopping amplitude, which is responsible for the tunneling of a boson from one site to another neighboring site. Since it is convenient to work in the grand canonical ensemble, a term with the chemical potential $\mu$ is incorporated into the Hamiltonian. For periodic lattice structure, it is convenient to diagonalize the Hamiltonian in Eq. (2) with the help of Fourier transform,

$$\begin{aligned} t_{ij}=t(|\mathbf{r}_i-\mathbf{r}_j|)=\frac{1}{N}\sum _{\mathbf{k}}t( \mathbf{k}) e^{-i\mathbf{k} \cdot (\mathbf{r}_i-\mathbf{r}_j)}, \end{aligned}$$

(3)

where $\mathbf{r}_i=(x_i,y_i,z_i)$ are lattice vectors and $\mathbf{k}=(k_x,k_y,k_z)$ is the boson quasi-momentum, while $t( \mathbf{k})$ stands for the energy band structure. The kinetic energy dispersion reads

$$\begin{aligned} t( \mathbf{k})=\sum _{i}t( \mathbf{R}_i) e^{-i\mathbf{k} \cdot \mathbf{R}_i} \end{aligned}$$

(4)

where we have defined $\mathbf{R}_i=\mathbf{r}_i-\mathbf{r}_j$. For in-plane $t_\parallel$ and inter-plane $t_\perp$, hopping amplitudes $t( \mathbf{k})$ are given by

$$\begin{aligned} t(\mathbf{k)}=2t_\parallel [\cos (k_x a)+\cos (k_y a)]+ 2t_\perp \cos (k_z a). \end{aligned}$$

(5)

To quantify the degree of the anisotropy defined by the hopping amplitudes, we write Eq. (5) as

$$\begin{aligned} t(\mathbf{k)}=2t_\parallel \left[ \cos (k_x a)+\cos (k_y a)+ \alpha \cos (k_z a)\right] , \end{aligned}$$

(6)

where the dimensionless parameter $\alpha =t_\perp /t_\parallel$ with $0<\alpha <1$ designates the (weakly, for $\alpha \ll 1$) coupled 2D planes scenario, where the natural energy scale is set by the in-plane parameter $t_\parallel$. Alternatively, by rewriting Eq. (5) as

$$\begin{aligned} t(\mathbf{k)}=2t_\perp \left\{ \alpha '\left[ \cos (k_x a)+\cos (k_y a)\right] + \cos (k_z a)\right\} \end{aligned}$$

(7)

with the ratio $\alpha '=t_\parallel /t_\perp$which defines the degree of anisotropy for the 1D chains (oriented along c-direction) coupled (weakly for $\alpha ' \ll 1$) via in-plain interactions.

2.2 The Thermodynamic Grand Potential

The thermodynamics of the system can be drawn from the Hamiltonian in Eq. ( 1) via the partition function $Z=\mathrm{Tr}\exp [-\beta (\hat{ H}-\mu \hat{\mathcal{N}})]$, where $\beta =1/k_\mathrm{B}T$ with T being the temperature. The grand potential $\Omega (\mu , T)=E-TS-\mu \mathcal{N}=-k_\mathrm{B}T\ln Z$, where E stands for the internal energy and S for the entropy, respectively. It is straightforward to show that—in case under study—the grand canonical potential per lattice site $\bar{\Omega }$ is given by

$$\begin{aligned} \bar{\Omega }\equiv & {} \frac{\Omega }{N} = \frac{1}{\beta N} \sum _{\mathbf{k}}^N \ln \left\{ 1-\exp \left[ -\beta \left( \mathcal{E}(\mathbf{k)} -\bar{\mu }\right) \right] \right\} , \end{aligned}$$

(8)

where we have defined the energy $\mathcal{E} (\mathbf{k)}$ that is measured with respect to the bandwidth energy $t(\mathbf{0})$,

$$\begin{aligned} \mathcal{E} (\mathbf{k)} \equiv t(\mathbf{0}) - t(\mathbf{k}). \end{aligned}$$

(9)

so that $\mathcal{E} (\mathbf{0)}=0$. It is evident that for small momenta (that is in the long wave length limit), the dispersion in Eq. (9) behaves as parabolic function with respect to the modulus of the wave vector: $\mathcal{E} (\mathbf{k})\sim \mathbf{k}^2$, which allows to make a link to the continuum “conventional” Bose system. Accordingly, the shifted chemical potential in Eq. (8) is given by

$$\begin{aligned} \bar{\mu } =t(\mathbf{0}) +\mu . \end{aligned}$$

(10)

As will be verified in the next section, the shifted chemical potential is constrained to $\bar{\mu }\le 0$ in order to keep the occupation of the ground state positive. The thermodynamic quantities of interest can be calculated from the thermodynamic potential $\bar{\Omega }$ in Eq. (8).

3 Superfluid Properties

3.1 Critical Temperature

Bose–Einstein condensation corresponds to the build-up of a macroscopic population in the ground state, which occurs below a certain temperature called the critical temperature $T_\mathrm{c}$. According to the rules of statistical mechanics, the total number of particles $\mathcal N$ (or the density of particles $n=\mathcal{N}/N$, or the number of bosons per lattice site) is conjugate to the chemical potential $\mu$. Accordingly, the Bose–Einstein distribution function is given by

$$\begin{aligned} n= & {} -\frac{\partial \bar{\Omega }}{\partial \mu } =-\frac{1}{N}\sum _{\mathbf{k}} \frac{1}{\exp \left[ \beta \left( \mathcal{E}(\mathbf{k)} -\bar{\mu }\right) \right] -1}. \end{aligned}$$

(11)

Fig. 1

Critical temperature $T_\mathrm{c}$ of BEC (normalized to the bandwidth energy $t(\mathbf{0})$) as a function of the anisotropy parameter $\eta$ (see, Eq.  28), quantifying the crossover from the 1D coupled chains $t_\perp >t_\parallel$ to the coupled 2D planes $t_\perp >t_\parallel$ (see, pictorial insets). For the 3D isotropic $t_\perp =t_\parallel$ case $\frac{k_\mathrm{B}T_\mathrm{c}}{t(\mathbf{0})}=0.931831$. The boson occupation number is fixed to $n=1$ (Color figure online)

Fig. 2

Evolution of the critical temperature $T_\mathrm{c}$ as a function of the number of bosons per lattice site n and the anisotropy parameter $\eta$ (see, Eq.  28) (Color figure online)

Below $T_\mathrm{c}$ out of total $\mathcal{N}$ particles, a macroscopic number $\mathcal{N}_0$ particles, where

$$\begin{aligned} \mathcal{N}_0=\frac{1}{\exp \left[ -\beta \left( t(\mathbf{0)}+\mu \right) \right] -1} \end{aligned}$$

(12)

occupy the quantum state with $\mathcal{E}(\mathbf{k=0)}=0$, so that for the total particle number one gets

$$\begin{aligned} \mathcal{N}=\mathcal{N}_0+\sum _{\mathbf{k \ne 0}} \frac{1}{\exp \left[ -\beta \left( t(\mathbf{k)}+\mu \right) \right] -1}. \end{aligned}$$

(13)

We can solve Eq. (12) with respect to the chemical potential $\mu$ by means of the inverse of the exponential function. We obtain

$$\begin{aligned} \mu (\mathcal{N}_0)=-k_\mathrm{B}T\ln \left( \frac{1}{\mathcal{N}_0}+1\right) -t(\mathbf{0)}. \end{aligned}$$

(14)

Therefore, a fixed fraction $n_0$ of the particles in the ground state $n_0\equiv \mathcal{N}_0/N$ in the thermodynamic limit $(N\rightarrow \infty )$ implies $\mathcal{N}_0=n_0N\rightarrow \infty$, so that

$$\begin{aligned} \lim _{\mathcal{N}_0 \rightarrow \infty }\mu (\mathcal{N}_0)=-t(\mathbf{0)} \equiv \mu _0, \end{aligned}$$

(15)

a constant chemical potential within the condensed state. The critical temperature $T_\mathrm{c}$ of BEC can be calculated by considering Eq. (13) in the limit where the number of condensed particles $\mathcal{N}_0$ vanishes, when approaching from the low-temperature ordered state with the chemical potential fixed to $\mu =\mu _0$. The critical temperature can be found from Eq. (11). In the thermodynamic limit $N\rightarrow \infty$, the energy levels are closely spaced, so that the multiple summation over quasi-momenta in Eq.  (11) is replaced with the integral according to

$$\begin{aligned} \lim _{N\rightarrow \infty }\frac{1}{N}\sum _\mathbf{}\ldots =\int \frac{\hbox {d}^3k}{(2 \pi /a)^3 }\ldots \end{aligned}$$

(16)

Therefore, $T_\mathrm{c}$ can be extracted from the equation for the density of bosons n by performing momentum integration.

$$\begin{aligned} n= \int \frac{\hbox {d}^3k}{(2 \pi /a)^3 } \frac{1}{\exp \left[ \beta _\mathrm{c}\left( t(\mathbf{0)}-t(\mathbf{k)}\right) \right] -1}, \end{aligned}$$

(17)

where $\beta _\mathrm{c}=1/k_\mathrm{B}T_\mathrm{c}$.

In order to obtain an analytic estimate of $T_\mathrm{c}$ as a function of hopping amplitudes $t_\perp$ and $t_\parallel$, it is convenient to introduce fugacity $z_0$ defined as

$$\begin{aligned} z_0=e^{\beta _\mathrm{c} t(\mathbf{0)} } \end{aligned}$$

(18)

and rewrite Eq. (17) in the form

$$\begin{aligned} n= \int \frac{\hbox {d}^3k}{(2 \pi /a)^3 } \frac{z_0^{-1}e^{\beta _\mathrm{c} t(\mathbf{k})}}{1-z_0^{-1}e^{\beta _\mathrm{c} t(\mathbf{k})}}. \end{aligned}$$

(19)

Furthermore, by expanding the rational term in Eq. (19) in power series the multiple integrals over quasi-momenta factorize giving

$$\begin{aligned} {n}= \sum _{p=1}^{\infty }z_0^{-p} I_0^2(2p\beta _\mathrm{c} t_\parallel )I_0( 2p\beta _\mathrm{c} t_\perp ), \end{aligned}$$

(20)

where $I_0(x)$ stands for the modified Bessel function of the 0th kind [37]. By observing that the arguments of the Bessel functions in Eq. (20) are dominated by large values, we can utilize the asymptotic form of $I_0(x)$ for $x \gg 1$ given by [37]

$$\begin{aligned} I_0(x) \approx \frac{e^x \left( \frac{1}{x}\right) ^{3/2} (8 x+1)}{8 \sqrt{2 \pi }} \end{aligned}$$

(21)

and with the help of Eq. (21) the Bessel functions in Eq. (20) can be replaced by the right-hand side of Eq. (21). Furthermore, the infinite series in the formula (20) can be summed up in a closed form using

$$\begin{aligned} \zeta (x)=\sum _{p=1}^{\infty }\left( \frac{1}{p}\right) ^x \end{aligned}$$

(22)

where $\zeta (x)$ refers to the Riemann zeta function [37]. We obtain

$$\begin{aligned}&\sum _{p=1}^{\infty } z_0^{-p} I_0^2(2p\beta _\mathrm{c} t_\parallel )I_0( 2p\beta _\mathrm{c} t_\perp ) = \frac{\sqrt{\frac{1}{ t_\perp \beta _\mathrm{c} }} \zeta \left( \frac{3}{2}\right) }{8 \pi ^{3/2} t_\parallel \beta _\mathrm{c} } \nonumber \\&\quad +\frac{\sqrt{\frac{1}{ t_\perp \beta _\mathrm{c} }} \zeta \left( \frac{5}{2}\right) }{64 \pi ^{3/2} t_\parallel ^2 \beta _\mathrm{c} ^2}+\frac{\sqrt{\frac{1}{ t_\perp \beta _\mathrm{c} }} \zeta \left( \frac{5}{2}\right) }{128 \pi ^{3/2} t_\perp t_\parallel \beta _\mathrm{c} ^2} +\frac{\sqrt{\frac{1}{ t_\perp \beta _\mathrm{c} }} \zeta \left( \frac{7}{2}\right) }{2048 \pi ^{3/2} t_\parallel ^3 \beta _\mathrm{c} ^3} \nonumber \\&\quad +\frac{\sqrt{\frac{1}{ t_\perp \beta _\mathrm{c} }} \zeta \left( \frac{7}{2}\right) }{1024 \pi ^{3/2} t_\perp t_\parallel ^2 \beta _\mathrm{c} ^3}+\frac{\sqrt{\frac{1}{ t_\perp \beta _\mathrm{c} }} \zeta \left( \frac{9}{2}\right) }{32768 \pi ^{3/2} t_\perp t_\parallel ^3 \beta _\mathrm{c} ^4} =n, \end{aligned}$$

(23)

Closed-form, analytic solution for the critical temperature $T_\mathrm{c}(n)$ can be obtained by solving the above equation for $n(\beta _\mathrm{c})$ with respect to T by means of the power series expansion. The final result reads

$$\begin{aligned}&k_\mathrm{B}T_\mathrm{c}(n)=\frac{4 n^{2/3} \pi t_\perp \left( \frac{ t_\parallel }{ t_\perp }\right) ^{2/3}}{\zeta \left( \frac{3}{2}\right) ^{2/3}}+n^{4/3} \left( -\frac{4 \pi ^2 t_\perp \root 3 \of {\frac{ t_\parallel }{ t_\perp }} \zeta \left( \frac{5}{2}\right) }{3 \zeta \left( \frac{3}{2}\right) ^{7/3}}-\frac{2 \pi ^2 t_\parallel \root 3 \of {\frac{ t_\parallel }{ t_\perp }} \zeta \left( \frac{5}{2}\right) }{3 \zeta \left( \frac{3}{2}\right) ^{7/3}}\right) \nonumber \\&\quad +n^2 \left( \frac{\pi ^3 t_\perp \zeta \left( \frac{5}{2}\right) ^2}{\zeta \left( \frac{3}{2}\right) ^4} \right. \left. +\frac{\pi ^3 t_\parallel \zeta \left( \frac{5}{2}\right) ^2}{\zeta \left( \frac{3}{2}\right) ^4}+\frac{\pi ^3 t_\parallel ^2 \zeta \left( \frac{5}{2}\right) ^2}{4 t_\perp \zeta \left( \frac{3}{2}\right) ^4}-\frac{\pi ^3 t_\perp \zeta \left( \frac{7}{2}\right) }{6 \zeta \left( \frac{3}{2}\right) ^3}-\frac{\pi ^3 t_\parallel \zeta \left( \frac{7}{2}\right) }{3 \zeta \left( \frac{3}{2}\right) ^3}\right) \nonumber \\&\quad +n^{8/3} \left( -\frac{77 \pi ^4 t_\perp \left( \frac{ t_\parallel }{ t_\perp }\right) ^{2/3} \zeta \left( \frac{5}{2}\right) ^3}{54 \zeta \left( \frac{3}{2}\right) ^{17/3}}-\frac{77 \pi ^4 t_\perp ^2 \left( \frac{ t_\parallel }{ t_\perp }\right) ^{2/3} \zeta \left( \frac{5}{2}\right) ^3}{81 t_\parallel \zeta \left( \frac{3}{2}\right) ^{17/3}} \right. \left. -\frac{77 \pi ^4 t_\parallel \left( \frac{ t_\parallel }{ t_\perp }\right) ^{2/3} \zeta \left( \frac{5}{2}\right) ^3}{108 \zeta \left( \frac{3}{2}\right) ^{17/3}}\right. \nonumber \\&\quad -\frac{77 \pi ^4 t_\parallel ^2 \left( \frac{ t_\parallel }{ t_\perp }\right) ^{2/3} \zeta \left( \frac{5}{2}\right) ^3}{648 t_\perp \zeta \left( \frac{3}{2}\right) ^{17/3}} +\frac{55 \pi ^4 t_\perp \left( \frac{ t_\parallel }{ t_\perp }\right) ^{2/3} \zeta \left( \frac{5}{2}\right) \zeta \left( \frac{7}{2}\right) }{72 \zeta \left( \frac{3}{2}\right) ^{14/3}} +\left. \frac{11 \pi ^4 t_\perp ^2 \left( \frac{ t_\parallel }{ t_\perp }\right) ^{2/3} \zeta \left( \frac{5}{2}\right) \zeta \left( \frac{7}{2}\right) }{36 t_\parallel \zeta \left( \frac{3}{2}\right) ^{14/3}} \right. \nonumber \\&\quad \left. +\frac{11 \pi ^4 t_\parallel \left( \frac{ t_\parallel }{ t_\perp }\right) ^{2/3} \zeta \left( \frac{5}{2}\right) \zeta \left( \frac{7}{2}\right) }{36 \zeta \left( \frac{3}{2}\right) ^{14/3}}-\frac{\pi ^4 t_\perp \left( \frac{ t_\parallel }{ t_\perp }\right) ^{2/3} \zeta \left( \frac{9}{2}\right) }{24 \zeta \left( \frac{3}{2}\right) ^{11/3}}\right) +O\left( n^{10/3}\right) . \end{aligned}$$

(24)

To have better insight into the density dependence of the critical temperature, we present the outcome in Eq. (24) with numerically evaluated coefficients of the expansion.

$$\begin{aligned} k_\mathrm{B}T_\mathrm{c}(n)= & {} 6.625 n^{2/3} t_\perp \left( \frac{ t_\parallel }{ t_\perp }\right) ^{2/3} +n^2 \left( 0.871467 t_\perp +0.54487 t_\parallel +\frac{0.299516 t_\parallel ^2}{ t_\perp }\right) \nonumber \\&+n^{4/3} \left( -1.8782 t_\perp \root 3 \of {\frac{ t_\parallel }{ t_\perp }}-0.9391 t_\parallel \root 3 \of {\frac{ t_\parallel }{ t_\perp }}\right) \nonumber \\&+n^{8/3} \left[ -0.306424 t_\perp \left( \frac{ t_\parallel }{ t_\perp }\right) ^{2/3}-\frac{0.459391 t_\perp ^2 \left( \frac{ t_\parallel }{ t_\perp }\right) ^{2/3}}{ t_\parallel } \right. \nonumber \\&\left. -0.217231 t_\parallel \left( \frac{ t_\parallel }{ t_\perp }\right) ^{2/3}-\frac{0.12108 t_\parallel ^2 \left( \frac{ t_\parallel }{ t_\perp }\right) ^{2/3}}{ t_\perp }\right] +O\left( n^{10/3}\right) . \end{aligned}$$

(25)

For the isotropic 3D system ( $t_\perp =t_\parallel \equiv t$), the formula for $T_\mathrm{c}$ simplifies to

$$\begin{aligned}&\frac{k_\mathrm{B} T_\mathrm{c}(n)}{t}=\frac{4 n^{2/3} \pi t}{\zeta \left( \frac{3}{2}\right) ^{2/3}}-\frac{2 n^{4/3} \pi ^2 t \zeta \left( \frac{5}{2}\right) }{\zeta \left( \frac{3}{2}\right) ^{7/3}} \nonumber \\&\quad +n^2 \left( \frac{9 \pi ^3 t \zeta \left( \frac{5}{2}\right) ^2}{4 \zeta \left( \frac{3}{2}\right) ^4}-\frac{\pi ^3 t \zeta \left( \frac{7}{2}\right) }{2 \zeta \left( \frac{3}{2}\right) ^3}\right) +n^{8/3} \left( -\frac{77 \pi ^4 t \zeta \left( \frac{5}{2}\right) ^3}{24 \zeta \left( \frac{3}{2}\right) ^{17/3}} \right. \nonumber \\&\quad \left. +\frac{11 \pi ^4 t \zeta \left( \frac{5}{2}\right) \zeta \left( \frac{7}{2}\right) }{8 \zeta \left( \frac{3}{2}\right) ^{14/3}}-\frac{\pi ^4 t \zeta \left( \frac{9}{2}\right) }{24 \zeta \left( \frac{3}{2}\right) ^{11/3}}\right) +O\left( n^{10/3}\right) , \end{aligned}$$

(26)

or in terms of numerical values for coefficients

$$\begin{aligned} \frac{k_\mathrm{B}T_\mathrm{c}(n)}{t}= & {} 6.625 n^{2/3}-2.8173 n^{4/3} +1.71585 n^2 \nonumber \\&-1.10413 n^{8/3}+O\left( n^{10/3}\right) . \end{aligned}$$

(27)

To examine the variation of the critical temperature of the condensation $T_\mathrm{c}$ as a function of lattice anisotropy with respect to the 3D isotropic system as a frame of reference, we set in Eq. (5) the hopping parameters as $t_\perp =2t\eta$ and $t_\parallel =2t(1-\eta )$ with the dimensionless parameter $0\le \eta \le 1$, where $\eta =t_\perp /(t_\perp +t_\parallel )$, so that

$$\begin{aligned} t(\mathbf{k)}=4t\{(1-\eta )[\cos (k_x a)+\cos (k_y a)]+ \eta \cos (k_z a)\}. \end{aligned}$$

(28)

In this way, we can smoothly interpolate between the case of decoupled planes ($\eta =0$), isotropic cubic lattice ($\eta =1/2$) and isolated chains ($\eta =1$). In Fig. 1, we plot $T_\mathrm{c}$ as a function of the parameter $\eta$; the critical temperature reaches its maximum for isotropic cubic lattice, while it drops to zero for 1D or 2D uniform case according to the observation that in a uniform gas, Bose–Einstein condensation cannot occur at a finite temperature in one or two dimensions because thermal fluctuations would destabilize the condensate. When the in-plane interactions are reduced, the critical temperature drops quite rapidly as a function of $t_\parallel /t_\perp$ compared to the isotropic system $(t_\parallel =t_\perp )$ and for $t_\parallel /t_\perp =10^{-5}$ (weakly coupled chains) is of 0.2% of the isotropic value of $T_\mathrm{c}$. However, in the case of reduced coupling along the c direction (weakly coupled planes), even for the anisotropy parameter as small as $t_\perp /t_\parallel =10^{-6}$ the critical temperature is still about 16% of the critical temperature for the isotropic 3D system. It appears that this feature of the 2D$\rightarrow$3D$\rightarrow$1D critical temperature crossover behavior may have interesting consequences, especially for physics of multilayered systems (Fig. 2).

3.2 Condensate Fraction

The condensate density $n_0$ of a BEC is a measure for the off-diagonal long-range order. The growth of the condensate fraction $\psi =n_0/n$ below $T_\mathrm{c}$ can be written by combining Eqs.(13) and (20) as

$$\begin{aligned} \frac{n_0}{n}=1- \frac{1}{n} \sum _{p=1}^{\infty } \frac{1}{z_0^p} I_0^2(2p\beta t_\parallel )I_0( 2p\beta t_\perp ). \end{aligned}$$

(29)

Making use of the asymptotic form of $I_0(x)$ and the power series expansion for the Bessel functions the infinite sum in Eq. (29) can be performed analytically. Furthermore, inverting the $T_\mathrm{c}(n)$ function in Eq. (24) by means of series expansion and substituting the resulting $n(T_\mathrm{c})$ dependence into Eq. (29), we arrive at order parameter $\psi$ of the form

$$\begin{aligned} \psi\simeq & {} 1-\left( \frac{T}{T_\mathrm{c}}\right) ^{3/2} \nonumber \\&\times \frac{1+\frac{ k_\mathrm{B} T \zeta \left( \frac{5}{2}\right) }{16 t_\perp \zeta \left( \frac{3}{2}\right) }+\frac{ k_\mathrm{B} T \zeta \left( \frac{5}{2}\right) }{8 t_\parallel \zeta \left( \frac{3}{2}\right) }+\frac{ k_\mathrm{B}^2 T^2 \zeta \left( \frac{7}{2}\right) }{256 t_\parallel ^2 \zeta \left( \frac{3}{2}\right) }+\frac{ k_\mathrm{B}^2 T^2 \zeta \left( \frac{7}{2}\right) }{128 t_\perp t_\parallel \zeta \left( \frac{3}{2}\right) }+\frac{ k_\mathrm{B}^3 T^3 \zeta \left( \frac{9}{2}\right) }{4096 t_\perp t_\parallel ^2 \zeta \left( \frac{3}{2}\right) }}{1+\frac{ k_\mathrm{B} T_\mathrm{c} \zeta \left( \frac{5}{2}\right) }{16 t_\perp \zeta \left( \frac{3}{2}\right) }+\frac{ k_\mathrm{B} T_\mathrm{c} \zeta \left( \frac{5}{2}\right) }{8 t_\parallel \zeta \left( \frac{3}{2}\right) }+\frac{ k_\mathrm{B}^2 T_\mathrm{c} ^2 \zeta \left( \frac{7}{2}\right) }{256 t_\parallel ^2 \zeta \left( \frac{3}{2}\right) }+\frac{ k_\mathrm{B}^2 T_\mathrm{c} ^2 \zeta \left( \frac{7}{2}\right) }{128 t_\perp t_\parallel \zeta \left( \frac{3}{2}\right) }+\frac{ k_\mathrm{B}^3 T_\mathrm{c} ^3 \zeta \left( \frac{9}{2}\right) }{4096 t_\perp t_\parallel ^2 \zeta \left( \frac{3}{2}\right) }}, \end{aligned}$$

(30)

which explicitly shows the modifications due the anisotropic lattice of the standard outcome for $\psi$ in the uniform BEC in 3D: $\psi =1-(T/T_\mathrm{c})^{3/2}$. For the isotropic case $t\equiv t_\perp =t_\parallel$, Eq. (30) reduces to

$$\begin{aligned} \frac{n_0}{n}\simeq 1-\left( \frac{T}{T_\mathrm{c}}\right) ^{\frac{3}{2}}\frac{1+\frac{3 k_\mathrm{B} T \zeta \left( \frac{5}{2}\right) }{16 t \zeta \left( \frac{3}{2}\right) }+\frac{3 k_\mathrm{B}^2 T^2 \zeta \left( \frac{7}{2}\right) }{256 t^2 \zeta \left( \frac{3}{2}\right) }+\frac{k_\mathrm{B}^3 T^3 \zeta \left( \frac{9}{2}\right) }{4096 t^3 \zeta \left( \frac{3}{2}\right) }}{1+\frac{3 k_\mathrm{B} T_\mathrm{c} \zeta \left( \frac{5}{2}\right) }{16 t \zeta \left( \frac{3}{2}\right) }+\frac{3 k_\mathrm{B}^2 T_\mathrm{c}^2 \zeta \left( \frac{7}{2}\right) }{256 t^2 \zeta \left( \frac{3}{2}\right) }+\frac{ k_\mathrm{B}^3 T_\mathrm{c}^3 \zeta \left( \frac{9}{2}\right) }{4096 t^3 \zeta \left( \frac{3}{2}\right) }} \end{aligned}$$

(31)

The plots in Fig. 3 show the superfluid fraction $n_0/n$ as a function of temperature for various anisotropies given by $t_\perp /t_\parallel$. As can be seen, plots of $\psi$ as a function of normalized temperature $T/T_\mathrm{c}$ fall approximately into the same curve, irrespective of the value of $t_\perp /t_\parallel$. This behavior indicates that $\psi (T)$ roughly behaves as a homogeneous function $\psi (T)\approx f(T/T_\mathrm{c})$—a feature that may be traced back in formula given by Eq. (30). It is known that the non-interacting continuous Bose model exhibits a scaling behavior [38, 39], where all the thermodynamic quantities can be expressed in terms of the critical temperature $T_\mathrm{c}$ and the reduced temperature $T/T_\mathrm{c}$, while the dimensionless quantities, like the condensate fraction $\psi$ or the entropy per particle, will depend only on $T/T_\mathrm{c}$. Therefore, for a lattice version of BEC the scaling behavior persists only approximately. However, it is desirable while deriving the temperature dependence of various thermodynamic quantities to study them as a function of $T/T_\mathrm{c}$ as well as of $T/t(\mathbf{0})$ in order to explore the role of the characteristic energy scales of the system. We will follow this prescription in the remainder of the paper.

Fig. 3

Order parameter $\psi$ representing the condensate fraction as a function of temperature T: In the left column, T is normalized to the critical temperature $T_\mathrm{c}$, whereas in the right column the normalization is given by the bandwidth energy $t(\mathbf{0})$. Panels in the top row show the coupled planes case (with the anisotropy parameters indicated); the insets depict occupation number n dependence of $\psi$ for the isotropic case ($t_\perp =t_\parallel$). In the bottom row, the plots are for the parameters corresponding to the coupled chains case. On main panels, the boson occupation number is fixed to $n=1$ (Color figure online)

3.3 Chemical Potential

The nontrivial aspect in the evaluation of thermodynamic properties of BEC is the determination of the chemical potential $\mu$, which describes the change in energy when one particle is added to the system at constant entropy and volume. In the context of a many-body systems, $\mu$ is the minimum energy required to add a particle to a system in thermal equilibrium with the reservoir. This condition requires that the grand potential in Eq. (8) reaches the minimum. Once the chemical potential is known, all the thermodynamic characteristics like total energy entropy and heat capacity follow directly via sums over the energy levels involving the occupation numbers. In Fig. 4, we plot numerically evaluated temperature dependence of $\mu$. It is interesting to observe the (quasi) scaling behavior of the chemical potential. When $\mu (T)$ is plotted against the bandwidth energy normalized temperature all $\mu (T)$ plots falls on the same curve for $T<T_\mathrm{c}$—which is a direct consequence of Eq. (15). Interestingly, a similar behavior persists (approximately) also in the high-temperature regime $T>T_\mathrm{c}$ for a broad range of lattice anisotropies. The evolution of the chemical potential $\mu$ as a function of number of bosons per lattice site n and temperature T is given in Fig. 5. For the BEC of quasi-particles, the chemical potential of confined magnons in a ferromagnetic nano-structured material has been measured to reflect the growth of the magnon density and to track the critical temperature of BEC in these materials [28].

Fig. 4

Chemical potential $\mu$ as a function of the temperature T for coupled planes case for $t_\parallel > t_\perp$ (upper row panels). Plots in the lower row corresponds to the coupled chain case . For plots in the left column, T is normalized to the in plane (chain) hoping amplitude temperature $t_\parallel$ ($t_\perp$), whereas for the plots in the right column the normalization is given by the bandwidth energy $t(\mathbf{0})$. The boson occupation number is fixed to $n=1$ (Color figure online)

Fig. 5

Evolution of the chemical potential $\mu$ as a function of number of bosons per lattice site n and temperature T normalized to the bandwidth energy $t(\mathbf{0})$ (Color figure online)

4 Thermodynamics

4.1 Internal Energy

With the definition in Eq. (9), the energy of the macroscopically occupied state (the ground state) at $\mathbf{k}=\mathbf{0}$ is taken to be zero. Therefore, only the excited states contribute to the total energy of the system. As a consequence, the internal energy per lattice site of the system $\bar{E}=E/N$ (which is useful for studying, e.g., isentropic processes), where $E=-(\partial /\partial \beta )\ln Z$, is given by

$$\begin{aligned} {\bar{E}}=\frac{1}{N}\sum _{\mathbf{k}}^N \frac{\mathcal{E}(\mathbf{k)}}{\exp \left[ \beta \left( \mathcal{E}(\mathbf{k)} -\bar{\mu }\right) \right] -1} \end{aligned}$$

(32)

Below $T_\mathrm{c}$, the shifted chemical potential $\bar{\mu }$ vanishes. By using the result of integral in Eq.  (23), the internal energy can be given in a form of the low-temperature expansion

$$\begin{aligned} {\bar{E}}= & {} \frac{3 k_\mathrm{B} ^{5/2} T^{5/2} \zeta \left( \frac{5}{2}\right) }{16 \pi ^{3/2} \sqrt{ t_\perp } t_\parallel }+\frac{ k_\mathrm{B} ^{7/2} T^{7/2} \zeta \left( \frac{7}{2}\right) }{64 \pi ^{3/2} \sqrt{ t_\perp } t_\parallel ^2} \nonumber \\&+ \frac{ k_\mathrm{B} ^{7/2} T^{7/2} \zeta \left( \frac{7}{2}\right) }{128 \pi ^{3/2} t_\perp ^{3/2} t_\parallel }+\frac{ k_\mathrm{B} ^{9/2} T^{9/2} \zeta \left( \frac{9}{2}\right) }{4096 \pi ^{3/2} \sqrt{ t_\perp } t_\parallel ^3} \nonumber \\&+\frac{ k_\mathrm{B} ^{9/2} T^{9/2} \zeta \left( \frac{9}{2}\right) }{2048 \pi ^{3/2} t_\perp ^{3/2} t_\parallel ^2}. \end{aligned}$$

(33)

When calculating $\bar{E}$ for $T>T_\mathrm{c}$ for fixed density n, we must combine Eq. (32) with equation for the chemical potential (11). Figures 6 and 7 show the internal energy $\bar{E}$ as a function of temperature (in units of $T_\mathrm{c}$) as well as of the bandwidth energy $t(\mathbf{0})$ combined with changes of anisotropy parameters and various occupation numbers.

Fig. 6

Internal energy ${\bar{E}}$ as a function of temperature normalized to the critical temperature $T_\mathrm{c}$ and the bandwidth energy $t(\mathbf{0})$ for various anisotropy parameters (as indicated in the plots). The boson occupation number is fixed to $n=1$ (Color figure online)

Fig. 7

Internal energy ${\bar{E}}$ as a function of normalized temperature (with respect to the critical temperature $T_\mathrm{c}$ and the bandwidth energy $t(\mathbf{0})$) for several values of occupation parameter n and anisotropy parameters (Color figure online)

4.2 Entropy

In thermodynamics, entropy S is commonly understood as a measure of disorder. The entropy S per lattice site can be determined from the thermodynamic potential as

$$\begin{aligned} S=\left( -\frac{\partial \bar{\Omega }}{\partial T}\right) _V. \end{aligned}$$

(34)

Accordingly, with the use of the lattice density of states (DOS), $\rho (\varepsilon )$, (see, “Appendix” for technical details) the 3D lattice integrals over momenta can be reduced to a single with the result

$$\begin{aligned} \frac{S}{k_\mathrm{B}}= & {} -\frac{\mu }{k_\mathrm{B}T} -\frac{1}{k_\mathrm{B}T} \int {\hbox {d}}\varepsilon \rho (\varepsilon ) \left( \varepsilon + \frac{\varepsilon +\mu }{e^{\frac{\varepsilon +\mu }{k_\mathrm{B}T}}-1} \right) \nonumber \\&- \int {\hbox {d}}\varepsilon \rho (\varepsilon ) \ln \left( 1- e^{\frac{\varepsilon +\mu }{k_\mathrm{B}T}} \right) . \end{aligned}$$

(35)

Fig. 8

Entropy S as a function of normalized temperature (with respect to the critical temperature $T_\mathrm{c}$ and the bandwidth energy $t(\mathbf{0})$) for several values of anisotropy parameters. Panels in the top row shows the coupled planes case with $t_\perp < t _\parallel$. In the lower row of panels, analogous plots for coupled chains scenario with $t_\parallel > t _\perp$. In each panel, the entropy S for the the isotropic case $t_\parallel = t _\perp$ is also depicted for comparison. The boson occupation number is fixed to $n=1$ (Color figure online)

In the low-temperature condensed phase, for $T<T_\mathrm{c}$ with fixed chemical potential, the entropy can be conveniently expressed in a form of a power series:

$$\begin{aligned} S= & {} \frac{5 k_\mathrm{B} ^{5/2} T^{3/2} \zeta \left( \frac{5}{2}\right) }{16 \pi ^{3/2} \sqrt{ t_\perp } t_\parallel }+\frac{7 k_\mathrm{B} ^{7/2} T^{5/2} (2 t_\perp + t_\parallel ) \zeta \left( \frac{7}{2}\right) }{640 \pi ^{3/2} t_\perp ^{3/2} t_\parallel ^2} \nonumber \\&+\frac{9 k_\mathrm{B} ^{9/2} T^{7/2} ( t_\perp +2 t_\parallel ) \zeta \left( \frac{9}{2}\right) }{28672 \pi ^{3/2} t_\perp ^{3/2} t_\parallel ^3 } +O\left( T^{9/2}\right) . \end{aligned}$$

(36)

Figures 8 and 9 show that the temperature dependence of the entropy S for various anisotropies and boson occupation numbers. The entropy curve monotonically decreases with decreasing temperature, and the drop-off of the entropy curve toward $T=0$ is in accordance with the third law of thermodynamics. It is interesting to note that for fixed density n, the curves corresponding to S as a function of the reduced temperature $T/t(\mathbf{0})$ to various anisotropy parameters lie closely together, exhibiting approximate scaling behavior—a feature mentioned while discussing the behavior of the chemical potential $\mu$. However, in contrast to the chemical potential, where the quasi-scaling was observed as a function of $T/T_\mathrm{c}$, the entropy underlines the role of the bandwidth $t(\mathbf{0})$ as the second relevant natural energy scale of the system.

Fig. 9

Entropy S as a function of normalized temperature (with respect to the critical temperature $T_\mathrm{c}$ and the bandwidth energy $t(\mathbf{0})$) for several values of boson occupation parameter n and anisotropy parameters (as indicated in the plots) (Color figure online)

4.3 Specific Heat

The heat capacity is one of the important quantities enclosing information about the nature of a Bose condensation phase transition, and predictions of the heat capacity for trapped atomic gases were made [39,40,41,42,43], and special techniques have been implemented to measure it in BEC [44]. The specific heat also gives access to other thermodynamic quantities, as the entropy, which can be obtained by the integral $S(T)= \int _0^T C_V(T')\hbox {d}T'$. The constant-volume specific heat $C_V$ is defined as

$$\begin{aligned} \frac{C_V}{k_\mathrm{B}}=\frac{1}{k_\mathrm{B}}\left( \frac{\partial {\bar{E}}}{\partial T}\right) _{\mathcal{N},V}. \end{aligned}$$

(37)

In order to understand in more detail what happens as we pass through the critical temperature $T_\mathrm{c}$, we will focus on how the heat capacity behaves on either side of the critical temperature. At constant particle number $\mathcal{N}$, the chemical potential $\mu$ becomes an implicit function of temperature T and must be determined accordingly from Eq. (11). The result is given by

Fig. 10

Plots of the specific heat $C_V$ as a function of normalized temperature (with respect to the critical temperature $T_\mathrm{c}$ and the bandwidth energy $t(\mathbf{0})$) for several values of anisotropy parameters. The boson occupation number is fixed to $n=1$ (Color figure online)

Fig. 11

Plots of the he specific heat $C_V$ as a function of normalized temperature (with respect to the critical temperature $T_\mathrm{c}$ and the bandwidth energy $t(\mathbf{0})$) for the several values of occupation parameter n. Panels in the top row show the isotropic case $t _\parallel / t_\perp =1$, in the mid-row, the coupled chains scenario $t _\parallel / t_\perp <1$, in the bottom row the coupled planes case with $t_\perp / t _\parallel >1$ (Color figure online)

Fig. 12

Plots showing the contrasting behavior of the specific heat $C_V$ for large values of 3D lattice anisotropy. In the right column, panels with plots for the coupled chains with $t _\parallel / t_\perp \ll 1$. In the left column, plots depicting the coupled planes scenario with $t_\perp /t _\parallel \ll 1$, (The actual values for $t _\parallel /t_\perp$ and $t_\perp /t _\parallel$ are indicated in each plot.) The boson occupation number is fixed to $n=1$ (Color figure online)

$$\begin{aligned} C_V=\frac{1}{{4 k_\mathrm{B} T^2}} \int \hbox {d} \varepsilon \rho (\varepsilon ) \frac{(\epsilon -t(\mathbf{0})) \left( \varepsilon +\mu (T)-T \mu '(T)\right) }{\sinh ^2\left( \frac{\varepsilon +\mu (T)}{2 k_\mathrm{B} T}\right) } , \end{aligned}$$

(38)

while the temperature derivative of the chemical potential calculated from Eq.  (11) is given by

$$\begin{aligned} \mu '(T) =\frac{1}{T} \left[ \frac{\int {\hbox {d}}\varepsilon \rho (\varepsilon ) \varepsilon \, \mathrm{csch}^2\left( \frac{\varepsilon +\mu (T)}{2 k_\mathrm{B} T}\right) }{ \int {\hbox {d}}\varepsilon \rho (\varepsilon ) \, \mathrm{csch}^2\left( \frac{\varepsilon +\mu (T)}{2 k_\mathrm{B} T}\right) }+\mu (T) \right] . \end{aligned}$$

(39)

Furthermore, we observe that $\mu (T_\mathrm{c}^+)=\mu (T_\mathrm{c}^-)\equiv -t(\mathbf{0})$ (see, Eq. 15) and as a consequence of constant chemical potential within the condensed phase we have $\mu '(T_\mathrm{c}^-)=0$. Defining the difference on both sides of the phase transition as

$$\begin{aligned} \Delta C_V\equiv C_V(T_\mathrm{c}^+)-C_V(T_\mathrm{c}^-), \end{aligned}$$

(40)

we obtain from Eqs.(38) and (39)

$$\begin{aligned} \Delta C_v=\frac{\mu '(T_\mathrm{c}^+)}{{4 k_\mathrm{B} T_\mathrm{c}}} \int {\hbox {d}} \varepsilon \rho (\varepsilon ) \frac{t(\mathbf{0})-\varepsilon }{\sinh ^2\left( \frac{t(\mathbf{0})-\varepsilon }{2k_\mathrm{B} T_\mathrm{c}}\right) } \end{aligned}$$

(41)

with

$$\begin{aligned} \mu '(T_\mathrm{c}^+) =\frac{1}{T_\mathrm{c}} \left[ \frac{\int {\hbox {d}}\varepsilon \rho (\varepsilon ) \varepsilon \, \mathrm{csch}^2\left( \frac{\varepsilon -t(\mathbf{0})}{2 k_\mathrm{B} T_\mathrm{c}}\right) }{ \int {\hbox {d}}\varepsilon \rho (\varepsilon ) \, \mathrm{csch}^2\left( \frac{\varepsilon -t(\mathbf{0})}{2 k_\mathrm{B} T_\mathrm{c}}\right) }-t(\mathbf{0}) \right] . \end{aligned}$$

(42)

which follows from Eq. (39). After some simple algebra, the expression for specific heat jump at $T_\mathrm{c}$ in Eq. (38) can be manipulated to the form

$$\begin{aligned} \Delta C_V=-\frac{1}{{4 k_\mathrm{B} T_\mathrm{c}^2}}\frac{\left\{ \int {\hbox {d}}\varepsilon \rho (\varepsilon )\left[ \varepsilon -t(\mathbf{0})\right] \, \mathrm{csch}^2\left( \frac{\varepsilon -t(\mathbf{0})}{2 k_\mathrm{B} T_\mathrm{c}}\right) \right\} ^2}{ \int {\hbox {d}}\varepsilon \rho (\varepsilon ) \, \mathrm{csch}^2\left( \frac{\varepsilon -t(\mathbf{0})}{2 k_\mathrm{B} T_\mathrm{c}}\right) }. \end{aligned}$$

(43)

Because the integral in the numerator of Eq. (43) stays finite while the one in the denominator diverges, there is no jump discontinuity in the specific heat. $C_V$ is exhibiting merely a cusp in its temperature dependence at $T_\mathrm{c}$. The main finding is that in the thermodynamic limit, at critical temperature $T_\mathrm{c}$ the heat capacity at constant number shows a cusp singularity, which is analogous to the $\lambda$-transition of liquid helium [45]. The low-temperature expansion for $T<T_\mathrm{c}$ of the heat capacity is given by

$$\begin{aligned} C_V= & {} \frac{15 k_\mathrm{B}^{5/2} T^{3/2} \zeta \left( \frac{5}{2}\right) }{32 \pi ^{3/2} \sqrt{t_\perp } t_\parallel }+\frac{7 k_\mathrm{B}^{7/2} T^{5/2} \zeta \left( \frac{7}{2}\right) }{128 \pi ^{3/2} \sqrt{t_\perp } t_\parallel ^2} \nonumber \\&+ \frac{7 k_\mathrm{B}^{7/2} T^{5/2} \zeta \left( \frac{7}{2}\right) }{256 \pi ^{3/2} t_\perp ^{3/2} t_\parallel } +\frac{9 k_\mathrm{B}^{9/2} T^{7/2} \zeta \left( \frac{9}{2}\right) }{8192 \pi ^{3/2} \sqrt{t_\perp } t_\parallel ^3} \nonumber \\&+\frac{9 k_\mathrm{B}^{9/2} T^{7/2} \zeta \left( \frac{9}{2}\right) }{4096 \pi ^{3/2} t_\perp ^{3/2} t_\parallel ^2}+O\left( T^{9/2}\right) . \end{aligned}$$

(44)

Numerical evaluation of the coefficients in the series results in

$$\begin{aligned} C_V= & {} \frac{0.000208094 k_\mathrm{B} ^{9/2} T^{7/2}}{\sqrt{ t_\perp } t_\parallel ^3}+\frac{0.000416188 k_\mathrm{B} ^{9/2} T^{7/2}}{ t_\perp ^{3/2} t_\parallel ^2} \nonumber \\&+ \frac{0.0110658 k_\mathrm{B} ^{7/2} T^{5/2}}{\sqrt{ t_\perp } t_\parallel ^2}+\frac{0.00553292 k_\mathrm{B} ^{7/2} T^{5/2}}{ t_\perp ^{3/2} t_\parallel } \nonumber \\&+\frac{0.112928 k_\mathrm{B} ^{5/2} T^{3/2}}{\sqrt{ t_\perp } t_\parallel }+O\left( T^{9/2}\right) . \end{aligned}$$

(45)

For the isotropic 3D case ($t_\perp =t_\parallel$), one obtains from Eq. (44)

$$\begin{aligned} C_V= & {} \frac{15 k_\mathrm{B} ^{5/2} T^{3/2} \zeta \left( \frac{5}{2}\right) }{32 \pi ^{3/2} t^{3/2}}+\frac{21 k_\mathrm{B} ^{7/2} T^{5/2} \zeta \left( \frac{7}{2}\right) }{256 \pi ^{3/2} t^{5/2}} \nonumber \\&+\frac{27 k_\mathrm{B} ^{9/2} T^{7/2} \zeta \left( \frac{9}{2}\right) }{8192 \pi ^{3/2} t^{7/2}} +O\left( T^{9/2}\right) , \end{aligned}$$

(46)

while the numerical evaluation for the above series gives

$$\begin{aligned} C_V= & {} \frac{0.112928 k_\mathrm{B} ^{5/2} T^{3/2}}{t^{3/2}}+\frac{0.0165988 k_\mathrm{B} ^{7/2} T^{5/2}}{t^{5/2}} \nonumber \\&+\frac{0.000624282 k_\mathrm{B} ^{9/2} T^{7/2}}{t^{7/2}}+O\left( T^{9/2}\right) . \end{aligned}$$

(47)

In Fig. 10, we plot the normalized heat capacity versus the reduced temperature (with respect to the critical temperature $T_\mathrm{c}$ and the bandwidth energy $t(\mathbf{0})$), for various anisotropies of the lattice. Again we note that for fixed density n and below $T_\mathrm{c}$, the curves for specific heat $C_V$, as a function of the reduced temperature $T/t(\mathbf{0})$, for various anisotropy parameters lie closely together, exhibiting approximate scaling behavior—a feature mentioned while discussing the behavior of the chemical potential $\mu$ and entropy S. However, in contrast to the chemical potential, where the quasi-scaling was observed as a function of $T/T_\mathrm{c}$, the specific heat underlines the role of the bandwidth $t(\mathbf{0})$ as the relevant energy scale of the system.

Figure 11 depicts the constant-volume specific heat for the several values of occupation parameter n. The lowering of the critical temperature for decreasing number of particle density is due to the fact that a system with smaller n has a larger available effective volume [40]. It is interesting to investigate the behavior of the specific heat in the extreme cases of the nearly isolated planes ($t_\perp /t_\parallel \ll 1$) or chains ($t_\parallel /t_\perp \ll 1$). In Fig. 12, we have presented the contrasting behavior of the specific heat $C_V$ for large values of 3D lattice anisotropy. While for coupled planes case, the behavior of $C_V$ is consistent in a broad range of anisotropies showing a sharp peak at the critical temperature $T_\mathrm{c}$, the large anisotropy coupled chains scenario looks quite different. Namely, the specific heat shows a characteristic “peak-dip-hump” feature (see Fig. 12). Starting from below $T_\mathrm{c}$ the specific heat raises rapidly exhibiting a sharp peak at the critical temperature $T_\mathrm{c}$. Subsequently, it diminishes in a narrow temperature range forming a small dip and rises again exhibiting a broad hump. It is interesting to note that a similar $C_V(T)$ dependence occurs in 1D magnetic chains [46], which might suggest that condensation of magnons dominates the underlying physics (Figs. 13, 14).

4.4 Compressibility

Thermodynamic quantities like isothermal compressibility are expected to become singular at the second-order phase transition. Now calculate the isothermal compressibility at constant number of particles, $\kappa _T$ defined as:

$$\begin{aligned} \kappa _T =-\frac{1}{V} \left( \frac{\partial V}{\partial p} \right) _{\mathcal{N},T} \end{aligned}$$

(48)

equivalently, using the relation $\Omega =-pV$, Eq. (48) can be manipulated into the form

$$\begin{aligned} \kappa _T =-\frac{1}{V} \frac{V^2}{\mathcal{N}^2}\left( \frac{\partial \mathcal{N}}{\partial \mu } \right) _{V,T} \end{aligned}$$

(49)

involving the differentiation with respect to the chemical potential instead of pressure. Furthermore, relating the volume V to the number of particles and introducing the density of bosons n via the relations $V=a^3N$, $\mathcal{N} = nN$ with $a^3$ being the volume of the unit lattice cell we arrive at

$$\begin{aligned} \kappa _T =- \frac{a^3}{n^2}\left( \frac{\partial n}{\partial \mu } \right) _{V,T}, \end{aligned}$$

(50)

or explicitly, using the lattice density of states $\rho (\varepsilon )$ (see, “Appendix”)

$$\begin{aligned} \kappa _T =\frac{ \beta a^3}{n^2} \int \hbox {d} \varepsilon \frac{ \rho (\varepsilon )e^{-\beta (\varepsilon +\mu )} }{\left[ 1-e^{-\beta (\varepsilon +\mu )}\right] ^2} \end{aligned}$$

(51)

Fig. 13

Isothermal compressibility $\kappa _T$ as a function of temperature (normalized to the critical temperature $T_\mathrm{c}$ and the bandwidth energy $t(\mathbf{0})$ for several values of the anisotropy parameters. The upper top row depicts the coupled chains case with $t_\parallel /t _\perp <1$. Panels in the lower row: coupled planes case $t _\perp /t_\parallel <1$ (with corresponding parameters placed in each plot). The boson occupation number is fixed to $n=1$(Color figure online)

Fig. 14

Compressibility $\kappa _T$ as a function of normalized temperature (with respect to the critical temperature $T_\mathrm{c}$ and the bandwidth energy $t(\mathbf{0})$) for several values of occupation parameter n. Plots in the top row show isotropic case $t _\parallel / t_\perp =1$, in the mid-row the coupled chains arrangement with $t _\parallel / t_\perp <1$, on the bottom the coupled planes case with $t_\perp / t _\parallel <1$. The actual values of parameters are indicated in each plot (Color figure online)

For $T<T_\mathrm{c}$, the compressibility is infinite, as it diverges while approaching the critical temperature from above. For high temperature, one finds, as expected, the Curie law-like behavior $\kappa _T\sim 1/T$ resembling the behavior of an ideal classical gas. The importance of the $\kappa _T$ parameter revealing nature of the phase transition of Bose systems has been recognized [47] as the most appropriate susceptibility indicating the phase transition of a harmonically confined Bose gas.

5 Conclusions

In this paper, we presented a theoretical study of condensation of bosons in tight-binding bands corresponding to the anisotropic three-dimensional lattices. The discussion focused on the non-interacting Bose gas where the theoretical predictions can be made without resorting to approximations, so that the connection of Bose–Einstein condensation and superfluidity can be studied at a more fundamental level. However, even in this limit to make the calculations feasible one, an efficient and precise method of handling multi-dimensional momentum lattice integrals is called for. This requires a combination of analytic solution for Watson-type integrals for lattice Green’s functions and numerical evaluation of the corresponding lattice density of states. The chief goal of this paper is to investigate how a periodic lattice characterized by the in-plane $t_\parallel$ and out of plane $t_\perp$ boson tunneling amplitudes modifies the condensate formation and influences the thermodynamic properties of the superfluid and normal phases. To this end, we have determined the temperature dependence of the internal energy, entropy, heat capacity and compressibility. These quantities were studied as a function of the reduced temperature, i.e., T normalized to the characteristic energy scales of the system, namely the energy at criticality $k_\mathrm{B}T_\mathrm{c}$, as well as the bandwidth energy $t(\mathbf{0})\equiv 4t_\parallel +2t_\perp$. In this way, we were able to elucidate the role of the anisotropy in a coupled planes (chains) scenario and explicate the full 2D$\rightarrow$3D$\rightarrow$1D dimensional crossover. One of the questions that we seek to answer in the present work is how the lattice anisotropy influences the condensation temperature. We show, inter alia, that an extremely small value of the inter-plane coupling along the c-axis direction renders the bulk critical temperature of the 3D anisotropic system at the level of substantial fraction of its isotropic counterpart. This may be relevant for layered materials, as high-$T_\mathrm{c}$ superconductors, where the intricate interplay of in-plane and bulk properties dominate the physics in these systems. Furthermore, we pointed out that the theoretical outcome of the present paper may be of interest to a number of situations in solid-state physics, where the relevant quasi-particles are of bosonic nature. Besides, our findings might be useful to the trapped Bose condensates loaded in artificially engineered three-dimensional optical lattices. With this variety of applications, we can conclude that the Bose condensation in anisotropic lattice settings may provide interesting scenarios for the realization of quantum many-body physics with potential links to science and technology.

Appendix: Evaluation of 3D Lattice Integrals

The calculations of a number thermodynamic quantities throughout the body of the present paper rest on the ability to perform 3D wave vector integrals I over the Brillouin zone of the type

$$\begin{aligned} I=\int _{-\pi /a}^{\pi /a}\frac{\hbox {d}^3k}{(2 \pi /a)^3 }F(t(\mathbf{k})), \end{aligned}$$

(A.1)

where the momentum dependence enters a given function F(x) solely via the dispersion relation $t(\mathbf{k})$. Direct numerical evaluation of multi-dimensional integrals like I in Eq. (A.1) is highly impractical and (computational) time-consuming. Proper handling of the singularities, especially for strong hopping anisotropies where the system approaches 1D or 2D lattice limit, requires nonuniform meshes that are error prone. For example, near the BEC transition, integrands become singular but integrals remain finite. Proper handling of the singularities demands nonuniform 3D meshes that are error prone. Moreover, for strong hopping anisotropies the integrals diverge logarithmically. A practical and efficient solution of the problem can be achieved via the lattice density of states $\rho (\varepsilon )$ is defined as

$$\begin{aligned} \rho (\varepsilon )= \int _{-\pi /a}^{\pi /a}\frac{\hbox {d}^3k}{(2 \pi /a)^3 } \delta (\varepsilon -t(\mathbf{k})) . \end{aligned}$$

(A.2)

With the knowledge of $\rho (\varepsilon )$, the computational effort of performing momentum integrals is greatly reduced, as Eq. (A.1) simplifies to a single integral of the form

$$\begin{aligned} I=\int _{-\infty }^{+\infty }\hbox {d}\varepsilon \rho (\varepsilon )F(\varepsilon ). \end{aligned}$$

(A.3)

Further simplifications follows from the observation (as in our case) that $\rho (\varepsilon )$ has a finite support, which renders I a proper integral. Using the well-known Cauchy formula, which pertains to the Dirac distribution $\delta (x)$,

$$\begin{aligned} \lim _{\gamma \rightarrow 0+} \frac{1}{x \pm i\gamma } = \mathcal{P} \frac{1}{x} \mp i\pi \delta (x), \end{aligned}$$

(A.4)

where the symbol $\mathcal{P}$ refers to the principal value, we can relate the lattice Green’s function

$$\begin{aligned} G(\varepsilon )= \int _{-\pi /a}^{\pi /a}\frac{\hbox {d}^3k}{(2 \pi /a)^3 }\frac{1}{\varepsilon -t(\mathbf{k})} \end{aligned}$$

(A.5)

to the density of states via the imaginary part of the complex function $G(\varepsilon - i \gamma )$ via

$$\begin{aligned} \rho (\varepsilon )= \frac{1}{\pi }\mathfrak {I}\lim _{\gamma \rightarrow 0+} G(\varepsilon - i \gamma ). \end{aligned}$$

(A.6)

With the definition of the anisotropy parameters $\alpha , \alpha '$ given in Eqs. (5), (6) and (7), we write the Lattice Green’s function in Eq. (A.5) as,

$$\begin{aligned} G_\parallel (\varepsilon )= \int _{0}^{\pi /a}\frac{{2t_\parallel }^{-1}(\pi /a)^{-3}\hbox {d}k_x\hbox {d}k_y\hbox {d}k_z}{\frac{\varepsilon }{2t_\parallel }-\cos (k_xa)-\cos (k_ya)- \alpha \cos (k_z a)} \end{aligned}$$

(A.7)

Fig. 15

Evolution of the 3D lattice density of states $\rho (\varepsilon )$ as a function of the anisotropy parameter $t _\parallel /t_\perp$depicting the crossover from the isotropic 3D cubic structure $t_\perp = t _\parallel$ to the 1D lattice (the case isolated chains $t_\parallel =0$) (Color figure online)

Fig. 16

Evolution of the 3D lattice density of states $\rho (\varepsilon )$ as a function of the anisotropy parameter $t_\perp /t_\parallel$ from the isotropic 3D cubic structure $t_\perp = t _\parallel$ to the 2D lattice (the case of isolated planes $t_\perp =0$) (Color figure online)

For the case of coupled planes and

$$\begin{aligned} G_\perp (\varepsilon )=\int _{0}^{\pi /a}\frac{(2t_\perp )^{-1} (\pi /a)^{-3}\hbox {d}k_x\hbox {d}k_y\hbox {d}k_z}{\frac{\varepsilon }{2t_\perp }-\alpha ' [\cos (k_x a)+\cos (k_y a)]-\cos (k_z a)} \end{aligned}$$

(A.8)

in the case of coupled chains scenario. The essential point is now that the expressions for $G_\parallel (\varepsilon )$ and $G_\perp (\varepsilon )$ fit into the pattern of the so-called Watson integral formula of the form

$$\begin{aligned} G_W(\xi ,\alpha )=\int _{0}^{\pi /a}\frac{(\pi /a)^{-3}dk_xdk_ydk_z}{\xi -\cos (k_xa)-\cos (k_ya)- \alpha \cos (k_z a)} \end{aligned}$$

(A.9)

Indeed, with the following substitutions of the arguments of $G_W(\xi ,\alpha )$ we obtain

$$\begin{aligned} G_\parallel (\varepsilon )= & {} \frac{1}{2t_\parallel } G_W\left( \frac{\varepsilon }{2t_\parallel },\alpha \right) ,\nonumber \\ G_\perp (\varepsilon )= & {} \frac{1}{2t_\perp \alpha '} G_W\left( \frac{\varepsilon }{2t_\perp \alpha '},\frac{1}{\alpha '}\right) . \end{aligned}$$

(A.10)

Historically, the integrals of the type as in Eq. (A.9) with $\alpha =1$ emerged when studying cubic structures formed by real crystals and where first considered by Watson [48], so they are now universally known by his name. However, quite recently Delves and Joyce in an exceptional work of analysis [49, 50] succeeded in finding a closed form for the integral in Eq. (A.11) for a general value of the parameter $\alpha$. The exact analytic expression for $G_W(\xi ,\alpha )$ assumes a compact form

$$\begin{aligned} G_W({\xi },\alpha )= \frac{8}{\pi ^2\xi } \frac{\mathbf{K}(k_+)\mathbf{K}(k_-)}{\sqrt{1-(2-\alpha )^2w}+\sqrt{1-(2+\alpha )^2w}} , \end{aligned}$$

(A.11)

where $w=1/ \xi ^2$, while $\mathbf{K}(k)$ stands for the complete elliptic integral of the first kind [37] with modulus k. Furthermore, the arguments of the elliptic integrals in Eq. (A.11) are given by

$$\begin{aligned} 2k_\pm ^2= & {} 1-\left[ \sqrt{1-(2-\alpha )^2w}+\sqrt{1-(2+\alpha )^2w}\right] ^{-3} \nonumber \\&\times \left[ \sqrt{1+(2-\alpha ) \sqrt{w}}\sqrt{1-(2+\alpha ) \sqrt{w}}+ \sqrt{1-(2-\alpha ) \sqrt{w}}\sqrt{1+(2+\alpha ) \sqrt{w}}\right] \nonumber \\\times & {} \left\{ \pm 16w + \sqrt{1-\alpha ^2w}\left[ \sqrt{1+(2-\alpha ) \sqrt{w}}\sqrt{1+(2+\alpha ) \sqrt{w}}\right. \right. \nonumber \\&+\left. \left. \sqrt{1-(2-\alpha ) \sqrt{w}}\sqrt{1-(2+\alpha ) \sqrt{w}}\right] ^2 \right\} \end{aligned}$$

(A.12)

With the help of Eq. (A.11), the lattice density of states $\rho (\varepsilon )$ in Eq. (A.6) can be computed numerically. To this end, we employ the routines that implement special functions in the complex domain. The outcome was than used in calculations of various quantities throughout the paper and is illustrated in Figs. 15 and 16, where the dimensional crossover for the behavior of lattice density of states function is depicted.