Renormalized contact interaction in degenerate unitary Bose gases

TL;DR

This paper develops a renormalized contact interaction for degenerate unitary Bose gases, enabling the study of their universal properties at large scattering lengths using mean-field and hyperspherical methods.

Contribution

It introduces a new renormalization of the contact interaction based on exact two-particle solutions, extending applicability to unitarity in Bose gases.

Findings

Universal energy relation: $E/N=12.67rac{\hbar^2}{2m}\langle n^{2/3} angle$

Universal contact relation: $C_{2}/N=11.8\langle n^{1/3} angle$

Applicable to large scattering lengths at unitarity

Abstract

We renormalize the two-body contact interaction based on the exact solution of two interacting particles in a harmonic trap. This renormalization extends the validity of the contact interaction to large scattering lengths. We apply this renormalized interaction to a degenerate unitary Bose gas to study its stationary properties and elementary excitations using the mean-field theory and the hyperspherical method. Since the scattering length is no longer a relevant length scale at unitarity, universal properties are obtained that depend only on the average particle density. Our treatment shows that the universal relations for the total energy and for the two-body contact are $E/N=12.67\hbar^2\langle n^{2/3}\rangle/2m$ and $C_{2}/N=11.8\langle n^{1/3}\rangle$ respectively.

Equations

$$V(\text{r})=\frac{4\pi {\hslash }^{2}a}{m}\delta (\text{r})\frac{\partial }{\partial r}r,$$

$$H_{2b}=-\frac{{\hslash }^2}{2m}({\nabla }_{{\text{r}}_1}^2+{\nabla }_{{\text{r}}_2}^2)+\frac{1}{2}m{\omega }_{ho}^2({\text{r}}_1^2+{\text{r}}_2^2)+V({\text{r}}_{12}),$$

$$\sqrt{2}\frac{\Gamma \left(-E_{rel}(a)/2\hslash {\omega }_{ho}+3/4\right)}{\Gamma \left(-E_{rel}(a)/2\hslash {\omega }_{ho}+1/4\right)}=\frac{l_{ho}}{a},$$

$$\tilde{V}(\text{r})=\frac{4\pi {\hslash }^2{a}_{eff}}{m}\delta (\text{r})=\frac{4\pi {\hslash }^2\zeta (k_{F}a)}{m{k}_F}\delta (\text{r}),$$

$$\xi \{\psi \}$$

$$E_{HF}(a_{eff})=\xi \{\psi \}{|}_{\delta \xi /\delta \psi =0}.$$

$$E_{HF}(a_{eff})=E_{exact}(a).$$

$$\zeta (x)=0.395-1.138\arctan (0.362-0.994x),$$

$$H=\sum _{i=1}^N\left(-\frac{{\hslash }^2}{2m}{\nabla }_i^2+\frac{1}{2}m{\omega }_{ho}^2{r}_i^2\right)+\sum _{i<j}^N\frac{4\pi {\hslash }^2\zeta (k_{F}a)}{m{k}_F}\delta ({\text{r}}_{ij}).$$

$$[-\frac{{\hslash }^2}{2m}{\nabla }^2+\frac{1}{2}m{\omega }_{ho}^2{r}^2+\frac{4\pi (N-1){\hslash }^2}{3m}$$

$$\times \left({\zeta }^{\prime }(k_{F}a)a+2\frac{\zeta (k_{F}a)}{k_F}\right)|\psi {|}^2]\psi =ϵ\psi .$$

$$E=\int [N\psi \left(-\frac{{\hslash }^2}{2m}{\nabla }^2+\frac{1}{2}m{\omega }_{ho}^2{r}^2\right)\psi$$

$$+\frac{N(N-1)}{2}\frac{4\pi {\hslash }^2\zeta (k_{F}a)}{m{k}_F}{\psi }^4]d^3\text{r}.$$

$${\psi }_{TF}(\text{r})={\left[\frac{3(6{\pi }^2{)}^{1/3}(R_{TF}^2-r^2)}{16\pi N^{2/3}\zeta (+\infty )l_{ho}^4}\right]}^{3/4},$$

