Inhomogeneous Bose Gas

File:Superfluid to Mott insulator transition-bec3-first-expt

The physics of a BEC happens not just in momentum space, but also in position space, and it is useful to analyze it accordingly. With a trapping potential applied, the Hamiltonian is written in terms of bosonic field operators ${\hat {\psi }},{\hat {\psi }}^{\dagger }$ (obeying $[{\hat {\psi }}(r),{\hat {\psi }}^{\dagger }(r')]=\delta (r-r'),[{\hat {\psi }}(r)^{(\cdot )/\dagger },{\hat {\psi }}^{(\cdot )/\dagger }(r')]=0$ )

{\hat {H}}=\int d^{3}r\,{\hat {\psi }}^{\dagger }(r){\Big [}{\frac {-\hbar ^{2}}{2m}}\nabla ^{2}+V_{trap}{\Big ]}{\hat {\psi }}(r)+{\frac {1}{2}}\int d^{3}r\int d^{3}r'{\hat {\psi }}^{\dagger }(r){\hat {\psi }}^{\dagger }(r^{\prime })U(r-r^{\prime }){\hat {\psi }}(r^{\prime }){\hat {\psi }}(r)

Note: this is a general field-quantized expression for a Hamiltonian with two-body interaction. If you are not familiar with second quantization, consult the first chapter of "Quantum Theory of Many-Particle Systems" by Fetter and Walecka, which reviews the mapping between first quantization and second quantization in detail.

We are interested in the time-evolution of the operator ${\hat {\psi }}$ (i.e. we work in the Heisenberg picture). The equation of motion is given by the Heisenberg equation

i\hbar {\frac {\partial }{\partial t}}{\hat {\psi }}=[{\hat {\psi }},{\hat {H}}]

Using the field operator commutation relation, we can write the right-hand side of the equation as

-{\frac {\hbar ^{2}}{2m}}\nabla ^{2}{\hat {\psi }}(r,t)+V_{trap}{\hat {\psi }}(r,t)+{\frac {1}{2}}{\Big (}\int d^{3}r'\,{\hat {\psi }}^{\dagger }(r',t)U(r-r'){\hat {\psi }}(r',t)+\int d^{3}r'\,{\hat {\psi }}^{\dagger }(r',t)U(r'-r){\hat {\psi }}(r',t){\Big )}{\hat {\psi }}(r,t)

This must be approximated, in the spirit of Bogoliubov's momentum space approximation, to obtain a useful solution. We thus replace

{\hat {\psi }}(r,t)=\psi (r,t)+{\hat {\delta \psi }}(r,t)\,,

where the complex number $\psi (r,t)$ is an expectation (mean field), and the operator ${\hat {\delta \psi }}(r,t)$ captures the quantum (+ thermal) fluctuations. We do not consider terms proportional to ${\hat {\delta \psi }}$ . We further assume that the interaction potential is a delta function, $U(r-r')=U_{0}\delta (r-r')$ (which is valid for s-wave scattering at short range). With these assumptions, we obtain a nonlinear Schrodinger equation:

i\hbar {\frac {\partial \psi }{\partial t}}=\left[{-{\frac {\hbar ^{2}}{2m}}\nabla ^{2}+V_{trap}+U_{0}|\psi (r,t)|^{2}}\right]\psi (r,t)

This is known as the Gross-Pitaevskii equation.

The $|\psi (r,t)|^{2}$ term adds an energy proportional to the density due to the interactions. That density-dependence implies that increasing the density comes with an energy cost and decreasing density lowers the energy. Since the density is multiplied by the wave function, there is a non-linear energy dependence on the density. Therefore, lowering the density in one region and raising it in another costs energy and so this term makes the condensate try to have a uniform density distribution (the total number of particles is fixed so the integrated density cannot change).

The density-dependent term has an interesting interplay with the other two terms. Consider its effect on a condensate in a box potential for instance. If it were not for the interactions, the BEC density distribution would simply be set by the ground state of the box potential. However, the central region of that distribution is very dense, so the density-dependent term will try to flatten out the distribution and push atoms out towards the wings. It can flatten out the distribution near the center very well, but the trap potential makes it impossible to push atoms out too far, so the atom distribution stays roughly flat over a finite region. At the edge of the flat region of the density distribution are the wings which must go to zero at the edges of the box. The density energy cost makes the cloud want to make this crossover region as small as possible. However, making it very small would take a lot of curvature in the wave function, causing the kinetic energy term to increase. Therefore the distribution at the edges reaches a compromise between interaction energy and kinetic energy curvature. The length scale of the crossover region set by this compromise (assuming repulsive interactions) is known as the healing length $\xi$ ,

\xi =(8\pi a_{s}n)^{-1/2}\,,

arising from

{\frac {\hbar ^{2}}{2m\xi ^{2}}}=U_{0}n={\frac {4\pi \hbar ^{2}}{m}}a_{s}n\,.

where $a_{s}$ is the s-wave scattering length.

We remarked that the Gross-Pitaevskii equation looks like a nonlinear Schrodinger equation. Then what should its time-independent version look like? In other words, what is the $E$ that appears in $\psi (t)=\psi \exp(-iEt/\hbar )$ ? For ordinary time-independent Schrodinger equation applied to wavefunctions, that would be the energy of an eigenstate. But the notion of an eigenfunction for a nonlinear equation is not clear. The answer is the chemical potential $\mu$ , and it can be reasoned as follows (from "Bose-Einstein Condensation and Superfluidity" by Pitaevskii and Stringari). Our assumption ${\hat {\psi }}\sim \langle {\hat {\psi }}\rangle \neq 0$ requires that the bra-state and the ket-state cannot have the same number of particles - the condensate acts like a reservoir. When the number of particles is large and the interaction is not too strong, adding one particle only makes correction to the state on the order of $1/N_{0}$ , where $N_{0}$ is the average number of particles. So in this limit, it is valid to interpret $\langle {\hat {\psi }}\rangle$ so that the bra-state has one number of particle less than the ket state. Now if we take an average over stationary states, whose time dependence is governed by $\exp(-iE(N)t/\hbar )$ , the time dependence of the mean value ${\hat {\psi }}=\psi$ is given by $\exp(-i[E(N)-E(N-1)]t/\hbar )\approx \exp(-i(\partial E/\partial N)t/\hbar )=\exp(-i\mu t/\hbar )$ . The quantity $\langle {\hat {\psi }}\rangle$ is known as the superfluid order parameter in statistical physics. The non-zero, complex value of the superfluid parameter distinguishes the superfluid phase from the normal phase. The term order parameter is used to describe a phase with broken symmetry, which we will explain later in the context of BEC.

In summary, the time-independent version of the GPE reads

\mu =\left[{-{\frac {\hbar ^{2}}{2m}}\nabla ^{2}+V_{trap}+U_{0}|\psi (r)|^{2}}\right]\psi (r)

It is worth emphasizing that the superfluid order parameter is not equivalent to an ordinary Schrodinger wavefunction (from Pitaevskii and Stringari). First of all, we already saw that the order parameter evolves in time with chemical potential, not energy. Also, because the GPE is non-linear, two solutions $\psi _{a}$ and $\psi _{b}$ with different chemical potentials are not necessarily orthogonal ( $N^{-1}\int d^{3}r\,\psi _{a}^{*}\psi _{b}\neq 0$ ). But quantum mechanics is linear and hence different eigenstates are orthogonal. How should these two facts be reconciled? In the weak-interacting limit, so that we can ignore correlations between particles (which the mean-field approximation implies), we can write the many-body wavefunction $\Phi$ corresponding to the field $\psi$ in a factorized form (Hartree-Fock approximation):

\Phi (r_{1},\ldots ,r_{N})=\prod _{i=1}^{N}{\frac {1}{\sqrt {N}}}\psi (r_{i})

The inner product of two different many-body wavefunctions corresponding to $\psi _{a}$ and $\psi _{b}$ is $(N^{-1}\int d^{3}r\psi _{a}^{*}\psi _{b})^{N}$ , which goes to zero unless $\psi _{a}=\psi _{b}$ when $N\to \infty$ . This also emphasizes that the GPE and the mean-field picture is only valid in the limit of large $N$ .

File:Superfluid to Mott insulator transition-bec3-healing

If the interactions are really strong, the kinetic energy term can be neglected, because the interactions will keep the density constant in its spatial distribution. This is a particularly good approximation in the flat center of the distribution where there is little wave function curvature. Such an approximation is the Thomas-Fermi approximation, giving an equation for the wavefunction,

[V_{ext}+U_{0}|\psi (r,t)|^{2}-\mu ]\psi (r)=0\,,

giving the solution

|\psi (r,t)|^{2}={\frac {1}{U_{0}}}\left(\mu -V(r)\right)\,.

The wavefunction is essentially just the potential filled up to the chemical potential level, inverted. For a quadaratic potential, $V(r)=m\omega ^{2}r^{2}/2$ , the chemical potential is

\mu ={\frac {\hbar \omega }{2}}x^{2/5}\,,

where $x=15Na/a_{osc}$ is a common term worth identifying, and $a_{osc}={\sqrt {\hbar ^{2}/m\omega }}$ is a characteristic length scale of the oscillator, its zero point motion. Defining $R=a_{osc}x^{1/5}$ , we may find $\mu =m\omega ^{2}R^{2}/2$ . This explains the profile of the condensate data obtained in experiments:

File:Superfluid to Mott insulator transition-bec3-bec-peak

Note that the size of the ground stat BEC is much larger than the zero-point motion of the harmonic oscillator. This is due to the pressure of the repulsive interactions. The Gross-Pitaevskii interaction gives not only the ground state wavefunction, but also the dynamics of the system. For example, it predicts soliton formation: stable wavefunctions with a size scale determined by a balance of the kinetic energy and the internal interactions. This requires, however, an attractive potential. Such soliton formation can nevertheless be seen in BEC's, with tight traps (see recent Paris experiments).

Length and energy scales in BEC

Size of atom: $a=3$ nm
Separation between atoms $n^{1/3}=200$ nm
Matter wavelength $\lambda _{dB}=1$ $\mu$ m
Size of confinement $a_{osc}=30$ $\mu$ m

Note that

a\ll n^{1/3}\leq \lambda _{dB}<2\pi \xi <a_{osc}\,.

For a gas, $a\ll n^{1/3}$ . For a BEC, $\lambda _{dB}\geq n^{1/3}$ in addition. The corresponding energy scales are also useful to identify. Let $\hbar \omega =a_{osc}$ . Then: