Critical behavior at the spatial boundary of a trapped inhomogeneous Bose-Einstein condensate

TL;DR

This paper studies the critical behavior at the boundary of the Bose-Einstein condensate in inhomogeneous trapped quantum gases, revealing universal scaling laws through numerical simulations of the 3D Bose-Hubbard model.

Contribution

It demonstrates the critical phenomena at the BEC boundary in inhomogeneous traps and links it to the universality class of the homogeneous transition, supported by numerical evidence.

Findings

Critical behavior occurs at the BEC boundary.

Scaling behavior is governed by the universality class of the homogeneous transition.

Numerical simulations confirm the theoretical predictions.

Abstract

We investigate some aspects of the Bose-Einstein condensation (BEC) of quantum gases in the presence of inhomogeneous conditions. We consider three-dimensional (3D) quantum gases trapped by an external potential when the temperature is sufficiently low to show a BEC phase region around the center of the trap. If the trap is sufficiently large, different phases may coexist in different space regions, when moving from the center of the trap. We show that the quantum gas develops a peculiar critical behavior at the boundary of the BEC region, whose scaling behavior is controlled by the universality class of the homogenous BEC transition. We provide numerical evidence of this phenomenon, for lattice atomic gases modeled by the 3D Bose-Hubbard Hamiltonian.

Equations

$$H$$

$$G(\mathbf{x},\mathbf{y};T)\equiv ⟨b_{\mathbf{x}}^{†}{b}_{\mathbf{y}}⟩\equiv \frac{\mathrm{Tr}{b}_{\mathbf{x}}^{†}{b}_{\mathbf{y}}{e}^{-H/T}}{\mathrm{Tr}{e}^{-H/T}}$$

$$G({\mathbf{x}}_1,{\mathbf{x}}_2;T)\approx {\xi }^{-1-\eta }\mathcal{G}(|{\mathbf{x}}_1-{\mathbf{x}}_2|/\xi ),$$

$$G({\mathbf{x}}_1,{\mathbf{x}}_2;T=T_c)\sim \frac{1}{|{\mathbf{x}}_1-{\mathbf{x}}_2{|}^{1+\eta }},$$

$$V(r)=v^2{r}^2,r\equiv |\mathbf{x}|,$$

$$H_t=H+\sum _{\mathbf{x}}V(r)n_{\mathbf{x}}.$$

$${\mu }_e(\mathbf{x})\equiv \mu -V(r).$$

$$l_t\equiv \sqrt{t}/v.$$

$${\mu }_e(\mathbf{x})\equiv {\mu }_e(\mu ;X)=\mu -X^2,X=r/l_t,$$

$$G_{\mathrm{sw}}(\mathbf{x}-\mathbf{y})=⟨e^{-i\phi (\mathbf{x})}{e}^{i\phi (\mathbf{y})}⟩$$

$$H_{\mathrm{sw}}=\int d^{3}x[1+V(r)](\nabla \phi {)}^2,$$

$$G_0(\mathbf{x};l_t)\equiv G(0,\mathbf{x};l_t)\approx g_0(X),X=r/l_t.$$

$$T_c[{\mu }_e(\mu ;X)]\approx T<T_c(\mu ).$$

$$T_c[{\mu }_e(\mu ;X_b)]=T,$$

$$X_b\equiv r_b/l_t=\sqrt{\mu -\bar{\mu }},T_c(\bar{\mu })=T.$$

$$V(r)=V(r_b)+2X_{b}z/l_t+z^2/l_t^2,z\equiv r-r_b.$$

$$V_l(z)=\frac{z}{l_b},l_b=\frac{l_t}{2X_b}.$$

$${\xi }_b\sim l_b^{\theta },$$

$$Y\equiv z/l_b^{\theta }.$$

$$\theta =\frac{\nu }{1+\nu }=0.40181(3),$$

$$G_r(r_1,r_2;l_t,T)\approx {\xi }_b^{-1-\eta }{\mathcal{G}}_r(Y_1,Y_2),$$

$$A(\mathbf{x},\mathbf{y})\equiv ⟨n_{\mathbf{x}}{n}_{\mathbf{y}}⟩-⟨n_{\mathbf{x}}⟩⟨n_{\mathbf{y}}⟩.$$

$$A_r(r_1,r_2;l_t,T)\approx {\xi }_b^{-2y_n}{\mathcal{A}}_r(Y_1,Y_2),r_1\ne r_2,$$

$$\rho (\mathbf{x})\approx {\rho }_{\mathrm{lda}}(\mathbf{x})={\rho }_h[{\mu }_e(\mu ;X),T].$$

$$y_n\theta =\frac{3\nu -1}{\nu +1}\approx \mathrm{0.607.}$$

$$\rho (r_b)={\rho }_{\mathrm{lda}}+O(l_b^{-y_n\theta }),$$

