Disordered Bose-Einstein condensate in hard walls trap

TL;DR

This paper investigates how quenched disorder affects a dilute Bose-Einstein condensate in a hard walls trap, analyzing the free energy, correlation functions, and density profile beyond mean-field approximation.

Contribution

It introduces a representation for the quenched free energy using integer moments of the partition function and explores effects of disorder and boundary conditions beyond mean-field.

Findings

Derived a series representation for quenched free energy.

Computed static two-point correlation functions in weak disorder.

Discussed the emergence of metastable and true ground states under strong disorder.

Abstract

We discuss the effects of quenched disorder in a dilute Bose-Einstein condensate confined in a hard walls trap. Starting from the disordered Gross-Pitaevskii functional, we obtain a representation for the quenched free energy as a series of integer moments of the partition function. Positive and negative disorder-dependent effective coupling constants appear in the integer moments. Going beyond the mean-field approximation, we compute the static two-point correlation functions at first-order in the positive effective coupling constants. We obtain the combined contributions of effects due to boundary conditions and disorder in this weakly disordered condensate. The ground state renormalized density profile of the condensate is presented. We also discuss the appearance of metastable and true ground states for strong disorder, when the effective coupling constants become negative.

Equations

$$[{\psi }^{†}(t,\mathbf{x}),{\psi }^{†}(t,{\mathbf{x}}^{\prime })]=[\psi (t,\mathbf{x}),\psi (t,{\mathbf{x}}^{\prime })]=0,$$

$$[\psi (t,\mathbf{x}),{\psi }^{†}(t,{\mathbf{x}}^{\prime })]={\delta }^d(\mathbf{x}-{\mathbf{x}}^{\prime }).$$

$$H({\psi }^{†},\psi )$$

$$N_0({\psi }^{†},\psi )=\int d^{d}x{\psi }^{†}(\mathbf{x})\psi (\mathbf{x}).$$

$$i\frac{\partial}{\partial\,t}\psi(t,\mathbf{x})=\bigl{[}\psi(t,\mathbf{x}),H(\psi^{\dagger},\psi)\bigr{]}=\left[-\frac{\Delta}{2m}+U(\mathbf{x})+\frac{g}{4}|\psi(t,\mathbf{x})|^{2}\right]\psi(t,\mathbf{x}),$$

$$\psi (t,\mathbf{x})=\left[\phi (\mathbf{x})+\vartheta (t,\mathbf{x})\right]e^{-i\mu t}$$

$$\left[-\frac{\Delta }{2m}+U(\mathbf{x})-\mu +\frac{g}{4}|\phi (\mathbf{x}){|}^2\right]\phi (\mathbf{x})=0.$$

$$Z(V)=\int[d\varphi][d\varphi^{*}]\exp\left[-\beta E(\varphi^{*},\varphi,V\bigr{)}\right],$$

$$E=\int d^{d}x\,\left[\varphi^{*}(\mathbf{x})\left(-\frac{\Delta}{2m}+V_{\text{trap}}(\mathbf{x})+V(\mathbf{x})-\mu\right)\varphi(\mathbf{x})+\,\frac{g}{4}\bigl{(}\varphi^{*}(\mathbf{x})\varphi(\mathbf{x})\bigr{)}^{2}\right].$$

$$\mathbb{E}[V(\mathbf{x})]=0\text{ and }\mathbb{E}[V(\mathbf{x})V({\mathbf{x}}^{\prime })]=\sigma {\delta }^d(\mathbf{x}-{\mathbf{x}}^{\prime }).$$

$$F=-\frac{1}{\beta }\int [dV]P(V)\ln Z(V),$$

$$\Phi (s)=\int [dV]P(V)\frac{1}{Z(V{)}^s},$$

$$F=\frac{1}{\beta }\left(\frac{d}{ds}\right)\Phi (s){|}_{s=0^{+}},\text{Re}(s)\ge 0,$$

$$\frac{1}{Z(V{)}^s}=\frac{1}{\Gamma (s)}{\int }_0^{\infty }dt{t}^{s-1}{e}^{-Z(V)t}.$$

$$F=\frac{1}{\beta}\Bigg{[}\sum_{k=1}^{\infty}\frac{(-1)^{k}a^{k}}{kk!}\mathbb{E}\left[Z^{k}\right]+\log{a}+\gamma-R(a)\Bigg{]},$$

$$\mathbb{E}\left[Z^k(V)\right]=\int [dV]P(V)Z^k(V),$$

$$R(a)=\int [dV]P(V){\int }_a^{\infty }\frac{dt}{t}{e}^{-Z(V)t}.$$

$$F=\frac{1}{\beta }\sum _{k=1}^{\infty }\frac{(-1{)}^k}{kk!}\mathbb{E}\left[Z^k\right],$$

$$\mathbb{E}\left[Z^k\right]=\int \prod _{i=1}^k\left[d{\phi }_i^{\ast (k)}\right]\left[d{\phi }_i^{(k)}\right]e^{-\beta E_{\text{eff}}^{(k)}\left({\phi }_i^{\ast (k)},{\phi }_i^{(k)}\right)},$$

$$E_{\text{eff}}^{(k)}(\varphi_{i}^{*(k)},\varphi_{i}^{(k)})=\int d^{d}x\sum_{i=1}^{k}\varphi_{i}^{*(k)}(\mathbf{x})\Biggl{[}-\frac{\Delta}{2m}+V_{\text{trap}}+\mu+\frac{1}{4}\sum_{j=1}^{k}g_{ij}|\varphi_{j}^{(k)}(\mathbf{x})|^{2}\Biggr{]}\varphi_{i}^{(k)}(\mathbf{x}),$$

$$g_{ij}=g{\delta }_{ij}-\beta \sigma .$$

$$\left(-\frac{△}{2m}+V_{\text{trap}}-\mu +\sum _{j=1}^k\frac{g_{ij}}{2}|{\phi }_j^{(k)}(\mathbf{x}){|}^2\right){\phi }_i^{(k)}(\mathbf{x})=0.$$