$$R_{TF}=N^{1/6}{\left(\frac{256\sqrt{2}}{9}\right)}^{1/6}{\left(\frac{\zeta (+\infty )}{\pi }\right)}^{1/4}{l}_{ho}.$$

$$ϵ=\frac{1}{2}m{\omega }_{ho}^2{R}_{TF}^2.$$

$$[-\frac{{\hslash }^2}{2m}{\nabla }^2+\frac{1}{2}m{\omega }_{ho}^2{r}^2+\frac{4\pi (N-1){\hslash }^2}{3m}$$

$$\times \left({\zeta }^{\prime }(k_{F}a)a+2\frac{\zeta (k_{F}a)}{k_F}\right)|\tilde{\psi }{|}^2]\tilde{\psi }=i\hslash \frac{\partial }{\partial t}\tilde{\psi }.$$

$${\tilde{\psi }}_{bog}(\text{r},t)=e^{-iϵt/\hslash }\left(\psi (\text{r})+u(\text{r})e^{-i\omega t}+v^{\ast }(\text{r})e^{i\omega t}\right).$$

$$\left(-\frac{{\hslash }^2}{2m}{\nabla }^2+\frac{1}{2}m{\omega }_{ho}^2{r}^2+f(N,a,\psi )-ϵ\right)u$$

$$+g(N,a,\psi )v=\hslash \omega u,$$

$$\left(-\frac{{\hslash }^2}{2m}{\nabla }^2+\frac{1}{2}m{\omega }_{ho}^2{r}^2+f(N,a,\psi )-ϵ\right)v$$

$$+g(N,a,\psi )u=-\hslash \omega v,$$

$$f\to \frac{8\pi N{\hslash }^{2}a}{m}{\psi }^2,g\to -\frac{4\pi N{\hslash }^{2}a}{m}{\psi }^2(a\to 0),$$

$$f\to \frac{40\pi N^{2/3}{\hslash }^2\zeta (+\infty )}{9(6{\pi }^2{)}^{1/3}m}{\psi }^{4/3},$$

$$g\to -\frac{16\pi N^{2/3}{\hslash }^2\zeta (+\infty )}{9(6{\pi }^2{)}^{1/3}m}{\psi }^{4/3}$$

$$R=\sqrt{\frac{1}{N}\sum _{i=1}^N{r}_i^2},$$

$$\tan {\alpha }_i=\frac{\sqrt{{\sum }_{j=1}^i{r}_j^2}}{r_{i+1}},(i=1,\cdots ,N-1)$$

$$H$$

$$V_{int}(R,\Omega )=\sum _{i<j}^N\frac{4\pi {\hslash }^2\zeta (k_{F}a)}{m{k}_F}\delta ({\text{r}}_{ij}).$$

$${\Lambda }^2{\Phi }_{\lambda }(\Omega )=\lambda (\lambda +3N-2){\Phi }_{\lambda }(\Omega ),$$

$$\left(-\frac{{\hslash }^2}{2M}\frac{d^2}{d{R}^2}+V_{eff}(R)-E\right)R^{(3N-1)/2}F(R)=0,$$

$$V_{eff}(R)=\frac{(3N-1)(3N-3){\hslash }^2}{8M{R}^2}+\frac{1}{2}M{\omega }_{ho}^2{R}^2$$

Full text

unitary Bose gas, renormalization, mean-field theory, hyperspherical coordinates

I Introduction

Strongly correlated systems near quantum degeneracy exhibit a wide range of intriguing phenomena. Paradigmatic examples include helium superfluidity and the fractional quantum Hall effect. In the atomic physics realm, the ultracold quantum gas, due to its simplicity, purity and high controllability, is an excellent candidate to be used to study strongly correlated systems. The interaction between cold atoms, which is typically characterized by the $s$-wave scattering length, can be readily controlled through Feshbach resonancesChin et al. (2010).

The Bose Einstein condensate (BEC) is a highly degenerate quantum system in which the interparticle interaction can also be tuned via a magnetic or other types of Feshbach resonanceInouye et al. (1998). When the scattering length in a BEC is much larger than any length scale of the system, the gas has reached the so-called unitary regimeCastin and Werner (2012). However, creating a BEC in the strongly interacting regime or even all the way to unitarity is extremely difficult. The major reason is that the three-body recombination rate at zero temperature in dilute gases is proportional to $a^{4}$Esry et al. (1999); Nielsen and Macek (1999); Bedaque et al. (2000); Braaten and Hammer (2006), which results in a very short lifetime of the strongly interacting Bose gas. This phenomenon contrasts sharply with the strongly interacting Fermi gas, for which the three-body recombination is suppressed by the Pauli exclusion principlePetrov et al. (2004). Because of the prohibitively high atom loss rate, it has been considered nearly impossible to access the unitary Bose gas adiabatically. However, a nonadiabatic approach to unitarity has been developed by the JILA groupMakotyn et al. (2014). In their experimental work, they studied the nonequilibrium dynamics of a degenerate unitary Bose gas and observed important universal properties of the system.

Although a few theories have been proposed to treat the degenerate unitary Bose gasCowell et al. (2002); Song and Zhou (2009); Diederix et al. (2011); Rossi et al. (2014); Zhou and Mashayekhi (2013); Werner and Castin (2006), no existing theories so far are capable of completely describing this system. Some theories involve complex derivations such as the renormlization group theoryZhou and Mashayekhi (2013) or extensive computations like Monte Carlo simulationRossi et al. (2014). Shortly after the JILA experiment on the unitary Bose gas, various theoretical descriptions were proposed in an effort to explain the experimental results, especially the momentum distributionSmith et al. (2014); Sykes et al. (2014); Laurent et al. (2014); Kira (2015); Comparin and Krauth (2016).

In this article, we introduce a renormalized contact potential similar to that in Ref. von Stecher and Greene (2007), to extend the validity of the zero-range potential to the strongly interacting regime. Then we employ this renormalized potential in company with traditional many-body theories to study the degenerate unitary Bose gas at zero temperature. The structure of this article is as follows: In section II, we elaborate the renormalization procedure and the physical ideas behind it. After that, we apply this renormalized potential to a few many-body theories and show how they are modified with the inclusion of the renormalization. In section III, we discuss the stationary properties and elementary excitations of a degenerate unitary Bose gas using our renormalization theory, and compare our results with other theoretical predictions. We particularly focus on a few important physical observables of the system. Finally, in section IV, we summarize our work and the most significant findings.

II Theory of renormalization

Our renormalization is similar to that in a two-component Fermi gasvon Stecher and Greene (2007). Here we summarize the procedure and the physical origin of this renormalization. For a uniform gas system, when the range of the two-body interaction is much smaller than both the scattering length $a$ and the average interparticle distance determined by the particle density $n$, the behavior of the system is characterized by the dimensionless parameter $na^{3}$. This is also equivalent to the dimensionless parameter $k_{F}a$, where $k_{F}=(6\pi^{2}n)^{1/3}$ is defined for the Bosonic system in a manner akin to the Fermi momentum. Our idea is to design an effective scattering length $a_{eff}$ that can replace the bare scattering length $a$ to describe the properties of the system. In this case, the dimensionless parameter becomes $k_{F}a_{eff}$. Therefore, there must be a correspondence between $k_{F}a_{eff}$ and $k_{F}a$ characterized by a renormalization function $k_{F}a_{eff}=\zeta(k_{F}a)$. The effective scattering length is designed specifically for a system of two interacting particles in a harmonic trap such that it can exactly describe the atomic ground state energy.

For two particles in a harmonic trap with a circular frequency $\omega_{ho}$ interacting with a regularized pseudo potential

[TABLE]

the Hamiltonian is given by

[TABLE]