$${\mathcal{G}}_r(0,Y)\sim |Y{|}^{-1-\eta },$$

$$G_r(r_b,r_b+z)\sim |z{|}^{-1-\eta }.$$

$$Z\equiv (L-r_b)/l_b^{\theta }.$$

$${\xi }_b(L{)}^2=\frac{{\sum }_{z=0}^{L-r_b}{z}^2{G}_r(r_b,r_b+z)}{{\sum }_{z=0}^{L-r_b}{G}_r(r_b,r_b+z)}$$

$${\xi }_b\equiv {\xi }_b(L\to \infty )\sim l_b^{\theta }.$$

Full text

I Introduction

Three-dimensional (3D) bosonic gases show the phenomenon of the Bose-Einstein condensation (BEC), which gives rise to a finite-temperature phase transition separating the high-temperature normal phase from the low-temperature superfluid BEC phase. The formation and properties of the BEC in quantum gases have been investigated by several experiments with cold atoms in harmonic traps, see, e.g., Refs. Andrews-etal-97 ; Stenger-etal-99 ; Hagley-etal-99 ; BHE-00 ; Cacciapuoti-etal-03 ; Hellweg-etal-03 ; Ritter-etal-07 ; BDZ-08 ; DGPS-99 ; PSW-01 ; Dettmer-etal-01 ; Hellweg-etal-02 . For generic 3D traps the coherence length turns out to be of the same size of the condensate (a substantial phase decoherence is only expected for a very elongated BEC, along the longer direction PSW-01 ; Dettmer-etal-01 ; Hellweg-etal-02 ; Mathey-etal-10 ; GCP-12 ; CDMV-16 ). The critical properties of 3D quantum gases at their BEC phase transition, i.e. when the condensate begins forming, have been investigated theoretically and experimentally, see, e.g., Refs. DRBOKE-07 ; DZZH-07 ; BB-09 ; CV-09 ; ZKKT-09 ; Trotzky-etal-10 ; HZ-10 ; PPS-10 ; NNCS-10 ; ZKKT-10 ; QSS-10 ; FCMCW-11 ; HM-11 ; Pollet-12 ; CR-12 ; CTV-13 ; CN-14 ; CNPV-15 ; CNPV-16 . In particular, the inhomogeneous conditions due to the spatially-dependent trapping potential give rise to a universal distortion of the critical behavior of the corresponding homogenous systems, i.e., without the external spatially-dependent trapping potential.

In this paper we consider 3D bosonic gases trapped by an external harmonic potential, when the temperature is sufficiently low to show a BEC-phase region around the center of the trap. Due to the trapping potential, the BEC region is generally spatially limited. Therefore, when moving from the center of the trap, the quantum gas passes from the BEC phase around the center of the trap (where space coherence is essentially described by spin waves) to a normal phase far from the center. We study the behavior of the correlation functions around the center of the trap and at the boundary of the BEC. We point out that the quantum many-body system develops a peculiar critical behavior at the boundary of the BEC region, with a nontrivial scaling behavior controlled by the universality class of the homogenous BEC transition. We provide some numerical evidence of this phenomenon in lattice atomic gases modeled by the Bose-Hubbard (BH) Hamiltonian FWGF-89 ; JBCGZ-98 .

The paper is organized as follows. Sec. II presents the 3D BH model, and the main features of its finite-temperature phase diagram. In Sec. III we discuss the effects of the inhomogeneous spatial conditions when the temperature is sufficiently low to show a BEC phase around the center of the trap; the scaling behaviors at the spatial boundary of the BEC are put forward. Sec. IV presents some numerical results for the 3D BH model, supporting the scaling ansatz at the spatial region between the BEC and normal-phase regions. Finally in Sec. V we draw our conclusions. The appendix reports the calculation of the critical exponent $\theta$ of the scaling theory describing the behavior of the one-particle correlation function at the boundary of the BEC.

II The 3D Bose-Hubbard model

The 3D BH model FWGF-89 is a realistic model of bosonic atoms in optical lattices JBCGZ-98 . Its Hamiltonian reads

[TABLE]

where $b_{\bf x}$ is a bosonic operator, $n_{\bf x}\equiv b_{\bf x}^{\dagger}b_{\bf x}$ is the particle density operator, the sums run over the bonds ${\langle{\bf x}{\bf y}\rangle}$ and the sites ${\bf x}$ of a cubic lattice, $a=1$ is the lattice spacing. The one-particle correlation function

[TABLE]

provides information on the phase coherence properties of the BEC.

The phase diagram of the 3D BH model and its critical behaviors have been much investigated, see e.g. Refs. FWGF-89 ; CPS-07 ; CR-12 ; CTV-13 ; CN-14 ; CNPV-15 . Their $T$-$\mu$ phase diagram presents a finite-temperature BEC transition line related to the formation of the condensate, i.e. the accumulation of a macroscopic number of atoms in a single quantum state. In Fig. 1 we sketch the phase diagram in the hard-core $U\to\infty$ limit. The BEC phase extends below the BEC transition line $T_{c}(\mu)$.

The condensate wave function provides the complex order parameter of the BEC transition. Thus its critical behavior belongs to the 3D XY universality class characterized by the spontaneous breaking of an Abelian U(1) symmetry PV-02 . In particular, around the transition point $T_{c}$, the one-particle correlation function of infinite homogenous systems behaves as

[TABLE]

where $\xi\sim\tau^{-\nu}$ is the diverging correlation length of the critical condensing modes, and $\tau\equiv T/T_{c}-1$ is the reduced temperature. The critical exponents $\nu$ and $\eta$ are those associated with the 3D XY universality class CHPV-06 , i.e., $\nu=0.6717(1)$ and $\eta=0.0381(2)$.

A common feature of the experiments with cold atoms BDZ-08 is the presence of an external potential $V({\bf x})$ coupled to the particle density, which traps the particles within a limited space region. In the experiments $V({\bf x})$ is usually effectively harmonic. We consider a rotationally-invariant harmonic potential (most results can be straightforwardly extended to more general cases)

[TABLE]

where $r$ is the distance from the center of the trap, which we locate at the origin ${\bf x}=0$. This trapping force is taken into account by adding a further term to the BH Hamiltonian (II), i.e.,

[TABLE]

Therefore, the external trapping potential coupled to the particle density turns out to be equivalent to an effective spatially-dependent chemical potential

[TABLE]

Far from the origin the potential $V(r)$ diverges, therefore the particle density $\rho({\bf x})\equiv\langle n_{\bf x}\rangle$ vanishes and the particles are trapped. We define the trap size by

[TABLE]

This definition naturally arises DSMS-96 ; BDZ-08 ; CV-10b when we consider the thermodynamic limit in a trap, which is generally defined by the limit $N,l_{t}\to\infty$ keeping $N/l_{t}^{3}$ fixed ($N$ is the number of particles). The value of the ratio $N/l_{t}^{3}$ is controlled by the chemical potential $\mu$: it remains constant when varying $l_{t}$ keeping $\mu$ fixed.

In the following, we set the hopping parameter $t=1$, so that all energies are expressed in units of $t$, and the Planck constant $\hslash=1$.

The inhomogeneity due to the trapping potential drastically changes, even qualitatively, the general features of the behavior of homogenous systems at a phase transition. For example, the correlation functions of the critical modes do not develop a diverging length scale in a finite trap. However, when the trap gets large the system shows a critical regime, although distorted by the presence of the trap. Around the transition point, i.e. when $T$ and $\mu$ are such that $T\approx T_{c}(\mu)$, the correlation functions around the center of the trap show power-law trap-size scaling behaviors with respect to the trap size $l_{t}$, controlled by the XY universality class of the critical behavior of the phase transition of the homogenous system CV-09 . The scaling behaviors of trapped BH models at the BEC transition are discussed in Refs. CTV-13 ; CN-14 , where also numerical analyses are reported.

III Scaling behaviors of the BEC

III.1 Coexisting phases in the presence of the trap

The inhomogeneous conditions due to the harmonic trapping potential, such as the BH model in Eq. (6), may give rise to separate regions where the quantum gas appears in the different BEC and normal phases.

For example, let us again consider the phase diagram sketched in Fig. 1, with the phase boundary $T_{c}(\mu)$. When $T<T_{c}(\mu)$ the quantum gas at the center of the trap is effectively in the BEC phase. Moving from the center, the effective chemical potential

[TABLE]

decreases. For sufficiently smooth changes, i.e. large trap size $l_{t}$, we can apply the local density approximation (LDA): the local observables around the spatial coordinate ${\bf x}$, and in particular the particle density, can be approximated by those of the homogeneous system at the same temperature and effective chemical potential $\mu_{e}({\bf x})$. When $r\equiv|{\bf x}|$ increases, $\mu_{e}=\mu-X^{2}$ decreases, until it reaches a region for which $T_{c}[\mu_{e}(\mu;X)]<T$, see Fig. 2. Thus, we pass from the BEC phase around the center of the trap to a region where the quantum gas is in the disordered normal phase. This scenario is sketched in Fig. 2. Analogous arguments apply to soft-core BH models with a finite on-site interaction parameter $U$.

III.2 Scaling behavior around the center of the trap