$$\left(-\frac{△}{2m}+V_{\text{trap}}-\mu +\frac{{\lambda }_k}{2}|{\phi }_i^{(k)}(\mathbf{x}){|}^2\right){\phi }_i^{(k)}(\mathbf{x})=0,$$

$$\{\begin{array}{l}{\phi }_i^{(k)}(\mathbf{x})={\phi }^{(k)}(\mathbf{x})\text{hfill}fork=1,2,...,N\\ {\phi }_i^{(k)}(\mathbf{x})=0fork>N.\end{array}$$

$$E_{\text{eff}}^{(k)}=k\int d^{d}x\left[\varphi^{*(k)}(\mathbf{x})\Bigl{(}-\frac{\Delta}{2m}+V_{\text{trap}}-\mu\Bigr{)}\varphi^{(k)}(\mathbf{x})+\frac{\lambda_{k}}{4}\varphi^{*(k)}(\mathbf{x})|\varphi^{(k)}(\mathbf{x})|^{2}\varphi^{(k)}(\mathbf{x})\right].$$

$$F=\frac{1}{\beta }\sum _{k=1}^{k_c}\frac{(-1{)}^k}{k!k}\mathbb{E}[Z^k]+\frac{1}{\beta }\sum _{k=k_c+1}^N\frac{(-1{)}^k}{k!k}\mathbb{E}[Z^k].$$

$$\mathbb{E}\left[Z^k\right]=\mathcal{N}\int \left[d\chi \right]\left[d{\chi }^{\ast }\right]e^{-\beta E_{\text{eff}}^{(k)}\left(\chi ,{\chi }^{\ast }\right)},$$

$$E_{\text{eff}}^{(k)}(\chi,\chi^{*})=k\int d^{d}x\,\left[\chi^{*}(\mathbf{x})\Bigl{(}-\frac{\Delta}{2m}+V_{\text{trap}}-\mu\Bigr{)}\chi(\mathbf{x})\frac{\lambda_{k}}{4}\chi^{*}(\mathbf{x})|\chi(\mathbf{x})|^{2}\chi(\mathbf{x})\right].$$

$$\chi(\mathbf{x})=\frac{1}{\sqrt{2}}\bigl{[}\phi_{1}(\mathbf{x})+i\phi_{2}(\mathbf{x})\bigr{]}$$

$$\chi^{*}(\mathbf{x})=\frac{1}{\sqrt{2}}\bigl{[}\phi_{1}(\mathbf{x})-i\phi_{2}(\mathbf{x})\bigr{]}.$$

$$E_{\text{eff}}^{(k)}({\varphi }_1,{\varphi }_2)$$

Full text

∗ E-mail address: racosta@if.ufrj.br

∗∗ E-mail address: carlos.zarro@if.ufrj.br

∘ E-mail address: gastao.krein@unesp.br

† E-mail address: asaldivar@cbpf.br

†† E-mail address: nfuxsvai@cbpf.br

1 Introduction

Effects due to boundary conditions in the dilute gases and condensates have been discussed by many authors [1, 2, 3]. In particular, the influence of hard walls type of boundary conditions on the statistical properties of an ideal Bose-Einstein condensate was discussed in Ref. [4]. In Ref. [5] the dynamics of an expanding condensate in a hard walls trap was investigated, finding that confinement generates structures in the condensate. Many other interesting situations can be created by coupling in a controlled manner disorder fields to a trapped Bose-Einstein condensate, allowing studies of the combined effect of boundary conditions and disorder.

Controlled disorder can be realized in ultracold Bose gases by using $e.g.$ laser speckle fields or incommensurate optical lattices [6, 7]. Under the influence of disorder, a bosonic system at low temperatures can contain many fragmented condensates [8, 9, 10, 11]. There are two main theoretical methods to investigate the effects of disorder in the Bose-Einstein condensate; methods based on Bogolibov theory [12, 13], and the replica method [14, 15, 16]. For low temperatures, small interaction strength and weak disorder, it is possible to use a perturbative approach. For strong disorder, the effects of the disorder must be treated non-perturbatively [17, 18, 19, 20]. To discuss the condensate in disordered environment, the replica method was used in Refs. [21, 22, 23].

Despite the successes in the application of the replica method in disordered systems, a mathematical rigorous derivation to support this procedure is still lacking [24, 25, 26, 27]. Intended to fill this gap, a mathematically rigorous technique for computing the average free energy of a system with quenched disorder in continuous field theory models was proposed in Refs. [28, 29, 30, 31]. The advantage of this approach is that contrary to the replica method, this technique neither involves derivatives of the $k$-th integer moments of the partition function $\mathbb{E}\left[Z^{k}\right]$ with respect to $k$, nor extension of these derivatives to non-integers values of $k$.

The purpose of this paper is to discuss the effects of disorder fields in a hard walls trapped dilute Bose-Einstein condensate. Here we are interested in the situation of a quenched disorder field coupled to the particle density of the system. Using mainly functional methods [32, 33, 34, 35], and with this new mathematical tool in hands, we investigate the density profile of the low temperature trapped dilute-gas Bose-Einstein condensate for weak disorder fields. Starting from a Gross-Pitaevskii energy functional, we employ the recently proposed method to evaluate the quenched free energy of disordered systems for arbitrary strengths of the disorder. A representation for the quenched free energy, described by a series of integer moments of the partition function is obtained. This representation is non-trivial and very rich in consequences. For some critical value of the strength of the disorder, an effective negative coupling constant appears in some integer moments. We call this situation the strong disorder limit of the model. Going beyond the mean-field approximation, we discuss the combined effect of the two-body interaction and a hard walls trap in the condensate for the situation of positive effective coupling constants, i.e., the weak disorder situation. From the static two-point correlation functions in the one-loop approximation, we obtain the combined contributions of effects due to boundary conditions and disorder on the density profile of the condensate. The renormalized density of the condensate which defines the ground state solution of the system is presented. To the best of our knowledge, these results are new. In addition, this is the first calculation of boundary effects in disordered Bose-Einstein condensates using field theoretical methods. For the cases in which effective negative coupling constants appear in some integer moments, there are metastable states and true ground states [36, 37, 38, 39, 40]. This is the scenario for first-order phase transitions induced by strong disorder fields.

This paper is organized as follows. In Sec. 2 we discuss the disordered Gross-Pitaevskii energy functional. In Sec. 3 we derive the disorder average free energy of the model for arbitrary strengths of the disorder. In Sec. 4 we discuss the effective Gross-Pitaevskii functional. In Sec. 5 we study the one-loop static two-point correlation functions in the presence of hard walls, in a generic integer moment of the partition function for the case where the effective coupling constant is positive. In Sec. 6 in this weak disorder limit, we obtain an expressions for the terms that contribute to the renormalized density of the condensate. In Sec. 7 the renormalized density of the condensate confined in a hard walls trap is presented. Conclusions are given in Sec. 8. In Appendix A we present some detailed calculations which were omitted through Sec. 6. To simplify the calculations we assume the units to be such that $\hbar=c=k_{B}=1$.

2 Gross-Pitaevskii energy functional

We consider a system of identical nonrelativistic spinless bosons in a $d$-dimensional spatial volume, where we assume that $d\geq 3$. For dilute systems, one can consider the bosons structureless and use the framework of nonrelativistic field theory. In the Heisenberg representation, the bosonic field operators $\psi^{\dagger}(t,\mathbf{x})$ and $\psi(t,\mathbf{x})$ satisfy the usual equal-time canonical commutation relations:

[TABLE]

For interacting Bose particles, the simplest situation is the one with a mutual two-body interaction represented by a instantaneous two-body potential $V_{0}(|\mathbf{x}-\mathbf{x}^{\prime}|)$, and an external potential $U(\mathbf{x})$. In this case, the Hamiltonian operator for the system is given by

[TABLE]

where $m$ is the mass of the bosonic particles and $\Delta$ is the Laplacian in $\mathbb{R}^{d}$. The operator for the total number of particles of the system is defined as

[TABLE]

In the absence of the external potential $U(\mathbf{x})$ the Hamiltonian operator is invariant under translations $(\mathbf{x}^{\prime}=\mathbf{x}+\mathbf{c})$ and SO(3) rotations in the three dimensional space. Also the Hamiltonian operator is invariant under a global rotation $\psi(\mathbf{x})\rightarrow\psi^{\prime}(\mathbf{x})=e^{i\alpha}\psi(\mathbf{x})$, for constant $\alpha$. The generator of the $U(1)$ gauge transformation is the particle number operator since $e^{-i\alpha\,N_{0}}\psi(\mathbf{x})\,e^{i\alpha\,N_{0}}=\psi(\mathbf{x})e^{i\alpha}$. The familiar Bogoliubov spectrum is the result of the spontaneous breaking of this continuous symmetry, a manifestation of the Goldstone theorem.

For dilute systems, one can replace the potential $V_{0}(|\mathbf{x}-\mathbf{x}^{\prime}|)$ by a contact interaction. In this case, from Eq. (3), one has that the Heisenberg equation of motion for the field operator $\psi(t,\mathbf{x})$ reads

[TABLE]

where $g$ is the strength of the contact interaction, related to the scattering length $a$ by $g={16\pi a}/{m}$. To study the stationary solutions of this Hamiltonian, it is usual to decompose the field operator into a condensate part, represented by a classical field $\varphi(\mathbf{x})$, the scalar order parameter, and a quantum fluctuating part, represented by ${\vartheta}(t,\mathbf{x})$ as:

[TABLE]

where $\mu$ is the chemical potential, which is determined by Eq. (4). When replacing this into Eq. (5) under the hypothesis that the effects of ${\vartheta}(t,\mathbf{x})$ are small, one obtains the time-independent Gross-Pitaevskii equation [41, 42, 43],

[TABLE]

As it is well known, this hypothesis is valid under certain conditions, namely, for a large number of particles and when the $s$-wave scattering length is much smaller than the mean separation between the particles. When treating the term ${\vartheta}(t,\mathbf{x})$ to first order in perturbation theory, one obtains the Bogoliubov spectrum of the quantum excitations of the condensate.

In many experimental setups the disorder couples to the particle density. In that case, one can write $U(\mathbf{x})=V_{\text{trap}}(\mathbf{x})+V(\mathbf{x})$, where $V_{\text{trap}}(\mathbf{x})$ is the potential representing the trap and $V(\mathbf{x})$ represents the disorder potential. In the following discussions, the shape of the trap potential is not important. We consider the Gross-Pitaesvikii framework, i.e. we neglect the fluctuating part of the field operator. For time-independent $V_{\text{trap}}$ and quenched disorder, we can employ the formalism of statistical field theory.

We remark that there is another approach to obtain a classical field theory. Suppose that we write the grand-canonical partition function as a path integral over complex field $\psi(\tau,\mathbf{x})$ and $\psi^{*}(\tau,\mathbf{x})$. Finite temperature effects are introduced where the energy integration is replaced by sum over the Matsubara frequencies. Therefore the complex field is a periodic function of the imaginary time with period $\beta$. After expand the field $\psi(\tau,\mathbf{x})$ and $\psi^{*}(\tau,\mathbf{x})$ in the path integral in terms of the Fourier components and restrict to the zero Matsubara component, the fields $\psi$ and $\psi^{*}$ can be considered as classical fields. The partition function is written in terms of a classical field theory where $\varphi_{0}(\mathbf{x})$ depends only on the spatial coordinates. This is called a classical field approximation.

Going back to our problem, the central object of interest is the disordered average free energy of the system, that is obtained from the logarithm of the disordered partition function:

[TABLE]

where $[d\varphi][d\varphi^{*}]$ is a measure in the space of all order-parameter configurations, with the disordered Gross-Pitaevskii energy functional $E(\varphi^{*},\varphi,V)$ given by

[TABLE]

To obtain the quenched free energy, $i.e.$ the average over the ensemble off all realizations of the quenched disorder, one needs to provide a probability distribution for the random potential $V(\mathbf{x})$, which we assumed as a delta-correlated Gaussian distribution with strenght $\sigma$:

[TABLE]

In the next section, we briefly review the method based on the distributional zeta function [28, 29, 30, 31] to compute the quenched free energy.

3 The quenched free energy

Given the probability measure, we are able to compute averages of physical quantities. Our main interest is to compute the quenched free energy:

[TABLE]

where $[d\,V]$ is a functional measure. Connected correlation functions are obtained by taking derivatives of $F$ with respect to an external complex source coupled to the complex order parameter in Eq. (9), as usual.

Inspired by situations in which one defines zeta functions in terms of countable collection of numbers, as for example prime numbers, non-trivial zeros of the Riemann zeta function, among others, we define the distributional zeta-function $\Phi(s)$ as

[TABLE]

for $s\in\mathbb{C}$, this function being well defined in an open connected subset of the complex plane. The quenched free energy is obtained from Eq. (12) as

[TABLE]

Next, use of the Euler’s integral representation for the gamma function $\Gamma(s)$ [44] allows us to write $1/Z(V)^{s}$ as

[TABLE]

Breaking this integral into two integrals, one from [math] to $a$ and another from $a$ to $\infty$ where $a$ is an arbitrary dimensionless real number, and expanding the exponential in powers of $t$ in the first integral, one can write the quenched free energy as:

[TABLE]

where $\mathbb{E}\left[Z^{k}(V)\right]$ is the $k$-th moment of the partition function

[TABLE]

$\gamma$ is the Euler constant and $R(a)$ is given by

[TABLE]

Notice that the disorder average free energy defined by Eq. (15) is $a$ independent.

The crucial question is the physical interpretation of Eq. (15). This equation indicates that the disordered system may be represented by an ensemble of subsystems, each of them being described by an integer moment of the partition function. In systems described by infinitely many degrees of freedom without disorder, the fundamental object is the generating functional of connected correlation functions, which are the moments of a probability measure. To describe disordered systems, the disordered-average generating functional of connected correlation functions is written as a series of the moments of the generating functional of correlation functions. As the result of such generalization, any renormalized quantity, as for example, the renormalized density of the condensate will be written also in a series of contributions coming from all moments of the partition function. This issue will be clarified further ahead in this paper.

Our next goal is to examine the fields in each integer moment of the partition function, for the concrete case of interest, the Gross-Pitaevskii energy functional of Eq.(9). This will be done in the next section.

4 The effective Gross-Pitaevskii energy functional

We have now at our disposal the Eq. (15), which, for large $a$, can be written as

[TABLE]

where $a^{k}$ factors have been absorbed in the functional measure of each integer moment of the partition function by redefining the fields which then amounts to redefining the parameters of the energy functional.

The object of interest is the $k$-th integer moment of the partition function, $\mathbb{E}\left[Z^{k}\right]$. From Eqs. (8), (9), (10) and Eq. (16) one can show that

[TABLE]

where

[TABLE]

and the generalized coupling constant $g_{ij}$ defined by

[TABLE]

For simplicity, we write $V_{\text{trap}}(\mathbf{x})=V_{\text{trap}}$. The saddle-point equations derived from the $k$-th moment read

[TABLE]

Using that $\varphi_{i}^{(k)}(\mathbf{x})=\varphi_{j}^{(k)}(\mathbf{x})$, the above equation becomes

[TABLE]

where we have defined the effective coupling constant $\lambda_{k}=g-\beta\sigma k$.

Consider the series in Eq. (18) where integer moment of the partition function is given by $\mathbb{E}\,[Z^{\,k}]$ — see also Eqs. (19) and (20). In our formalism, we have introduced auxiliary fields to describe the disordered condensate. As we discussed before, the consequence of such approach is a totally new picture for disordered systems where the ensemble of subsystems are not individually observed. Now we can imagine that there is a mapping between the number of particles of the condensate $n_{0}$, and the number of auxiliary fields. Since the number of bosonic particles in the condensate is finite, the number of fields also must be finite. We are claiming that the truncation of the series does not affect our ability to describe the disordered condensate. The consequence of such choice, is that the series given by Eq. (18) is truncated for some $N$. Hence, the structure of the fields in each integer moment of the partition function can be defined as

[TABLE]

This leads to the effective Gross-Pitaevskii energy functional:

[TABLE]

Due to the disorder, to describe the system we have to take into account contributions from fields governed by different effective coupling constants $\lambda_{k}$. It is important to stress that from these effective coupling constants it is possible to find the renormalized boson-boson interactions in the one-loop approximation using all the terms of the series for the average free energy. This issue will be investigated in the next sections where we will obtain the condensate density from the two-point correlations functions of the model.

For fixed $N$, let us define the weak disorder limit as the case when the effective coupling constants $\lambda_{k}$ in all $k$-th moments are positive. In the more general situation, there are two sets of moments: the integer moments where effective coupling constants are positive and a second set where the effective coupling constants become negative for some critical value for the strength of disorder. We called this situation the strong disorder limit of the model. Let us assume that this happens for some critical $k_{c}$. In the series representation for the quenched free energy there is a critical $k_{c}=\lfloor\frac{g}{\beta\sigma}\rfloor$. In this case the quenched free energy can be written as

[TABLE]

Let us investigate a generic integer moment of the series with $k<k_{c}$. Defining $\chi=\varphi^{(k)}$, we have

[TABLE]

where $\left[d\chi\right]\left[d\chi^{*}\right]$ is a functional measure and $\cal{N}$ is a normalization factor, and

[TABLE]

We are interested in calculating the ground states of this energy functional in the weak disorder limit. For accessing them, we work with the Cartesian representation for the complex field $\chi(\mathbf{x})$. We define the real fields $\phi_{1}(\mathbf{x})$ and $\phi_{2}(\mathbf{x})$ such that

[TABLE]

and

[TABLE]

The effective energy functional for these two real fields can be written as

[TABLE]

Let us define the functional potencial, $\Phi^{(k)}(\phi_{1},\phi_{2})$, which is bounded from bellow, as

[TABLE]

The $O(2)$ symmetry corresponds to an invariance under rotations in the real order parameters $(\phi_{1},\phi_{2})$ plane. The extremum of the potential functional in Eq. (32), for $\mu>0$ and the effective coupling constants $\lambda_{k}>0$, is at \bigl{(}\phi_{1}^{2}(\mathbf{x})+\phi_{2}^{2}(\mathbf{x})\bigr{)}=v^{2}, where $v^{2}=4\mu/\lambda_{k}$. The points on the circle with squared radius $v^{2}$ in the field $(\phi_{1},\phi_{2})$ plane are the minima of the potential. There is an infinite number of degenerate ground states. This is the standard situation of spontaneous symmetry breaking. Defining $\phi^{\prime}_{1}(\mathbf{x})=\phi_{1}(\mathbf{x})-v$ and $\phi^{\prime}_{2}(\mathbf{x})=\phi_{2}(\mathbf{x})$, the new energy functional $E_{\text{eff}}(\phi^{\prime}_{1},\phi_{2})$ is given by

[TABLE]

Two points are worth stressing. First, observables describing a stationary state in which the condensate structure is time independent can be expressed in terms of the static connected correlation functions. Second, in a standard field theory scenario, to go beyond the saddle-point approximation, one calculates the effective potential. This approach is suitable for systems without boundaries and a different technique needs to be considered in the present case of a trapped system, where translational invariance is broken.

For instance, in a field theory scenario it is difficult to study systematically finite size effects. The geometry of the system and boundary conditions determine the correlation functions of the model. Our aim is to investigate the combined effect of a hard walls trap and the interactions on the density profile of the disorder condensate. The density can be expressed in terms of the two-point static correlation functions and to go beyond the saddle-point approximation we consider the one-loop corrections to such a zeroth-order result. This approach is able to exhibit the ground state solution of the Bose-Einstein condensate for weak disorder and the presence of the trap. The ground state solution for the trapped Bose-Einstein condensate without disorder, for repulsive boson-boson interaction in $d=2$ and $d=3$ was discussed in Ref. [45].

A further step is taken if we remember that the number of particles in the condensate is given in terms of the density $\rho(\mathbf{x})$ as

[TABLE]

Defining a series representation for $\rho(\mathbf{x})$ through

[TABLE]

one can write

[TABLE]

Using Eqs. (26), (27) and (4), we get

[TABLE]

where $G_{1}^{(k)}(\mathbf{x},\mathbf{x})$ is the static two-point correlation function associated to the field $\phi^{\prime}_{1}(\mathbf{x})$ at coincident points, obtained from $k$-th integer moment of the partition function. Using dimensional analysis, the factor $m$ in the above equation must be introduced. In principle it is possible to use a resummation procedure [46, 47, 48], summing certain classes of diagrams from all orders of perturbation theory. Here we are using the one-loop approximation. In the following we will discuss the zeroth-order contribution and first-order corrections in the effective coupling constants to derive the two-point correlation function associated to the field $\phi^{\prime}_{1}(\mathbf{x})$.

5 The two-point static correlation functions in the presence of the trap

We consider a hard walls potential in a generic $d$-dimensional space with $d-1$ unbounded coordinates $\mathbf{r}=(x^{1},x^{2},\cdots,x^{d-1})$ and one bounded coordinate $0\leq z\leq L$. Therefore we have $\mathbf{x}=(\mathbf{r},z)$. The hard walls potential is defined by $V_{\text{trap}}=0$, for $0<z<L$ and $V_{\text{trap}}=\infty$ for $z\,\leq 0$ and $z\,\geq L$. We impose Dirichlet boundary conditions on the fields at the planes $z=0$ and $z=L$:

[TABLE]

and

[TABLE]

It is more convenient to use the notation $\varphi(\mathbf{x})$ for both fields $\phi^{\prime}_{1}(\mathbf{x})$ and $\phi_{2}(\mathbf{x})$. The fields that satisfy Dirichlet-Dirichlet boundary conditions on $z=0$ and $z=L$ can be expanded in Fourier series as:

[TABLE]

where $\mathbf{p}$ is the continuum momentum and $u_{n}(z)$ stands for the eigenfunctions of the operator $-\frac{d^{2}}{dz^{2}}$:

[TABLE]

satisfying the completeness and orthonormality relations. Due to the boundary condition we have $k_{n}=\frac{n\pi}{L}$, $n=1,2,3,...$.

From Eq. (40), the static two-point correlation functions associated with the $\phi^{\prime}_{1}(\mathbf{x})$ field, $G_{1}^{(0)}(\mathbf{r},z,\mathbf{r}^{\prime},z^{\prime},\mu)$ can be expressed as:

[TABLE]

The static two-point correlation functions associated with the field $\phi_{2}(\mathbf{x})$ is $G_{2}^{(0)}(\mathbf{r},z,\mathbf{r}^{\prime},z^{\prime})=G_{1}^{(0)}(\mathbf{r},z,\mathbf{r}^{\prime},z^{\prime},\mu)|_{\mu=0}$. Inspection of Eq. (4) reveals that there are three- and four-point vertices. Let us denote $G_{1}^{(1)}(\mathbf{r}_{1},z_{1},\mathbf{r}_{2},z_{2},\mu)$, the one-loop correction to $G_{1}^{(0)}(\mathbf{r},z,\mathbf{r}^{\prime},z^{\prime},\mu)$ of order $\lambda_{k}$ and $G_{1}^{(2)}(\mathbf{r}_{1},z_{1},\mathbf{r}_{2},z_{2},\mu)$, for corrections of order $(\lambda_{k}v)^{2}$, which is actually of order $\lambda_{k}$. The one-loop correction to the two-point correlation function that we defined as $G_{1}^{(1)}(\mathbf{r}_{1},z_{1},\mathbf{r}_{2},z_{2},\mu)$ is given by

[TABLE]

The second one, given by $G_{1}^{(2)}(\mathbf{r}_{1},z_{1},\mathbf{r}_{2},z_{2},\mu)$ has two contributions. We can write that $G_{1}^{(2)}(\mathbf{r}_{1},z_{1},\mathbf{r}_{2},z_{2},\mu)=G_{1a}^{(2)}(\mathbf{r}_{1},z_{1},\mathbf{r}_{2},z_{2},\mu)+G_{1b}^{(2)}(\mathbf{r}_{1},z_{1},\mathbf{r}_{2},z_{2},\mu)$. These one-loop contributions are given respectively by

[TABLE]

and

[TABLE]

In the next section we will obtain the correction to the density profile of the condensate due to the disorder and the presence of the hard walls trap in first-order approximation of the effective coupling constant. A detailed study of these two-point correlation functions of the fields $\phi^{\prime}_{1}(\mathbf{x})$ and $\phi_{2}(\mathbf{x})$ in the context of a Euclidean field theory with boundaries can be found in Refs. [49, 50, 51, 52]. Since we are taking the continuum limit of this field theory, all the ultraviolet divergences can be controlled using an analytic regularization procedure [53] with a renormalization scheme. It provides an adequate description of the effective long-distance physics in scalar field theory models.

6 The density profile of the condensate due to the disorder and hard walls trap

The aim of this section is to present, after a regularization and a renormalization procedure, the contribution to the density profile of the condensate coming from some integer moment of the partition function, produced by the disorder and the presence of a hard walls trap. As discussed before, we study this correction in first-order approximation of the effective coupling constant in the weak disorder limit. From a generic $k$-th moment of the partition function and using Eqs. (37), (42), (43) and (44), $\rho^{(k)}(z)$ is expressed as:

[TABLE]

where $\Delta\rho^{(0)}(z,\mu)=\Delta\rho_{1}^{(0)}(z,\mu)$, $\Delta\rho^{(1)}(z,\mu,\lambda_{k})$ and $\Delta\rho^{(2)}(z,\mu,\lambda_{k})$ are the contributions to the density profile from the two-point static correlation functions. To proceed we write $\Delta\rho^{(1)}(z,\mu,\lambda_{k})$ and $\Delta\rho^{(2)}(z,\mu,\lambda_{k})$ as

[TABLE]

and

[TABLE]

respectively. In the expressions above, the superscripts on the $\Delta\rho$ denote the power of the effective coupling constant in each contribution, and the subscripts denote that such contributions come from only the field $\phi^{\prime}_{1}(\mathbf{x})$ or both fields $\phi^{\prime}_{1}(\mathbf{x})$ and $\phi_{2}(\mathbf{x})$.

6.1 Zeroth-order result to the density profile

In order to get the first contribution to the density profile from the two-point static correlation functions, which is of zeroth-order in $\lambda_{k}$, let us define $\Delta\rho_{1}^{(0)}(z,\mu)$ as

[TABLE]

To compute this quantity, we introduce $T(z,\mu)\equiv m^{-1}G_{1}^{(0)}(\mathbf{r},z,\mathbf{r},z,\mu)$ as

[TABLE]

Using the fact that the mode solutions of Eq. (41) are $u_{n}(z)=\sqrt{\frac{2}{L}}\sin\left(\frac{n\pi\,z}{L}\right)$ and the following identity:

[TABLE]

the quantity $T(z,\mu)$ can be split into two contributions:

[TABLE]

the first one, $f_{1}(\mu)$, is independent of the distance to the hard wall and $f_{2}(z,\mu)$ depends on that distance. They are given by

[TABLE]

and

[TABLE]

In the infinite cut-off limit, it appears volume divergences in $f_{1}(\mu)$ and surface divergences in $f_{2}(z,\mu)$. Surface divergences are the ones that depend on the distance to the hard walls. We will discuss these divergences later on. Using the well known identity of dimensional regularization [54, 55, 56], i.e.,

[TABLE]

we can write $f_{1}(\mu)$ as

[TABLE]

The regularization of the above equation can be done using the definition of the modified Epstein zeta function. The modified Epstein zeta function $\xi(s,a)$ is a function of the complex variable $s=\sigma+it$, where $\sigma,t\,\,\in{\mathbb{R}}$, for $a^{2}>0$ defined by the absolutely convergent series as

[TABLE]

and in the whole complex plane $\mathbb{C}$ by analytic continuation, which is given by

[TABLE]

where $K_{\nu}(z)$ is the modified Bessel function of second kind and $\Gamma(z)$ is the gamma function. Using this analytic continuation, one obtains a size-independent polar contribution plus a size-dependent contribution. We find

[TABLE]

For odd-dimensional spaces, the gamma function is finite. For even $d$, one can use that for $n=0,1,2...$, $\Gamma(-n+\varepsilon)=\frac{(-1)^{n}}{n!}\left(\frac{1}{\varepsilon}+\psi(n+1)+\mathcal{(}\varepsilon)\right)$, where $\psi(z)$ is the digamma function. The introduction of bulk counterterms eliminate this divergence.

Let us discuss $f_{2}(z,\mu)$. We begin by an angular integration and $\int{}d\Omega_{d}=\frac{2\pi^{\frac{d}{2}}}{\Gamma(\frac{d}{2})}$ that leads to the following expression for $f_{2}(z,\mu)$, namely

[TABLE]

Using the change of variables $s=\sqrt{\mathbf{p}^{2}+4\mu m}$ in the above expression yields the following formula for $f_{2}(z,\mu)$:

[TABLE]

where $h(d)=\left[2(2{\sqrt{\pi}})^{d-1}\Gamma\left(\frac{d-1}{2}\right)\right]^{-1}$.

The function $f_{2}(z,\mu)$ has a regular contribution plus a singular one near the boundaries. To discuss this singular quantity, one can show that its singular structure is the same as in the case where $\mu=0$. A complete analysis of these divergences are found in Refs. [49, 50]. Therefore we study the case where $\mu=0$. In fact, we are particularly interested in examining the limits ($z\rightarrow{}0^{+}$ and $z\rightarrow{}L^{-}$) which obviously contain the information about the effect of a trap. In order to fulfill this goal we introduce two new variables $x=Ls$ and $q=zs$, in terms of which we can write $f_{2}(z,\mu)|_{\mu=0}$ as

[TABLE]

The first term of Eq. (62) gives a divergent contribution if we approach the boundary. The other term of Eq. (62) behaves as $\frac{1}{(L-z)^{d-2}}$. To see this let us investigate the behavior of the first integral of $f_{2}(z,\mu)|_{\mu=0}$ near the boundary, i.e., when $z\rightarrow{}0^{+}$ and $z\rightarrow{}L^{-}$. In order to do this, we make use of two formulas involving the definition for the gamma function, and also another well known integral representation for the product of the gamma function times the Hurwitz zeta function given by

[TABLE]

$Re(\beta)>0$, $Re(\alpha)>1$. The Hurwitz zeta function $\zeta(s,a)$ is a function of the complex variable $s=\sigma+it$, where $\sigma,t\,\,\in{\mathbb{R}}$, defined by the absolutely convergent series

[TABLE]

and in the whole complex plane $\mathbb{C}$ by analytic continuation. From the definition of the gamma function and using Eq. (63) in Eq. (62) we may write the following closed expression

[TABLE]

From this last expression and using the definition of the Hurwitz zeta function given by Eq. (64), it is evident that the regularized $f_{2}(z,\mu)|_{\mu=0}$ has two poles of order $(d-2)$, one at $z=0$ and another at $z=L$. Therefore $\int_{0}^{L}dz\,f_{2}(z,\mu)|_{\mu=0}$ diverges at $z=0$ and $z=L$. The standard way to circumvent this singular behavior is to renormalize such divergences by introducing surface counterterms at $z=0$ and $z=L$ [57]. For a careful discussion of the physical meaning of such divergences in quantum field theory can be found in Ref. [58].

Hence, the zeroth-order contribution to the density of the condensate is given by

[TABLE]

There is a uniform correction to the density due to the presence of the hard walls associated with the $\phi^{\prime}_{1}(\mathbf{x})$ field given by $m^{2}f_{1}(\mu)$ and a distortion given by $m^{2}f_{2}(z,\mu)$. The above calculation shows that the density profile of the condensate is highly sensitive to the presence of the hard walls potential.

6.2 First-order corrections to the density profile

We shall discuss the first contribution that depends on the combined effect of the trap, disorder and interactions. The second correction to the density profile is $\Delta\rho^{(1)}(z_{1},\mu,\lambda_{k})$. We are going to investigate the first term in Eq. (47), $\Delta\rho_{1}^{(1)}(z_{1},\mu,\lambda_{k})$, which is given by

[TABLE]

The other contribution to $\Delta\rho^{(1)}(z_{1},\mu,\lambda_{k})$ reads:

[TABLE]

The computation of these quantities are presented in details in Appendix A. Both contributions show the presence of spatial inhomogeneities.

6.3 Second type of first-order contributions to the density profile

Now, we are going to deal with the last group of first-order corrections to the density profile, which contain three-point vertices. Using the Eqs. (37), (44), (45) and (48), and defining $\Delta\rho_{1}^{(2)}(z_{1},\mu,\lambda_{k})=\Delta\rho_{1a}^{(2)}(z_{1},\mu,\lambda_{k})+\Delta\rho_{1b}^{(2)}(z_{1},\mu,\lambda_{k})$ one finds:

[TABLE]

and

[TABLE]

In the same way we define $\Delta\rho_{12}^{(2)}(z_{1},\mu,\lambda_{k})=\Delta\rho_{12a}^{(2)}(z_{1},\mu,\lambda_{k})+\Delta\rho_{12b}^{(2)}(z_{1},\mu,\lambda_{k})$ where $\Delta\rho_{12a}^{(2)}(z_{1},\mu,\lambda_{k})$ reads

[TABLE]

Finally the last term $\Delta\rho_{12b}^{(2)}(z_{1},\mu,\lambda_{k})$ reads

[TABLE]

Inserting Eq. (42) into the last group of first-order corrections to the density profile, one can find the contribution from $\Delta\rho_{1}^{(2)}$ and $\Delta\rho_{12}^{(2)}$. Also, the computation of these quantities can be obtained using the same method discussed in A. However, they will not be presented here as the sums cannot be decoupled.

7 The renormalized density of the condensate

From the above discussions, it is possible to write the renormalized density $\rho_{R}(z)$ in the case where all the terms of the series that represents the quenched free energy has effective positive coupling constants. We have

[TABLE]

where the constant $C_{N}$ is

[TABLE]

There are two contributions to the renormalized density. The first one depends on the presence of the hard wall. Disorder does not affect $\Delta\rho^{(0)}_{R}(z,\mu)$. This contribution is the same for the system without disorder. It is not difficult to find the contribution of $\Delta\rho^{(0)}_{R}(\mu)$ without the trap. The second one combines the effects of the trap, disorder and the two-body interaction, which were discussed before. These terms indicate the presence of inhomogeneities in the disordered condensate. One can estimate the effects on the condensate, due to disorder, as follows: the case without disorder is given by $\lambda_{k}=g=\text{constant}$, so $\Delta\rho^{(1)}_{R}(z,\mu,\lambda_{k})=\Delta\rho^{(1)}_{R}(z,\mu,g)$, which does not depend on $k$. The difference of the renormalized density with and without disorder is

[TABLE]

In the same way that we obtained from the disordered-average generating functional of connected correlation functions the renormalized density $\rho_{R}(z)$ in the first-order approximation of the effective coupling constant, it is possible to get the renormalized coupling constant between bosonic particles also in this first-order approximation. For attractive boson-boson interactions, the disorder produces an effect similar to that that of a Feshbach resonance, in that it changes the sign of the two-body interaction between the atoms [59]. Therefore, the scattering length becomes negative due to the effects of the disorder in this case.

For a weak disorder, one can increase the number of particles to obtain a negative effective coupling constant. Or one can also obtain a negative coupling constant increasing the strength of the disorder. The next problem which arises in the case where the disorder is strong. There is a critical $k$ such that $k_{c}=\left\lfloor\frac{g}{\beta\sigma}\right\rfloor$. In this case the quenched free energy has two contributions. The first one, for $k\leq k_{c}$ was analyzed. In the second one, the effective coupling constants are negative. In the conclusions we discuss briefly this situation.

8 Conclusions and Outlook

In this work we study the effects of disorder in a dilute Bose-Einstein condensate confined in a hard walls trap. We analyze the behavior of the system in equilibrium at temperatures below the critical condensation temperature for different values of the strength of the disorder. Physical quantities are calculated from the quenched free energy.

We start working at the mean-field level, where the random system is described by the disordered Gross-Pitaevskii energy functional. Then, the quenched free energy is written in a series of moments of the partition function of the system. This series representation describes the system as an ensemble of subsystems, each of them being characterized by a disorder-dependent effective coupling constant that can be positive or negative. In a generic moment of the partition function, the effective coupling constant becomes negative for some critical value of the strength of the disorder. We have implemented an analysis for a generic integer moment, in which the system is described by $k$ identical and non-interacting complex fields. Since all the fields are equal, one can study the contribution from a single field. The ground state of such system is discussed for different strengths of the disorder.

In the first-order approximation of the effective coupling constant, we discuss the combined effect of the hard walls and the two-body interaction. This can be done by studying the behavior of the two-point static correlation functions. Studying a generic integer moment with positive effective coupling constant, which corresponds to weak disorder, we present a general expression for its contribution to the density profile of the condensate. There are two contributions to the renormalized density. The first one depends on the presence of the hard walls. The second one combines effects of the trap, disorder and the two-body interaction.

Since it is possible to control the strength of the disorder thereby changing the magnitude of the effective coupling constants, it is always possible to produce a scenario for effective negative coupling constant in some integer moment of the partition function. In such a situation, it is necessary to introduce a three-particle interaction contribution [60] for defining a potential $\Phi(\phi_{1},\phi_{2})$ with a well-defined extremum. The new effective Gross-Pitaevskii energy functional, $E^{\prime\prime}_{\text{eff}}(\phi_{1},\phi_{2})$, can be written as

[TABLE]

where $E^{\prime}_{\text{eff}}(\phi_{1},\phi_{2})$ is given by Eq. (4). Note that we could have anticipated the need of such a term and included it from the very beginning in the Hamiltonian of Eq. (3), as one would do within an effective field theory framework [61]. Bose-Einstein condensate with a three-particle interaction term was also discussed in Ref. [62, 63]. The analysis of the existence of first and second order phase transitions for such functional with three-particle interaction contribution were discussed in Ref. [64]. This first-order phase transition in the condensate is induced by disorder fields.

A natural continuation of this paper is to use the series representation for the quenched free energy to discuss the behavior of the Bose-Einstein condensate with disorder in three different situations. First, still in the scenario of a stationary situation, discuss the strong disorder limit by taking into account the second contribution to the series defined by Eq. (26), and find the ground state of the system for strong-disorder fields. Second, going back to the weak disorder limit introduce dynamics in the problem, i.e., using Eq. (6), instead of the static Gross-Pitaevskii equation. In this second situation one is interested to show, for instance how the disorder modifies the sound velocity in the condensate. Finally the strong disorder limit with dynamics must be analyzed. These subjects are under investigation by the authors.

Acknowledgments

We would like to thank G. Menezes and for useful discussions. This work was partially suported by Conselho Nacional de Desenvolvimento Científico e Tecnológico - CNPq, 309982/2018-9 (C.A.D.Z.), 305894/2009-9 (G.K.), 303436/2015-8 (N.F.S.), INCT Física Nuclear e Aplicações, 464898/2014-5 (G.K) and Fundação de Amparo à Pesquisa do Estado de São Paulo - FAPESP, 2013/01907-0 (G.K).

Appendix A Computations of $\Delta\rho^{1}$

In this Appendix, we will present a detailed derivation of $\Delta\rho^{1}$. From Eqs. (43) and(48) one finds

[TABLE]

Inserting Eq. (42) into the above equation, one has

[TABLE]

The integral in $\mathbf{r}$ can be performed using the relation

[TABLE]

From the above equation the integral in $\mathbf{p}_{3}$ can be also evaluated. It remains to compute the integral in $z$. The product of sines can be written as:

[TABLE]

The $z$ integral can be computed using the completeness of $u_{n}(z)$:

[TABLE]

Thefore Eq. (78) can be written as

[TABLE]

where

[TABLE]

and

[TABLE]

The quantity $\Delta\rho_{12}^{(1)}$ is given by

[TABLE]

This quantity can be computed by the same steps done to calculate $\Delta\rho_{1}^{(1)}$ as the unique change consists in taking the limit $\mu\to 0$ in the $\mathbf{p}_{2}$ integral.

In order to renormalize $\Delta\rho^{(1)}_{\textbf{(a)}}$, we will follow the same steps as before to renormalize $\Delta\rho^{(0)}$. There are contributions similar to $f_{1}(\mu)$ and $f_{2}(\mu,z)$, however the complete analysis of that terms due to its complexity will be omitted here. The second term, $\Delta\rho^{(1)}_{\textbf{(b)}}$ cannot be computed directly as the sums are not decoupled. However, one can integrate in $z_{1}$ to circumvented this problem. Again, we will omit that computation. The same considerations apply to $\Delta\rho_{12\textbf{(a)}}^{(1)}$ and $\Delta\rho_{12\textbf{(b)}}^{(1)}$. These quantities can be written as

[TABLE]

where $\Delta\rho_{12\textbf{(a)}}^{(1)}$ is given by

[TABLE]

and $\Delta\rho_{12\textbf{(b)}}^{(1)}$ is

[TABLE]