where $\textbf{r}_{12}=\textbf{r}_{1}-\textbf{r}_{2}$ is the relative coordinate between two particles. The wave function is separable in the center of mass motion and the relative motion, that is, $\Psi=\psi_{cm}(\textbf{R}_{cm})\psi_{rel}(\textbf{r}_{12})$, where $\textbf{R}_{cm}=(\textbf{r}_{1}+\textbf{r}_{2})/2$ is the center of mass coordinate. The exact solution to the corresponding Schrödinger equation $H_{2b}\Psi=E_{exact}(a)\Psi$ has been discussed in Ref. Busch et al. (1998). The eigen-energy can be written as $E_{exact}(a)=E_{cm}+E_{rel}(a)$, where $E_{cm}=(n_{cm}+3/2)\hbar\omega_{ho}$ corresponds to the center of mass motion in the harmonic trap. The eigen-energy for the relative motion satisfies the following condition:

[TABLE]

where $l_{ho}=\sqrt{\hbar/m\omega_{ho}}$ is the harmonic trap length.

On the other hand, when the two particles interact with the renormalized contact potential given by

[TABLE]

we assume the total wave function has a Hartree-Fock (HF) expression $\tilde{\Psi}=\psi(\textbf{r}_{1})\psi(\textbf{r}_{2})$. Consequently, the energy expectation value of the system is given by

[TABLE]

Since there is no Fermi momentum in few body systems, it is natural to replace $k_{F}$ in Eq. (4) by its average value $\langle k_{F}\rangle=\int[6\pi^{2}2\psi(\textbf{r})^{2}]^{1/3}\psi(\textbf{r})^{2}d^{3}\textbf{r}$ in this two body case. Minimizing $\xi\{\psi\}$ with the normalization constraint $\langle\psi|\psi\rangle=1$ yields the ground state energy that depends on the effective scattering length:

[TABLE]

In order to make the effective scattering length and the bare scattering length equivalent for this trapped two-particle system, we match the HF energy to the exact energy:

[TABLE]

Before matching these two energies, we should note that the exact energy has many branches including a molecular branch for $a>0$. Because the HF approximation describes the lowest atomic gas state, which corresponds to the branch with $n_{cm}=0$ and $(1/2)\hbar\omega_{ho}<E_{rel}<(5/2)\hbar\omega_{ho}$, we match the HF energy to the exact energy in this particular branch. Since Eq. (7) yields a pointwise correspondence between $\langle k_{F}\rangle a_{eff}$ and $\langle k_{F}\rangle a$, we can numerically interpolate the renormalization function $\langle k_{F}\rangle a_{eff}=\zeta(\langle k_{F}a\rangle a)$. This interpolation can be excellently fitted by an analytical expression

[TABLE]

which satisfies the asymptotic conditions $\zeta(+\infty)=2.182$, $\zeta(-\infty)=-1.392$ and $\zeta(k_{F}a)\rightarrow k_{F}a$ for $|k_{F}a|\ll 1$.

From the renormalization procedure above, we can see that the renormalized contact potential Eq. (4) reproduces the exact energy solution for the system of two interacting particles in a trap. The next step is to apply such a renormalized contact potential to many body systems. Since the many body Hamiltonian cannot be diagonalized exactly due to the huge number of degrees of freedom, we must make some aggressive but reasonable approximations, as we will discuss in the following subsections.

II.1 Mean-field approach

One natural and intuitive idea is to generalize the HF approximation employed above along with the renormalized interaction to many body systems. Such an approximation for bosons is also called the mean-field approximation.

With the inclusion of renormalized interactions, the $N$-body Hamiltonian now becomes

[TABLE]

The mean-field theory assumes the $N$-body ground state wave function to be $\Psi=\prod_{i=1}^{N}\psi(\textbf{r}_{i})$. By taking the variation $\delta H/\delta\psi=0$ under the normalization condition $\langle\psi|\psi\rangle=1$, we can obtain a renormalized $N$-body Gross-Pitaevskii (GP) equation:

[TABLE]

where $\epsilon$ is the Lagrange multiplier enforcing normalization in the variation procedure, and is also identified as the orbital energy. $\zeta^{\prime}(x)$ means the derivative of $\zeta$ with respect to the variable $x$. $k_{F}$ is the local Fermi momentum. With such a mean-field approximation and in the framework of the local density approximation (LDA), the local Fermi momentum is given by $k_{F}=(6\pi^{2}N|\psi|^{2})^{1/3}$. After solving Eq. (10), we can evaluate the total energy of the system as

[TABLE]