The space-coherence properties of the quantum gas within the BEC phase are essentially associated with spin-wave modes. For $T\ll T_{c}$, when the density $n_{0}$ of the condensate is much larger than the density of the noncondensed atoms, the particle-field operator of homogenous systems can be approximated PSW-01 by $b({\bf x})\approx\sqrt{n_{0}}e^{i\varphi({\bf x})}$. Then, the long-distance modes of the phase correlations are expected to be described by an effective spin-wave theory for the phase field $\varphi({\bf x})$, which is invariant under a global shift $\varphi({\bf x})\to\varphi({\bf x})+\alpha$. The simplest spin-wave Hamiltonian is given by $H_{\rm sw}=\int d^{3}x\;(\nabla\varphi)^{2}$. Actually, the region where the spin-wave theory effectively describes the long-distance phase correlations extends to the whole BEC phase PSW-01 ; CDMV-16 , i.e. for $T\lesssim T_{c}$, excluding only the relatively small critical region close to $T_{c}$. Therefore, the two-point spin-wave function

[TABLE]

is expected to describe the long-range phase-coherence properties of homogenous particle systems in the whole BEC phase. Computations of $G_{\rm sw}$ for various shapes and boundary conditions have been reported in Ref. CDMV-16 .

In the spin-wave limit, the inhomogeneity arising from the trapping potential $V(r)$ may be taken into account by considering the spin-wave Hamiltonian CV-11

[TABLE]

where the potential modulates the coefficient related to the superfluid density. Due to the Gaussian nature of the spin-wave theory, the scaling behavior of the correlation functions can be inferred by a straightforward dimensional analysis CV-11 . Therefore, the relevant spatial scaling variable is expected to be the ratio $X\equiv r/l_{t}$. For example, if we consider the one-particle correlation function, cf. Eq. (2), between the center of the trap and any point around the center, we expect the scaling behavior

[TABLE]

These considerations can be straightforwardly extended to correlations between generic points around the center, and to other observables.

III.3 Scaling behavior at the boundary of the BEC

We now argue that, for smooth trapping potentials, the quantum gas develops a peculiar critical behavior around the spatial surface separating the BEC and normal-phase regions, where

[TABLE]

For example, consider the hard-core BH lattice gas (6) for $\mu<0$ and $T<T_{c}(\mu)$, see Fig. 2. Since $T_{c}(\mu)$ decreases with decreasing $\mu$, a spherical surface exists at distance $r_{b}$ such that

[TABLE]

This surface separates the BEC region from the normal-fluid region. Expanding the rotationally-invariant potential $V(r)$ around $r_{b}$, we obtain

[TABLE]

Therefore, the quantum gas around the surface at distance $r_{b}$ effectively behaves as a BH model with chemical potential $\bar{\mu}=\mu_{e}(\mu;X_{b})$ such that $T=T_{c}(\bar{\mu})$, in the presence of an effectively linear potential $V_{l}(r)$ along the radial direction, given by

[TABLE]

Since $X_{b}(\mu,T,l_{t})>0$ is assumed finite and fixed, $l_{b}\sim l_{t}$.

We now derive some scaling ansatz meant to describe the asymptotic behavior of the observables around the critical surface, in the limit of large $l_{b}$. For this purpose, we note that the relevant spatial variable is expected to be the distance from the critical surface at $r=r_{b}$, i.e. $z=r-r_{b}$. This value should be compared with the length scale $\xi_{b}$ of the critical correlation functions around the surface $r=r_{b}$. Like general critical phenomena, see, e.g., Ref. PV-02 , the asymptotic scaling behavior of the length scale is expected to be characterized by a power law:

[TABLE]

where $\theta$ is an appropriate exponent. Therefore, the relevant scaling variable is expected to be given by the ratio

[TABLE]

According to Eq. (18), $\theta$ is the exponent that controls the relation between the length scale $\xi_{b}$ around the transition point $r_{b}$ and the length scale $l_{b}$ of the linear variation of the potential around $r_{b}$. The exponent $\theta$ can be determined by a scaling analysis of the perturbation associated with the external linear potential coupled to the particle density. Details are reported in the appendix. The exponent $\theta$ turns out to be related to the correlation-length exponent $\nu$ of the universality class of the critical behavior of the homogeneous BEC transition, i.e.,

[TABLE]

where $\nu=0.6717(1)$ is the correlation-length exponent of the 3D XY universality class.

We may apply the above considerations to the correlation function $G_{r}(r_{1},r_{2})$ along an arbitrary radial direction, i.e. along the points with spherical coordinates ${\bf x}=(r,\theta,\varphi)$ keeping the angle $\theta$, $\varphi$ fixed. We put forward the scaling behavior

[TABLE]

where $Y_{i}=z_{i}/l_{b}^{\theta}$ with $z_{i}\equiv r_{i}-r_{b}$, and $r_{1}\neq r_{2}$ (for equal points the scaling behavior differs, since $G_{r}(r,r)=\rho(r)$, see below). The prefactor is analogous to that reported in Eq. (3) for the critical behavior of homogenous systems. The difference among the different radial directions, due to the lattice structure of the BH model, is expected to be suppressed in the large-$l_{t}$ limit. Analogous scaling relations can be straightforwardly derived for correlations between generic points around the boundary between the two phases.