At very large particle numbers, the kinetic energy composes a very small portion of the total energy of the system. In the meanwhile, the particle density varies slowly in the trap except for the edge of the cloud, which is an evidence of the validity of the LDA. In this condition, we can solve Eq. (10) using the Thomas-Fermi approximation, which neglects the kinetic energy term and thereby converts Eq. (10) to a regular algebraic equation. In the unitary regime, if we neglect edge effects, which means assuming $k_{F}a\rightarrow+\infty$ at any position of the cloud, we can obtain an analytical expression for the wave function, which is given by

[TABLE]

where $R_{TF}$ is the Thomas-Fermi radius indicating the size of the cloud. It is given by

[TABLE]

Consequently, the orbital energy is given by

[TABLE]

The fact that $\psi_{TF}$ does not depend on $a$ is another signature that the scattering length is no longer a relevant length scale in the unitary regime.

In order to calculate the dynamics of a degenerate Bose gas, it is natural to convert Eq. (10) to a time-dependent GP equation:

[TABLE]

We should notice that now the local Fermi momentum $k_{F}$ also becomes time-dependent. Although such a non-linear and time-dependent Shrödinger equation can be solved directly using a brute force time-evolution method, we can obtain a clearer physical picture if appropriate approximations are made to Eq. (15). One typical method is the Bogoliubov approximation, which is commonly used to predict elementary excitations of a BEC. The Bogoliubov approximation assumes the time-dependent wave function to be

[TABLE]

Inserting this wave function into Eq. (15) and linearizing the equation to the first order in $u(\textbf{r})$ and $v(\textbf{r})$, we can obtain a pair of coupled differential equations:

[TABLE]

[TABLE]

where $\omega$ corresponds to the eigen mode frequency. $f$ and $g$ are in general complicated functions of $N$, $a$, and $\psi$, but they have simple forms in the asymptotic limits $a\rightarrow 0$ and $a\rightarrow\infty$:

[TABLE]

and

[TABLE]

The lowest and most significant eigen-mode is called the breathing mode, which corresponds to the oscillation of the overall size of the cloud with a fixed geometry. We will discuss this mode in later sections.

II.2 Hyperspherical description

The hyperspherical coordinate system is a powerful toolkit to treat few-body problemsvon Stecher et al. (2009); Wang et al. (2012a, b); Rittenhouse et al. (2011). To generalize this toolkit to many body systems, aggressive approximations must be made to significantly reduce the dimension of the problem. The hyperspherical description of a single component weakly-interacting BEC and a degenerate Fermi gas have been studied with the bare Fermi pseudo potential, and important physics has been predicted even with crude approximationsBohn et al. (1998); Rittenhouse et al. (2006). Similar to the procedure in Ref. Bohn et al. (1998), we formulate the hyperspherical theory in a degenerate unitary Bose gas with the renormalized interaction.

For $N$ particles in a 3D space, which contains $3N$ degrees of freedom, the hyperspherical coordinates are constructed as follows: The hyperradius $R$, which is a collective coordinate and indicates the overall size of the cloud, is given by

[TABLE]

The remaining $3N-1$ coordinates are called hyperangles. $2N$ of them are the defined as the regular spherical angles of the $N$ particles, that is, $\{\theta_{1},\phi_{1},\theta_{2},\phi_{2},\cdots,\theta_{N},\phi_{N}\}$. The remaining $N-1$ hyperangles can be defined as

[TABLE]

With this set of hyperspherical coordinates, the Hamiltonian in Eq.(9) can be rewritten as

[TABLE]

where $M=Nm$ is the total mass of the $N$ particles. $V_{int}(R,\Omega)$ denotes the renormalized interactions written in hyperspherical coordinates, that is,

[TABLE]

$\Lambda$ is called the grand angular momentum operator. The eigen-functions of the operator $\Lambda^{2}$, denoted by $\Phi_{\lambda}$, are called hyperspherical harmonicsAvery (1989). They satisfy the equation

[TABLE]