A scaling ansatz can be analogously derived for the connected density-density correlation function

[TABLE]

Along the radial direction,

[TABLE]

where $y_{n}=3-1/\nu$ is the renormalization-group dimension of the density operator at the XY fixed point CTV-13 .

The above scaling results can be straightforwardly extended to different geometries of the trap and generic space dependences of the external potential, when BEC and normal phases appear simultaneously due to the inhomogeneous conditions. Indeed, the linear approximation at the spatial surface between the two phases is quite general. It arises from general spatial dependences of the potentials, because a linear term generally appears in the expansion of the external potential around the critical surface. Using scaling arguments one can show that the quadratic (and in general any higher-order) term of the expansion of the potential, cf. Eq. (16), gives rise to suppressed contributions with respect to the asymptotic behavior (21). Indeed, close to the BEC boundary, when considering the scaling behavior at fixed ratio $Y\equiv z/l_{b}^{\theta}$ as in Eq. (21), the asymptotic behavior is determined by the smallest $O(z)$ power law, while the quadratic $O(z^{2})$ term gives only rise to $O(l_{b}^{-1+\theta})$ corrections to the leading behavior (21). Therefore, the precise form of the trapping potential turns out to be irrelevant for our considerations on the critical behavior around the spatial separation of the different phases.

Actually we expect that there are other power-law corrections to the asymptotic behavior, i.e., those related to the critical behavior of the corresponding homogenous system, and in particular with the scaling-correction exponent $\omega$ of the 3D XY universality class CHPV-06 ; GZ-98 ; PV-02 , i.e. $\omega=0.785(20)$. They also induce corrections in the presence of a linear external potential. They should be $O(\xi_{b}^{-\omega})$, thus $O(l_{b}^{-\omega\theta})$. Since $\omega\theta=0.315(8)$ is smaller than $1-\theta=0.59819(3)$, the scaling corrections arising from the leading irrelevant perturbation of the XY universality class provide the leading asymptotic corrections to the problem at hand.

The particle density $\rho({\bf x})\equiv\langle b_{\bf x}^{\dagger}b_{\bf x}\rangle$ is expected to vary smoothly across the transition surface. Its $l_{t}\to\infty$ asymptotic behavior approaches a finite value, even at the transition surface. This should be given by its LDA, i.e., by the particle density $\rho_{h}$ of the homogenous system at the corresponding values of the chemical potential and temperature, i.e.,

[TABLE]

Thus the $l_{t}\to\infty$ limit of the particle density is expected to be a function of the ratio $X\equiv r/l_{t}$. Analogously to the particle density of homogenous systems at criticality CTV-13 , the critical behavior at the boundary of the BEC region is expected to give rise to a nonanalytic subleading $O(\xi_{b}^{-y_{n}})=O(l_{b}^{-y_{n}\theta})$ contribution, where

[TABLE]

An analogous behavior has been also put forward, and numerically checked, in the case of the 2D Ising model in the presence of a temperature gradient BDV-14 .

We also note that more complicated scenarios may appear in the presence of inhomogeneous trap. For example, assuming the phase diagram of Fig. 1 for the hard-core $U\to\infty$ limit, when $T<T_{c}(\mu=0)$, if $\mu>0$ is larger than the critical value $\mu_{c}(T)>0$ (inverse function of $T_{c}(\mu)$ for $\mu>0$), then the gas at the center of the trap is in the disordered normal phase. At the distance such that $\mu_{e}(r)=\mu_{c}(T)$ it passes to a superfluid phase until $\mu_{e}(r)=-\mu_{c}(T)$, where it passes again to a normal phase. In this case we have two transition surfaces. The critical behaviors at such surfaces, in the regime of smooth external potential, is expected to show scaling behaviors analogous to Eq. (21). The exponent $\theta$ is the same of that given in Eq. (20), because at these transition surfaces the spatial variation of the external potential is effectively linear.

Finally, we mention that quantum zero-temperature transitions between spatially separated phases have been discussed in Refs. CV-10b ; CTV-12 for low-dimensional quantum gases, between the superfluid and vacuum phases. They present similar scaling behaviors, although the quantum nature of their critical modes makes these phenomena substantially different from those at finite temperature, as those discussed in this paper.

IV Numerical results

In order to check the scaling behaviors put forward in the previous section, we consider the 3D BH model in the hard-core $U\to\infty$ limit, whose phase diagram is sketched in Fig. 1. We present numerical results for the hard-core BH model at $\mu=-3$ and $T=1.482$ in the presence of a harmonic trap, cf. Eq. (6). The chosen value of the temperature corresponds to the critical temperature for $\mu=-4$, indeed CTV-13 $T_{c}(\mu=-4)=1.4820(2)$. This choice allows us to determine the critical surface with great accuracy, see below, avoiding the numerical uncertainty related to its location. Of course, whenever $T<T_{c}(\mu)$, the main features of the BEC scaling behaviors of trapped BH models are expected to be largely independent of this particular choice of parameters.

The trap is centered in the middle of a cubic $M^{3}$ lattice with open boundary conditions. We consider odd sizes $M=2L+1$, so that the lattice-site coordinates are $x_{i}=-L,...,L$, and the trap center is at ${\bf x}={\bf 0}$. Numerical results are obtained by quantum Monte Carlo (QMC) simulations using the directed operator-loop algorithm SK-91 ; SS-02 ; DT-01 (see Ref. CTV-12 for further details). We report results for several values of the trap size $l_{t}$ (up to $l_{t}=12$), for sufficiently large lattice size $L$ (up to $L=22$) to obtain the infinite-$L$ limit of the quantities we are interested in, within the typical statistical errors of our computations (this is easily checked by comparing data for increasing values of $L$, see also below).

According to the LDA discussed in the previous section, for sufficiently large trap size $l_{t}$, the particle density $\rho({\bf x})\equiv\langle b_{\bf x}^{\dagger}b_{\bf x}\rangle$ is expected to be a function of the local effective chemical potential $\mu_{e}({\bf x})=\mu-X^{2}$ where $X\equiv r/l_{t}$, thus a function of $X$. This is confirmed by the data for the density profile shown in Fig. 3. They clearly approach a function of $X$ with increasing $l_{t}$. We expect that their convergence is generally characterized by $O(l_{t}^{-1})$ corrections (related to the effective linear variation of the potential around points at $r>0$), except at the boundary of the BEC region, where it is expected to be slower (due to the critical modes, see below), and at center of the trap where the convergence is expected to be $O(l_{t}^{-2})$ (here the spatial dependence of the potential is quadratic). The numerical data turn out to be consistent. For example, for $X=1/2$, where the effective chemical potential has a linear slope, the data nicely fit $\rho_{1/2}+cl_{t}^{-1}$ (we obtain $\rho_{1/2}=0.2411(3)$ and $\chi^{2}/{\rm d.o.f}\approx 0.5$ using the data for $l_{t}\geq 6$). At the center of the trap, the data favor the behavior $\rho_{0}+cl_{t}^{-2}$ (we obtain $\rho_{0}=0.2592(4)$ and $\chi^{2}/{\rm d.o.f}\approx 0.1$ from the data for $l_{t}\geq 6)$.

The particle density $\rho({\bf x})$ shows a smooth behavior also across the transition surface at $r_{b}=l_{t}$, see Fig. 3. Indeed, like the energy density at a continuous transition PV-02 , its large-$l_{b}$ limit is expected to be a nonuniversal constant given by its LDA, while critical modes give rise to subleading $O(\xi_{b}^{-y_{n}})$, thus $O(l_{b}^{-y_{n}\theta})$, contributions, cf. Eq. (25). Therefore, we expect

[TABLE]

where $\rho_{\rm lda}$ is the value of the particle density of the homogenous system at $\mu=-4$ and $T_{c}(\mu=-4)=1.482$, i.e. CTV-13 $\rho_{\rm lda}=\rho_{c}(\mu=-4)=0.16187(1)$. The data for the particle density $\rho(r_{b})$ are consistent, see Fig. 3; they approach $\rho_{\rm lda}$, indeed $\rho(r_{b})=0.1600(2),\,0.1608(2),\,0.1613(2)$ for $l_{t}=8,\,10,\,12$ respectively.

Figure 4 shows data for the one-particle correlation function $G_{0}(r)$ between the center and a generic point, cf. Eq. (12). Their behavior is substantially consistent with the large-$l_{t}$ scaling behavior (12) inferred by spin-wave arguments, i.e., $G_{0}(r)\approx g_{0}(X)$ in the large trap-size limit (with a $X$-dependent convergence rate).

We now consider the behavior at the spatial boundary of the BEC, where the quantum gas is expected to develop a critical behavior. According to the scenario put forward in Sec. III.1, the spatial boundary between the BEC and normal phase regions is located at the surface $r=r_{b}$ corresponding to the effective chemical potential $\mu_{e}(\mu;X_{b})$ were $X_{b}=r_{b}/l_{t}$, such that $T=T_{c}[\mu_{e}(\mu;X_{b})]$. For our choice of parameters, $\mu_{e}(\mu;X_{b})=-4$, thus $r_{b}=l_{t}$. Therefore, the effective potential (17) around the critical surface is characterized by the length scale $l_{b}=l_{t}/2$.