where $\Omega$ represents all hyperangles. For a given $\lambda$, $\Phi_{\lambda}$ usually has huge degeneracy especially for large $\lambda$. The eigen-functions $\Phi_{\lambda}$ form a basis in the hyperangular Hilbert space. An aggressive approximation we make here is to only retain one hyperangular momentum eigen-state out of this huge basis set. This is also known as the $K$-harmonics approximation in nuclear theoriesSmirnov and Shitikova (1977) and it becomes exact in the unitary limit and in the non-interacting limitWerner and Castin (2006); Blume et al. (2008). The natural choice of this eigen state for bosons would be the lowest eigen state of $\Lambda^{2}$, denoted by $\Phi_{0}$, which in fact is a constant. Such a choice of hyperangular wave function also freezes the geometry of the atomic cloud into that of the non-interacting case, while the interparticle interactions modify the overall size of the cloud, which is reflected in the hyperradial wave function. With this $K$-harmonics approximation, the total $N$-body wave function can be separated as $\Psi(R,\Omega)=F(R)\Phi_{0}(\Omega)$. Inserting this expression into the time-independent Schrödinger equation and integrating over all hyperangles yields a hyperradial Schrödinger equation:

[TABLE]

where $V_{eff}(R)$ is an effective hyperradial potential written as

[TABLE]

where $\langle\cdots\rangle_{\Omega}$ denotes integration over all hyperangles. The evaluation of the last term in $V_{eff}$ has been elaborated in Ref. Bohn et al. (1998). In the large $N$ limit which we are interested in, $\langle\Phi_{0}|V_{int}(R,\Omega)|\Phi_{0}\rangle_{\Omega}$ can be readily calculated numerically. We can even obtain analytical expressions in the unitary limit and the weakly interacting limit: For $a\rightarrow\infty$,

[TABLE]

and for $a\rightarrow 0$,

[TABLE]

Fig. 1 shows the effective hyperradial potential curves at different scattering lengths. At large hyperradius, $V_{eff}(R)$ is always dominated by the $R^{2}$ term representing the confinement of the harmonic trap. At a small hyperradius, the system feels two repulsive forces originating from the quantum pressure and the two-body interaction. It is interesting that the two-body interaction term transitions from $R^{-3}$ to $R^{-2}$ as the scattering length increases to infinity.

Eq. (26) is equivalent to the Schrödinger equation of a particle moving in a 1D potential. Moreover, in the large $N$ limit, the mass of this “particle” is so huge that it can be treated classically. Consequently, the minimum of $V_{eff}$ corresponds to the total energy of the system at equilibrium, that is,

[TABLE]

where $R_{0}$ denotes the equilibrium position. Furthermore, as the hyperangular wave function is kept frozen, the oscillation of this massive “particle” in the hyperradial potential indicates the oscillation of the overall size of the system, which corresponds to the breathing mode. The breathing mode frequency is associated with the coefficient of the second order expansion of $V_{eff}$ at $R_{0}$, that is,

[TABLE]

III Results and discussions

At first, we discuss the total energy of a degenerate Bose gas. Fig. 2 shows the average energy per particle for $N=10^{4}$ particles as a function of the scattering length. In the weakly interacting regime where $\langle n\rangle a^{3}\ll 1$, the mean-field energy obtained by solving the GP equation with renormalization agrees excellently with that without renormalization. These two results start to separate near $a/l_{ho}=0.2$, which corresponds to $\langle n\rangle a^{3}=0.25$. The mean-field energy diverges at unitarity without renormalization. The energy obtained using the Thomas-Fermi approximation and the mean-field energy differ in the weakly interacting regime where the kinetic energy is a significant contribution to the total energy. However, at large scattering lengths, the total energy of the system is dominated by the strong interactions between particles and thereby the Thomas-Fermi approximation agrees excellently with the mean-field result. The energy obtained using the hyperspherical method agrees qualitatively with the mean-field result, though it is slightly smaller. Overall, the total energy of the system saturates at large scattering lengths with the inclusion of the interaction renormalization, which indicates that the scattering length is no longer a relevant length scale of the system near unitarity.