In order to check the scaling ansatz (21), we compute the one-particle correlation function $G_{r}(r_{b},r_{b}+z)$ along one of the main axis of the lattice, around the critical point $r=r_{b}=l_{t}$. The data are shown in Fig. 5. They nicely support the scaling ansatz (21). Indeed the data for the product $l_{b}^{(1+\eta)\theta}G_{r}(r_{b},r_{b}+z)$ appears to approach a scaling function ${\cal G}_{r}(0,Y)$ of $Y\equiv z/l_{b}^{\theta}$ with increasing $l_{b}$ [except at $z=0$ where such a scaling is not expected, because $G_{r}(r_{b},r_{b})=\rho(r_{b})$, see Eq. (26)]. Note that the data for small values of $|Y|$ are substantially consistent with the asymptotic power-law behavior

[TABLE]

matching the small-distance critical behavior (4) of homogenous systems. Therefore, in the neighborhood of the critical surface, the one-particle correlation function shows the typical singular behavior of critical phenomena, cf. Eq. (4), i.e. when $|z|\ll l_{b}^{\theta}$

[TABLE]

In Fig. 5 the corrections to the asymptotic scaling function ${\cal G}_{r}(0,Y)$ appear larger for $Y<0$, corresponding to the BEC phase region, while the scaling approach turns out to be quite rapid for $Y>0$ corresponding to the normal phase region. As argued in the previous section, the leading corrections to the asymptotic behavior are expected to be $O(l_{b}^{-\omega\theta})$ due to the leading irrelevant perturbation of the 3D XY universality class. The effects of the higher-order terms (beyond the linear one) in the expansion of the potential around the critical surface, cf. Eq. (16), are expected to be more suppressed, i.e., they are $O(l_{b}^{-1+\theta})$.

Analogous results are expected along the other directions, indeed the effects of the breaking of the rotational symmetry due to the lattice structure are expected to be $O(l_{t}^{-1})$ suppressed in the large-$l_{t}$ limit.

We finally discuss the finite-size effects arising from the finite size $M=2L+1$ of the lattice. The above results have been obtained for sufficiently large $L$ to make finite-size effects negligible. This condition can be made more precise. Since the length scale of the critical modes around the critical sphere is $\xi_{b}\sim l_{b}^{\theta}$, we expect that finite-size effects get suppressed when $(L-r_{b})/\xi_{b}\gg 1$, thus $(L-r_{b})/l_{b}^{\theta}\gg 1$, where $L-r_{b}$ is the distance of the critical sphere from the boundary of the lattice. Actually, we expect that the interplay between $L$ and $l_{b}$ leads to finite-size scaling (FSS) behaviors analogous to FSS in homogenous models Cardy-88 . In the case at hand the relevant scaling variable is expected to be given by the ratio

[TABLE]

In order to show this further feature of the critical behavior at the boundary of the superfluid region, we consider the correlation length $\xi_{b}$ of the critical modes around the critical surface, defined as

[TABLE]

where the sum is meant along one of the main directions of the lattice. Eq. (21) implies that

[TABLE]

This is confirmed by the data of $\xi_{b}$ shown in the inset of Fig. 6, obtained for sufficiently large lattice sizes $L>r_{b}=l_{t}$ to provide a good approximation of $\xi_{b}(L\to\infty)$ (within an accuracy of one per cent, as checked by increasing $L$ at fixed $l_{t}$). Indeed, the data from $l_{b}=2$ to $l_{b}=6$ nicely fit the ansatz $\xi_{b}=c\,l_{b}^{\theta}$ with $\theta$ given by Eq. (20), i.e. $\theta=0.40181(3)$, see the inset of Fig. 6. In particular, considering $\theta$ as a free parameter, we obtain a consistent value with an accuracy of a few per cent, i.e. $\theta=0.407(3)$ using all available data, with an acceptable $\chi^{2}/{\rm d.o.f.}\approx 0.1$.

Then, we compute the ratio

[TABLE]

for various values of $l_{t}$ and $L$. The FSS hypothesis implies that $s_{\xi}$ must approach a scaling function of $Z$, cf. Eq. (29). This is nicely supported by the data reported in Fig. 6. This FSS analysis shows that the condition $Z\gtrsim 4$ is sufficient to provide the $L\to\infty$ limit of $\xi_{b}$ within approximately one per cent.

We finally note that the problem that we consider, i.e. the critical behavior at the boundary of the BEC region, differs from the phenomena of boundary critical behaviors at the boundaries of statistical systems Binder-83 ; Diehl-86 . In the latter case the inhomogeneous behavior arises from the presence of the boundaries. In the case that we study the inhomogeneous critical behavior is controlled by the external effectively linear potential, giving rise to a critical region of size $\xi\sim l_{b}^{\theta}$ separating the BEC and normal-phase regions.

V Conclusions

We investigate some aspects of the BEC of quantum gases in the presence of inhomogeneous conditions, such as cold-atom systems in harmonic traps, which are typically realized in experiments, see, e.g., Refs. Andrews-etal-97 ; Stenger-etal-99 ; Hagley-etal-99 ; BHE-00 ; Cacciapuoti-etal-03 ; Hellweg-etal-03 ; Ritter-etal-07 ; BDZ-08 ; DGPS-99 ; PSW-01 ; Dettmer-etal-01 ; Hellweg-etal-02 . We consider 3D bosonic gases trapped by an external potential, when the temperature is sufficiently low to show a superfluid phase region around the center of the trap. In particular, we consider 3D BH models in the presence of a harmonic trap, cf. Eq. (6), which model realistic gases of bosonic atoms in optical lattices JBCGZ-98 .

We point out that, if the trap is sufficiently large, different phases may coexist in different space regions, when moving from the center of the trap. For example, we may pass from the superfluid BEC phase around the center of the trap (where space coherence is essentially described by spin waves) to a normal phase far from the center. Between the superfluid and normal phase regions the quantum gas experiences a spatial phase transition. The system develops a peculiar critical behavior at the surface separating the different phases, characterized by the presence of an external potential with an effectively linear space dependence. Using scaling arguments, we put forward the asymptotic behavior in the limit of smooth external potential, i.e., large length scale $l_{b}$ associated with the effective linear variation of the potential at the critical surface, cf. Eq. (17). We argue that this peculiar scaling behavior is controlled by the universality class of the homogenous BEC transition, i.e., the 3D XY universality class whose critical exponents are known with great accuracy. The length scale $\xi_{b}$ of the critical modes around the critical surface diverges as $\xi_{b}\sim l_{b}^{\theta}$ with $\theta=\nu/(1+\nu)$, where CHPV-06 $\nu=0.6717(1)$ is the correlation-length critical exponent of the 3D XY universality class.

We provide numerical evidence of this phenomenon for the 3D BH model in the hard-core $U\to\infty$ limit (whose phase diagram is sketched in Fig. 1), in the presence of a rotationally-invariant harmonic trap. The numerical results, obtained by QMC simulations, nicely support our scaling predictions, in particular those at the boundary of the BEC. An analogous scenario is expected for soft-core BH models with finite on-site interaction parameters $U$.

The scaling theory at the boundary of the BEC can be straightforwardly extended to different geometries of the trap, and/or generic space dependences of the external potential. This is essentially due to the fact that the linear approximation at the spatial surface between the two phases is quite general. Moreover, high-order terms are expected to be irrelevant: they give only rise to $O(l_{b}^{-1+\theta})$ suppressed contribution in the $l_{b}\to\infty$ limit.

We remark that the lattice structure of the BH model does not play any particular role in our scaling arguments concerning the spatial phase transition at the boundary of the BEC, and the derivation of the corresponding power laws. Indeed the microscopic details of the model are irrelevant in the asymptotic $l_{b}\to\infty$ limit. Therefore, our scaling predictions at the boundary of the BEC apply to a wide class of models, i.e. to any 3D inhomogeneous interacting bosonic systems at the spatial boundary of their BEC phase.

The experimental evidence of the critical behavior at the boundary of the BEC in trapped quantum gases requires the measurement of the one-particle correlation functions, or higher-order correlations. These quantities are not easily accessible. However several examples of such space-coherence measurements have been reported in the literature, see, e.g., Refs. Hagley-etal-99 ; BHE-00 ; DRBOKE-07 ; NGSH-15 . On the other hand, the more accessible particle density, which can be measured by in situ density image techniques GZHC-09 ; BGPFG-09 , can hardly provide evidence of this critical phenomenon, since the critical modes give only rise to subleading contributions to the particle density, cf. Eq. (26).

Appendix A Derivation of the exponent $\theta$

The exponent $\theta$ can be inferred by a scaling analysis of the perturbation $P_{V}$ representing the external linear potential $V_{l}(r)=ur$ coupled to the particle density. We follow the field-theoretical approach of Refs. CV-09 ; CV-10 , that is we consider the 3D $\Phi^{4}$ quantum field theory which represents the 3D XY universality class, see e.g. Ref. ZJ-book ,

[TABLE]

where $\psi$ is the complex field associated with the order parameter, and $r,u$ are coupling constants. Since the particle density corresponds to the energy operator $|\psi|^{2}$, we can write the perturbation $P_{V}$ as

[TABLE]

Introducing the renormalization-group dimension $y_{u}$ of the constant $u$ of the linear potential, we derive the scaling relation $y_{u}-1+y_{n}=3$, where $y_{n}=3-1/\nu$ is the RG dimension of the density/energy operator $|\psi|^{2}$ (we recall that $\nu$ is the length-scale critical exponent of the 3D XY universality class). We eventually obtain $\theta\equiv 1/y_{u}$, and therefore Eq. (20).