We now discuss the energy of a unitary Bose gas using the Thomas-Fermi approximation since it is very accurate in the strongly interacting regime. One advantage of the Thomas-Fermi approximation is that we can obtain analytical expressions for many physical quantities, which offers us a clear picture of the unitary Bose gas. With the Thomas-Fermi approximation, the ground state energy is given by

[TABLE]

For a unitary gas in a uniform space, the only relevant length scale of the system is the average interparticle distance determined by $n^{-1/3}$, where $n$ is the particle density. This also defines the only energy scale of the system $\hbar^{2}n^{2/3}/2m$. When the gas is inhomogeneous while the density varies slowly in space, the local density approximation can be applied to the system and thereby $n^{2/3}$ is replaced by its average value $\langle n^{2/3}\rangle$. For a unitary Bose gas in a harmonic trap, the average value is given by

[TABLE]

Therefore, from Eq. (32) and Eq. (33), we can obtain the universal relation of the energy of a unitary Bose gas, which is given by

[TABLE]

This universal relation is close to the value $E/N=13.33\hbar^{2}n^{2/3}/2m$ reported in Ref. Cowell et al. (2002).

We show in Fig. 3 the average energy per particle, in units of $\hbar^{2}\langle n^{2/3}\rangle/2m$, as a function of the average density for different scattering lengths to verify this universal relation. The results are obtained by solving the renormalized GP equation Eq. (10). At small scattering length $a/l_{ho}=0.1$, the value of $2mE/N\hbar^{2}\langle n^{2/3}\rangle$ varies significantly with the average density. As the scattering length increases, $2mE/N\hbar^{2}\langle n^{2/3}\rangle$ has weaker dependence on the average density and the value approaches the universal constant in Eq. (34). At $a/l_{ho}=10$, where gas has reached the unitary regime, the result agrees excellently with the universal relation Eq. (34) for $\langle n\rangle l_{ho}^{3}>5$, The small deviation from the universality at small densities may be due to the inaccuracy of LDA when the interparticle distance is comparable to the trapping length.

Besides the total energy, there are many interesting physical quantities worthy of investigation in unitary Bose gases. An important quantity that bridges the two-body correlations and the thermodynamics of a many-body system is called the two-body contact or Tan’s contact. It was first introduced by Shina Tan to study the universal properties of a two-component Fermi gas with $s$-wave contact interactionsTan (2008a, b). Universal relations determined by the two-body contact have also been identified in systems consisting of identical bosonsBraaten et al. (2011). The two-body contact in bosons has been measured using rf spectroscopyWild et al. (2012). The two-body contact is determined by the derivative of the total energy of the system with respect to the scattering length:

[TABLE]

It is an extensive thermodynamic quantity of the system. Another intrinsic quantity, which is commonly used in homogeneous systems, is the contact density $\mathcal{C}_{2}$, which can be obtained from the limit of the high momentum tail: $\mathcal{C}_{2}=\lim_{k\rightarrow\infty}k^{4}n_{\textbf{k}}$, where $n_{\textbf{k}}$ is the number of particles in the k momentum state. Since the interparticle distance is the only relevant length scale of a homogeneous system at unitarity, the two-body contact density must scale as

[TABLE]

where $\alpha$ is an universal dimensionless coefficient. Such a universal relation can be generalized to a trapped system under LDA, which is given by

[TABLE]

Fig. 4 shows the average two-body contact per particle, in units of $\langle n^{1/3}\rangle$, as a function of the scattering length. The results are obtained by directly solving the renormalized GP equation Eq. (10). The value of $C_{2}/(N\langle n^{1/3}\rangle)$ has similar behavior for different numbers of particles: It increases drastically in the weakly interacting regime as the scattering length grows; it then attains a maximum value before saturating in the unitary regime. For different numbers of particles $C_{2}/(N\langle n^{1/3}\rangle)$ saturates at the same value, demonstrating the universality of $\alpha$.

From our renormalized mean-field model, we can derive the universal relation Eq. (37) and determine the universal coefficient $\alpha$ analytically with appropriate approximations. To calculate the two-body contact, we take the derivative of the renormalized mean-field energy Eq. (11) with respect to the scattering length. We should note that at unitarity the scattering length dependence of the wave function $\psi$ is much weaker than the renormalization function $\zeta(k_{F}a)$. Thus, by approximating $\partial\psi/\partial a|_{a\rightarrow\infty}\approx 0$ and neglecting the edge effect, we can readily obtain Eq. (37) and the coefficient $\alpha$ is given by

[TABLE]

Other theoretical works have reported the the universal coefficient to be $\alpha=10.3$Diederix et al. (2011), $9.04$Rossi et al. (2014) and $12$Sykes et al. (2014), which agree qualitatively with our prediction from the renormalized mean-field approach.

To further verify our prediction of the universal relation Eq. (37) and the value of the coefficient $\alpha$, we show in Fig. 5 the average two-body contact per particle, in units of $\langle n^{1/3}\rangle$, as a function of the average particle density for different scattering lengths. At small scattering length $a/l_{ho}=0.1$, $C_{2}/(N\langle n^{1/3}\rangle)$ increases significantly with the density. The value of $C_{2}/(N\langle n^{1/3}\rangle)$ at $a/l_{ho}=1$ is larger than its values for both $a/l_{ho}=0.1$ case and $a/l_{ho}=10$ case because it attains a maximum with the increment of the scattering length, as shown in Fig. 4. In the unitary regime, $C_{2}/(N\langle n^{1/3}\rangle)$ converges to the universal coefficient $\alpha$, which corresponds to $a/l_{ho}=10$ case. The value of $\alpha$ from exact calculation is slightly smaller than the the analytical approximation $\alpha\approx 11.8$, which might be due to the edge effect and the small scattering length dependence of the wave function.

Finally, we discuss the elementary excitations of the degenerate Bose gas in a harmonic trap. Specifically, we focus on the lowest radial excitation, which corresponds to the breathing mode. The breathing mode frequency can be determined from the hyperspherical method, the Bogoliubov approximation, or by directly solving the time-dependent GP equation, as discussed in the theory of renormalization section above. Fig. 6 shows the breathing mode frequency as a function of the scattering length for $N=10^{4}$ particles. These three different methods show overall consistency with each other. The time-evolution results are obtained by solving Eq. (15) with an interaction quench. At vanishing scattering length, the breathing mode frequency is apparently $\omega=2\omega_{ho}$, which corresponds to the beating between two adjacent harmonic levels with zero angular momentum. The breathing mode frequency increases gradually with the scattering length before it reaches a maximum value. It is interesting that the breathing mode frequency regresses to the non-interacting value $\omega=2\omega_{ho}$ at unitarity. Similar unitary behavior has also been predicted in Ref. Werner and Castin (2006).

IV Conclusion

In summary, we introduced a renormalized contact potential to study degenerate Bose gases with large scattering lengths. Such a renormalized interaction is designed by matching the Hartree-Fock energy to the exact energy of two interacting particles in a trap. We employ this renormalized contact potential in company with the mean-field theory and the hyperspherical theory to study the stationary properties and elementary excitations of a degenerate Bose gas, especially in the unitary regime. In the framework of the local density approximation, the only relevant length scale of a degenerate unitary Bose gas is the interparticle spacing $n^{-1/3}$, where $n$ is the particle density. This length scale also determines the only energy scale $\hbar^{2}n^{2/3}/2m$ and the only two-body contact scale $n^{1/3}$ of the system. Our renormalization theory offers us a much more clear and convenient approach to obtain the universal relations for energy and two-body contact at unitarity, which are given by $E/N=12.67\hbar^{2}\langle n^{2/3}\rangle/2m$ and $C_{2}/N=11.8\langle n^{1/3}\rangle$ respectively. Our results are in consistent with other theoretical predictions. Moreover, we studied the lowest radial excitation of a degenerate Bose gas, which can be induced by an interaction quench. This excitation is also known as the breathing mode. Our theory shows an interesting phenomenon that the breathing mode frequency at unitarity returns to the value of a non-interacting gas.

Acknowledgements.

The authors thank Q. Zhou and F. Robicheaux for helpful discussions, and M. Eiles for careful reading of the manuscript. This work is supported by National Science Foundation and US-Israel Binational Science Foundation.