text:
# Localized Bosonic Modes in Superconductors

## Abstract

We show that a localized bosonic mode acts as a new type of “defect” in $`s`$\- and $`d_{x^2y^2}`$-wave superconductors. The mode induces bound or resonance states, whose spectral signature are peaks in the superconductor’s density of states (DOS). We study the peaks’ shape and energy as a function of temperature and the mode’s frequency and lifetime. We identify several characteristic signatures of the localized mode that qualitatively distinguishes its effects from those of magnetic or non-magnetic impurities.

The study of local defects or impurities in superconductors has attracted significant experimental Yaz97 ; exp  and theoretical theo1 ; theo2 ; Bal96  interest over the last few years. These studies have proven particularly important for further elucidating the nature of the superconducting (SC) pairing mechanism in unconventional superconductors. In particular, recent scanning tunnelling microscopy (STM) experiments provided a detailed picture of the frequency and spatial dependence of defect induced resonance states in the high-temperature superconductor (HTSC) Bi2Sr2CaCu2O8+δ exp , and the triplet superconductor Sr2RuO4 Davispc . Several theoretical scenarios for the physical origin of these impurity states have been proposed, ranging from electronic scattering off classical impurities theo1 ; Bal96  to the onset of Kondo-screening theo2 .

In this Letter we investigate a new type of “defect”, a localized bosonic mode, and study its effect on the local electronic structure in $`s`$\- and $`d_{x^2y^2}`$-wave superconductors. The study of localized modes, and in particular phonon modes in normal metals, has been of great interest in the context of molecular spectroscopy Jak66 ; review ; Sca67  (or inelastic electron tunnelling). It was recently argued that a localized magnetic mode could also arise from unscreened magnetic moments in a $`d_{x^2y^2}`$-wave superconductor Bal02 . The nature of the mode is only relevant for the present study to the extent that it determines the form of its coupling to the electronic degrees of freedom. While theoretical studies Sca67 ; Bal02  have so far only focused on the Born-limit of weak electronic scattering off the localized mode, we show that the strong-scattering, unitary limit gives rise to an abundance of novel phenomena. In particular, we show that the presence of a localized mode leads to the appearance of fermionic bound or resonance states whose spectral signature are peaks in the SC DOS. We study the peaks’ energy and shape as a function of the mode’s characteristic frequency $`\omega _0`$, lifetime $`\mathrm{\Gamma }^1`$, and temperature. It is of particular interest that the peaks move to higher energies with increasing $`\omega _0`$ or $`\mathrm{\Gamma }`$, but are shifted to lower energies with increasing temperature, $`T`$. In addition, the temperature evolution of the peaks differs qualitatively if its frequency is smaller or larger than $`\mathrm{\Delta }_0\omega _0`$, where $`\mathrm{\Delta }_0`$ is the maximum SC gap. Moreover, we show the mode also induces a “dip” in the DOS at frequencies $`\pm (\mathrm{\Delta }_0+\omega _0)`$. This dip, together with the temperature evolution of the DOS, qualitatively distinguishes the effects of a localized mode from those of non-magnetic or static magnetic impurities theo1 ; Bal96 . These characteristic differences therefore provide an important tool for future experiments to further study and clarify the nature of “defects”.

Starting point for our calculations is the Hamiltonian

$`H`$ $`=`$ $`{\displaystyle \underset{𝐤}{}}ϵ_𝐤c_{𝐤,\sigma }^{}c_{𝐤,\sigma }+{\displaystyle \underset{𝐤}{}}\mathrm{\Delta }_𝐤c_{𝐤,}^{}c_{𝐤,}^{}+h.c.`$ (1)
$`+\omega _0b_0^{}b_0+{\displaystyle \underset{𝐤,𝐪,\sigma }{}}g_\sigma \left(b_0^{}+b_0\right)c_{𝐤,\sigma }^{}c_{𝐪,\sigma }.`$
Here, $`c_k^{},b_0^{}`$ are the fermionic and bosonic creation operators, respectively. We consider a two-dimensional electronic system with normal state dispersion $`ϵ_𝐤`$ and SC gap $`\mathrm{\Delta }_𝐤`$. The bosonic mode with frequency $`\omega _0`$ is localized at $`𝐫=0`$. We allow for a spin dependence of the scattering vertex, $`g_\sigma `$, and for simplicity consider only an on-site interaction. In the following, we study a single bosonic mode with no internal degrees of freedom, i.e., spin $`S=0`$. The extension of our results to $`S1`$, in which case electronic spin-flip scattering is allowed, is straightforward and will be discussed in a future publication Nyb03 . The retarded bosonic Greens function $`D_R(\omega )`$ in frequency space is given by

$$D_R(\omega )=\frac{\alpha }{(\omega +i\mathrm{\Gamma }/2)^2\omega _0^2}.$$
(2)
where $`\alpha `$ depends on $`\mathrm{\Gamma }`$ and $`\omega _0`$, and is determined via

$$1+2n_b(\omega _0)=_{\mathrm{}}^{\mathrm{}}\frac{d\omega }{\pi }n_b(\omega _0)\mathrm{Im}D_R(\omega ).$$
(3)
We assume that $`\omega _0`$ and the mode’s lifetime, $`\mathrm{\Gamma }^1`$, are determined intrinsically such as in the case of phonon modes Jak66 . Possible feedback effects on $`\omega _0`$ and $`\mathrm{\Gamma }`$ due to the fermion-boson coupling are accounted for (to the extent that they do not induce a frequency dependence of $`\mathrm{\Gamma }`$) by considering $`\omega _0`$ and $`\mathrm{\Gamma }`$ to be effective temperature independent input parameters of our theory. The unperturbed (clean) fermionic Greens function is given by $`\widehat{G}_0^1(𝐤,i\omega _n)=\left[i\omega _n\tau _0ϵ_𝐤\tau _3\right]\sigma _0+\mathrm{\Delta }_𝐤\tau _2\sigma _2,`$ where $`\sigma _i`$ and $`\tau _i`$ are the Pauli-matrices in spin and Nambu-space, respectively. Changes in the fermionic Greens function due to the scattering of the localized mode are accounted for by summing the infinite series of diagrams shown in Fig. 1a which yields

$`\widehat{G}(r,r^{},\omega _n)`$ $`=`$ $`\widehat{G}_0(r,r^{},\omega _n)`$ (4)
$`+\widehat{G}_0(r,0,\omega _n)\widehat{S}(\omega _n)\widehat{G}_0(0,r^{},\omega _n).`$
Here

$$\widehat{S}(\omega _n)=\widehat{\mathrm{\Sigma }}(\omega _n)\left[\widehat{1}\widehat{G}_0(0,0,\omega _n)\widehat{\mathrm{\Sigma }}(\omega _n)\right]^1$$
(5)
and $`\widehat{\mathrm{\Sigma }}`$ is the local fermionic self-energy given by

$$\widehat{\mathrm{\Sigma }}(\omega _n)=T\underset{m}{}\tau _3\widehat{g}\widehat{G}_0(0,0,\omega _n\nu _m)D(\nu _m)\tau _3\widehat{g}.$$
(6)
with $`\widehat{g}`$ being determined by the spin dependence of $`g_\alpha `$. We show later that vertex corrections (see Fig. 1b) are negligible for the range of parameters considered here. Based on earlier studies Bal96 ; Bal02 , we expect that the inclusion of non-crossing diagrams in the calculation of $`\widehat{G}`$ does not qualitatively change our results and therefore reserve its discussion for a future publication Nyb03 . The DOS, $`N=A_{11}+A_{22}`$ with $`A_{ii}(𝐫,\omega )=2\mathrm{I}\mathrm{m}\widehat{G}_{ii}(𝐫,\omega +i\delta )`$ is obtained numerically from Eqs.(4)-(6) with $`\delta =0.2`$ meV.

As a frame of reference, we first consider a localized mode in an $`s`$-wave superconductor with normal state dispersion $`ϵ_k=k^2/2m\mu `$, $`m^1/\mathrm{\Delta }_0=15`$, Fermi momentum $`k_F=\pi /2`$, and a momentum independent SC gap $`\mathrm{\Delta }_0`$. For a spin-independent coupling, $`\widehat{g}=g\widehat{1}`$, no bound states are induced, similar to the case of a non-magnetic impurity in an $`s`$-wave superconductor Shiba68 . In contrast, if the coupling is spin-dependent the mode induces a bound state, with poles in the $`\widehat{S}`$-matrix leading to two peaks in the DOS. For simplicity we take $`g_{}=g`$, and $`g_{}=g`$, i.e., $`\widehat{g}=g\sigma _3`$, and in Fig. 2 present the DOS at $`𝐫=0`$ for $`\mathrm{\Gamma }=0^+`$ as a function of $`\omega _0`$ and temperature.

For $`T=0,\omega _0=0`$ (see Fig. 2a) the frequencies $`\pm \mathrm{\Omega }_b`$ of the bound state peaks are given by $`\mathrm{\Omega }_b/\mathrm{\Delta }_0=D_{}/D_+`$ where $`D_\pm =1\pm (gm/4)^2`$. This result is similar to the energy of a bound state induced by a static magnetic impurity Shiba68 . Note that increasing $`g`$, i.e., stronger scattering, leads to a decrease of $`\mathrm{\Omega }_b`$. For $`\omega _0,T0`$ the bound state frequencies are shifted from $`\mathrm{\Omega }_b`$, and in the limit $`T\mathrm{\Delta }_0`$ are determined by a transcendent equation. It is too cumbersome to be presented here, but can be simplified in the limit $`\mathrm{\Omega }_b,\omega _0\mathrm{\Delta }_0`$, where we find to leading order in $`\omega _0/\mathrm{\Delta }_0`$ and $`n_B(\omega _0)`$

$$\mathrm{\Omega }_b(T,\omega _0)=\mathrm{\Omega }_b+\omega _0/\pi n_B(\omega _0)\mathrm{\Delta }_0.$$
(7)
Note, that the localized mode effectively behaves as an “oscillating vector” with Ising symmetry (for a single bosonic degree of freedom) and frequency $`\omega _0`$. When $`\omega _0`$ increases (and in particular, when it becomes comparable to the rate of electronic scattering, $`\mu /\mathrm{}`$) subsequent scattering events begin to cancel each other; the effective scattering becomes weaker, and the peaks are shifted to higher energies. In particular, for $`\omega _0\mathrm{}`$ the bound state vanishes and the DOS becomes that of the unperturbed systems Nyb03 . This picture provides an explanation for our analytical results in Eq.(7) and fully agrees with the numerically obtained DOS shown in Fig. 2a for $`T=0`$ and several values of $`\omega _0`$. A non-zero $`\omega _0`$ also leads to a “dip” in the DOS at $`|\mathrm{\Omega }|=\mathrm{\Delta }_0+\omega _0`$, as shown in the inset of Fig. 2a. This dip is a characteristic signature of the localized mode, which distinguishes it from a static magnetic impurity. The origin of the dip lies in the poles of the retarded self-energy, $`\widehat{\mathrm{\Sigma }}_R(k,\mathrm{\Omega })`$ at $`|\mathrm{\Omega }|=E_k+\omega _0`$. The imaginary part of the retarded local self-energy, $`\widehat{\mathrm{\Sigma }}_{loc}`$, then possesses a threshold energy, $`\mathrm{\Omega }_c^+=\mathrm{\Delta }_0+\omega _0`$, with Im$`\widehat{\mathrm{\Sigma }}_{loc}1/\sqrt{(|\mathrm{\Omega }|\omega _0)^2\mathrm{\Delta }_0^2}`$ for $`|\mathrm{\Omega }|\mathrm{\Omega }_c^+`$, and Im$`\widehat{\mathrm{\Sigma }}_{loc}0`$ for $`|\mathrm{\Omega }|<\mathrm{\Omega }_c^+`$. As a result, Re$`\widehat{\mathrm{\Sigma }}_{loc}`$ exhibits a square-root singularity for $`|\mathrm{\Omega }|\mathrm{\Omega }_c^+0^+`$, which in the denominator of $`\widehat{S}`$ leads to a suppression, i.e., dip, in the DOS. This dip is therefore similar in nature to that observed in angle-resolved photoemission experiments in the HTSC arpes . Finally, we find that the amplitude of the bound state peaks exhibits spatial oscillations described by $`\mathrm{cos}(2k_Fr)e^{r/\xi }/r`$, where $`\xi =k_F/m\sqrt{\mathrm{\Delta }_0^2\mathrm{\Omega }_b^2}`$ Nyb03 .

For non-zero temperatures the $`\omega =+\omega _0`$ branch of the bosonic mode becomes populated, opening a new scattering channel and yielding significant changes in the DOS, as shown in Fig. 2b. We find that the specific temperature evolution of the bound state peaks depends on whether $`\mathrm{\Omega }_b`$ at $`T=0`$ is smaller or larger than $`\mathrm{\Omega }_c^{}=\mathrm{\Delta }_0\omega _0`$. This distinction arises since $`\widehat{\mathrm{\Sigma }}_R(k,\mathrm{\Omega })`$ now not only possesses poles at $`|\mathrm{\Omega }|=E_k+\omega _0`$, with weight $`1+n_B(\omega _0)`$ (for $`T,\omega _0\mathrm{\Delta }_0`$), but also poles at $`|\mathrm{\Omega }|=E_k\omega _0`$ with weight $`n_B(\omega _0)`$. Hence, Im$`\widehat{\mathrm{\Sigma }}_{loc}0`$ already for $`|\mathrm{\Omega }|\mathrm{\Omega }_c^{}`$, in contrast to the case at $`T=0`$. Thus, a bound state peak with $`\mathrm{\Omega }_c^{}<\mathrm{\Omega }_b<\mathrm{\Delta }_0`$ at $`T=0`$ becomes damped with increasing temperature, as shown in Fig. 2b for $`T=0.15\mathrm{\Delta }_0`$. At the same time, a new peak appears in the DOS at $`|\mathrm{\Omega }|=\mathrm{\Omega }_c^{}`$, since Im$`\widehat{\mathrm{\Sigma }}_{loc}`$ now also possesses square root singularities at $`\pm \mathrm{\Omega }_c^{}`$. With increasing temperature the damping of the bound state peak grows ($`T=0.4\mathrm{\Delta }_0`$) and its amplitude consequently decreases, while the newly emerged peak moves towards lower energies. In contrast, if $`\mathrm{\Omega }_b\mathrm{\Omega }_c^{}`$ (not shown) the bound state peak simply moves towards lower energies with increasing temperature, in agreement with our analytical results in Eq.(7). This downshift is expected since the population of the mode’s branches grows with increasing temperature, thus leading to stronger scattering. Thus, varying $`T`$ or $`\omega _0`$ has an opposite effect on the energies of the bound state peaks in the DOS.

We next consider a localized mode in a $`d_{x^2y^2}`$wave superconductor with SC gap $`\mathrm{\Delta }_𝐤=\mathrm{\Delta }_0(\mathrm{cos}k_x\mathrm{cos}k_y)/2`$ and $`\mathrm{\Delta }_0=25`$ meV. To simplify the discussion we consider a particle-hole symmetric dispersion, $`ϵ_k=2t(\mathrm{cos}k_x+\mathrm{cos}ky)`$, with $`t=300`$ meV. A band dispersion with strong particle-hole asymmetry, characteristic of the HTSC, yields qualitatively similar results Nyb03 . In contrast to the $`s`$-wave superconductor, a localized mode induces a resonance state even with a spin-independent coupling, $`g_\alpha =g`$. These resonances are a hallmark of the strong coupling theory, and are not observed in the weak scattering limit Bal02 . In Fig. 3 we present the DOS at $`𝐫=0`$ as a function of $`\omega _0`$ and temperature ($`\mathrm{\Gamma }=0^+`$). For $`T,\omega _0=0`$ (see Fig. 3a), the energy of the resonance peaks, $`\mathrm{\Omega }_r`$, is given by the solution of $`\mathrm{\Omega }_r=\pm 2\pi \mathrm{\Delta }_0/\left\{g\mathrm{log}\left[|\mathrm{\Omega }_r|/(4\mathrm{\Delta }_0)\right]\right\}`$, similar to the energies of resonance states induced by a non-magnetic impurity Bal96  (this result was also independently obtained by Si QSipc ). For $`|\mathrm{\Omega }|\omega _00`$,

$$\mathrm{Im}\mathrm{\Sigma }_{loc}(\mathrm{\Omega })=\frac{1}{\pi }\frac{|\mathrm{\Omega }|\omega _0}{\mathrm{\Delta }_0}K\left(\frac{|\mathrm{\Omega }|\omega _0}{\mathrm{\Delta }_0}\right)$$
(8)
where $`K`$ is the complete elliptic integral of the first kind, while Im$`\mathrm{\Sigma }_{loc}0`$ for $`|\mathrm{\Omega }|<\omega _0`$. Thus there exist a real gap in Im$`\mathrm{\Sigma }_{loc}`$, in contrast to Im$`G_{loc}`$. Performing a Kramers-Kronig transform, we find for $`\mathrm{\Omega }/\omega _01`$, Re$`\mathrm{\Sigma }_{loc}A\mathrm{\Omega }/\omega _0`$ where $`A=\frac{2}{\pi }_1^{\mathrm{}}𝑑x\mathrm{Im}\mathrm{\Sigma }_{loc}(x)/x^2`$. In the limit $`\mathrm{\Omega }_r\mathrm{\Delta }_0`$, the resulting equation for $`\mathrm{\Omega }_r`$ is

$$\mathrm{\Omega }_r^2=\frac{4\pi \mathrm{\Delta }_0\mathrm{\Omega }}{g^2}\left[A\mathrm{log}\left(\frac{4\mathrm{\Delta }_0}{|\mathrm{\Omega }_r|}\right)\right]^1$$
(9)
As expected, increasing $`\omega _0`$ leads to a shift of $`\mathrm{\Omega }_r`$ to higher energies, in agreement with our numerical results shown in Fig. 3a. Moreover, for $`\omega _0\mathrm{}`$ we again recover the unperturbed DOS (see, e.g., $`\omega _0=200`$ meV). The logarithmic divergence in $`\mathrm{\Sigma }_{loc}`$ at $`|\mathrm{\Omega }|=\mathrm{\Delta }_0+\omega _0`$ leads to a dip in the DOS (see arrows in the inset of Fig. 3a). A much weaker dip emerges in the weak scattering limit Bal02 .

The specific temperature evolution of the resonance peaks again depends on the relative order of $`|\mathrm{\Omega }_r|`$ at $`T=0`$ and $`\mathrm{\Omega }_c^{}=\mathrm{\Delta }\omega _0`$, similar to the $`s`$-wave case. For $`|\mathrm{\Omega }_r|<\mathrm{\Omega }_c^{}`$ (see Fig. 3b) the resonance peak is shifted to lower energies with increasing temperature. In addition, the opening of a second scattering channel for $`T0`$ leads to a logarithmic divergence in $`\mathrm{\Sigma }_{loc}`$, and consequently a dip in the DOS at $`|\mathrm{\Omega }|=\mathrm{\Omega }_c^{}=15`$ meV, as indicated by the arrows. In contrast, for $`\mathrm{\Omega }_c^{}<|\mathrm{\Omega }_r|<\mathrm{\Delta }_0`$ (see Fig. 3c), the same opening of a second scattering channel leads to an increase in Im$`\mathrm{\Sigma }_{loc}`$ for $`|\mathrm{\Omega }|>\omega _0`$, which in turn shifts the peaks to higher energies and suppresses their amplitude. The smaller peaks at $`|\mathrm{\Omega }|=\mathrm{\Omega }_c^{}=5`$ meV (see arrows) arise from the opening of the second scattering channel, and the resulting logarithmic divergence of Im$`\mathrm{\Sigma }_{loc}`$ at $`\mathrm{\Omega }_c^{}`$. Finally, we find that due to the momentum dependence of the SC gap, the spatial dependence of the resonance peaks (not shown) Nyb03  is qualitatively similar to that observed experimentally near impurities in the HTSC exp .

To study the effects of a mode with finite lifetime, $`\mathrm{\Gamma }0`$, we present in Fig. 4, the DOS in a $`d_{x^2y^2}`$-wave and $`s`$-wave (see inset) superconductor at $`T=0`$ and for several values of $`\mathrm{\Gamma }`$.

In both superconductors the peaks shift to higher energies with increasing $`\mathrm{\Gamma }`$. Note that a non-zero $`\mathrm{\Gamma }`$ does not directly change the lifetime of the induced state, i.e., the width of the peaks. This is particularly evident for the $`s`$-wave case (see inset). A non-zero $`\mathrm{\Gamma }`$ lowers the threshold frequency for Im$`\widehat{\mathrm{\Sigma }}_{loc}0`$ from $`\mathrm{\Omega }_c^+=\mathrm{\Delta }_0+\omega _0`$ (for $`\mathrm{\Gamma }=0^+`$) to $`\mathrm{\Omega }_c=\mathrm{\Delta }_0`$. However, for $`|\mathrm{\Omega }|<\mathrm{\Delta }_0`$, we still have Im$`\widehat{\mathrm{\Sigma }}_{loc}=0`$ and the width of the bound or resonance peaks remains unaffected by $`\mathrm{\Gamma }0`$.

Finally, we consider the lowest order vertex correction, shown in Fig. 1b. For a quadratic normal state dispersion, the general form of this correction is $`\delta g/g=\left(\frac{gm}{4}\right)^2F(\mathrm{\Omega }_m,\nu _n,\omega _0,\mathrm{\Delta }_0)`$ where $`\mathrm{\Omega }_m,\nu _n`$ are the bosonic and fermionic Matsubara frequencies, respectively, and $`F`$ is a bounded function of $`O(1)`$. In particular, for an $`s`$-wave superconductor, we have $`F=+1`$ for $`\mathrm{\Omega }_m,\nu _n,\omega _0\mathrm{\Delta }_0`$, and $`F=1`$ for $`\mathrm{\Omega }_m,\nu _n,\omega _0\mathrm{\Delta }_0`$, implying that the bare vertex is enhanced (suppressed) for frequencies smaller (larger) than the SC gap. In contrast, in a $`d_{x^2y^2}`$-wave superconductor we obtain $`F0`$ for $`\mathrm{\Omega }_m,\nu _n,\omega _0\mathrm{\Delta }_0`$, and vertex corrections become irrelevant. Thus, vertex corrections can in general be neglected for $`gm/4<1`$, which applies to all cases considered above.

In conclusion, we consider the effects of a localized bosonic mode on the electronic structure in $`s`$\- and $`d_{x^2y^2}`$-wave superconductors. The mode acts as a new type of “defect”, leading to the emergence of bound or resonance states. We identify several characteristic features, such as a “dip” in the DOS and the temperature dependence of the mode induced peaks, that qualitatively distinguishes the effects of a localized mode from those of static impurities. This result provides further insight into the nature of impurities observed in STM experiments.

We would like to thank A. Abanov, A. Balatsky, J.C. Davis, A. de Lozanne and Q. Si for stimulating discussions.


format: text

text:
# Pseudopotential model of ultracold atomic collisions in quasi-one- and two-dimensional traps

## I Introduction

Unprecedented control has been gained over atomic collisions through cooling to nanoKelvin temperatures and selection of internal hyperfine states. Further control of these systems is now being gained by manipulating the atoms with external fields, including magnetic fields and optical lattices. An opportune example where both collisions and the external trapping potential are essential is the superfluid-Mott insulator transition in an optical lattice. This transition was recently demonstrated by beginning with an atomic Bose-Einstein condensate and adiabatically turning on an optical lattice Greiner02 . The final Mott insulator state has a fixed number of atoms per lattice site. Such a system represents an ideal ensemble for measuring scattering properties of the atoms. The free-space scattering amplitudes determine the energy eigenvalues and loss rates of the system and vice versa.

There are many current and proposed applications of cold atoms in optical lattices. One is high-precision measurement of atomic potentials, by determining positions of weakly-bound molecular states for example Bolda02a . Another is the implementation of quantum logic gates with neutral atoms. Several proposed implementations combine optical lattices with internal-state-dependent cold collisions Brennen99 ; Jaksch99 . One can also consider the combination of an optical lattice along one or two directions and a weak dipole trap in the remaining directions. Such highly anisotropic trapping configurations were already used in experiments on Bose-Einstein condensate number-squeezing Orzel01  and controlled loading HeckerDenschlag02 , as well as looking for dimensional effects on the condensate density Goerlitz01 . Optical Bongs01  or magnetic Thywissen99  waveguides have been used for guiding cold atoms. These quasi-one- and two- dimensional configurations can be combined with a tunable scattering length, such as from magnetic Feshbach resonance, to lead to new physical regimes. Workers on Bose-Einstein condensate experiments are observing increased phase-fluctuation in the crossover to one dimension Dettmer01 ; Schreck01 . Theoretically, the (one-dimensional) Tonks gas regime can occur at low density when bosons are tightly confined along two directions and weakly confined along the third Olshanii98 ; Dunjko01 . One-half anyon statistics and the fractional quantum Hall effect may be possible for bosons confined tightly along one direction and weakly along the other two Paredes01 . Another proposal is the Kosterlitz-Thouless transition for interacting bosons in two-dimensions Stoof93 . For all these applications, one first needs to understand the relation between cold collisions in free space and in a trap with arbitrary geometry. Then the many-body physics can be treated on the basis of effective low-dimensional interactions. In this article, we show how to compute the eigenvalues from the free-space scattering data and trap frequencies for arbitrary axisymmetric harmonic traps, emphasizing the one- and two-dimensional regimes. We also describe a novel bound state induced by the trap in both low-dimensional regimes when the scattering length is large in magnitude.

Previous theoretical work has included exact solutions of collisions for special interatomic potentials in an isotropic harmonic trap Busch98a ; Block02 , and comparison with results including realistic ground-state interatomic potentials. We have shown how an effective-scattering-length combined with a Fermi pseudopotential can be used to calculate the energy eigenvalues for collisions in an isotropic harmonic trap Blume02 ; Bolda02a . We compared our model with full interaction potential results for both a single-channel collision and a multi-channel collision with a magnetically tunable Feshbach resonance. In both cases, the model can accurately treat tight traps, as long as the trap size is larger than the van der Waals scale length. On the strength of this evidence we propose to also apply the effective-scattering length model to treat ground state collisions in axisymmetric harmonic traps.

The effect of confinement along only one or two dimensions has also been considered theoretically by others. For one dimension, the scattering could be represented by a one-dimensional delta function provided the coupling strength was renormalized by the confining trap Olshanii98 . Quasi-two-dimensional scattering solutions were found for two dimensions, where the renormalization is of the two-dimensional scattering length Petrov00b . We will consider the applicability of these two results in the case of very prolate and oblate traps respectively.

The outline of our paper is as follows. In the next section, we state the problem of atoms colliding in an axisymmetric trap and the regime of applicability of our effective-scattering-length model. In Sec. III we show how this problem can be solved in a particular basis. Sec. IV contains the results of the calculation for the one-dimensional regime and a comparison with a one-dimensional scattering theory; Sec. V contains the same but for the two-dimensional regime. In Sec. VI we discuss a novel trap-induced bound state which appears in both the one- and two- dimensional regimes for large positive or negative scattering length. We give an example of magnetically tunable Feshbach resonance for Na atoms in quasi-one-dimension in Sec. VII. We conclude the paper in Sec. VIII. In the appendix we derive the matrix elements of the anisotropic potential required for the calculation.

## II Two atoms colliding in an anisotropic harmonic trap

We assume ultracold atoms are trapped in an optical lattice detuned far off resonance. Specific optical lattice potentials for different angular and polarization configurations are calculated in Petsas94 . For our purposes we assume that two atoms in specific internal sublevels remain near a local minimum of the potential. With the assumption of local azimuthal symmetry about an axis through a potential minimum, we approximate the anisotropic potential near a particular site by

$`V_t(𝐫_j)`$ $`=`$ $`{\displaystyle \frac{1}{2}}m\left[\omega _{}^2(x_j^2+y_j^2)+\omega _z^2z_j^2\right],`$ (1)
where $`𝐫_j`$ is the positions of atom $`j=1`$ or $`2`$, $`m`$ is the atomic mass, and $`\omega _{}`$ and $`\omega _z`$ are trapping frequencies. We define the trap anisotropy

$$A=\frac{\omega _z}{\omega _{}}$$
(2)
so that one-dimensional physics is approached for the oblate trap $`A1`$ (”cigar”) and two-dimensional physics for the prolate trap $`A1`$(”pancake”). The length scales associated with the transverse and longitudinal trap directions are

$$d_{}=\sqrt{\frac{\mathrm{}}{\mu \omega _{}}},d_z=\sqrt{\frac{\mathrm{}}{\mu \omega _z}},$$
(3)
where $`\mu =m/2`$ is the reduced mass of the atom pair.

As in the case of an isotropic harmonic trap, the center of mass and relative motion are separable. The center of mass motion is independent of the interatomic potential, $`V_{\mathrm{int}}(r)`$, and has the usual harmonic oscillator solutions. The relative coordinate Hamiltonian is given in spherical coordinates by

$$\widehat{H}=\frac{\mathrm{}^2}{2\mu r^2}\frac{}{r}\left(r^2\frac{}{r}\right)+\frac{\widehat{L}^2}{2\mu r^2}+V(r,\theta )+V_{\mathrm{int}}(r)$$
(4)
where $`𝐫=𝐫_1𝐫_2`$ and $`r=|𝐫|`$. The interatomic orbital angular momentum operator$`\widehat{L}`$ gives the partial wave quantum numbers $`l=`$ 0, 1, 2, … for $`s`$-, $`p`$-, $`d`$-, … waves. The potential due to the trap written in terms of the spherical harmonic $`Y_{20}`$ is

$$V(r,\theta )=\frac{1}{2}\mu \omega ^2r^2\left[1+\sqrt{\frac{16\pi }{5}}\mathrm{\Lambda }Y_{20}(\theta ,0)\right],$$
(5)
where the mean-square trap frequency appears as

$$\omega =\sqrt{\frac{2\omega _{}^2+\omega _z^2}{3}}$$
(6)
and

$$\mathrm{\Lambda }=\frac{\omega _z^2\omega _{}^2}{\omega _z^2+2\omega _{}^2}=\frac{A^21}{A^2+2}.$$
(7)
The term proportional to $`\mathrm{\Lambda }`$ in Eq. (5) defines the anisotropic part $`\widehat{H}^{(1)}`$ of the Hamiltonian $`\widehat{H}=\widehat{H}^{(0)}+\widehat{H}^{(1)}`$. The length scale associated with the mean-square frequency is defined to be

$$d=\sqrt{\frac{\mathrm{}}{\mu \omega }}.$$
(8)
The interatomic $`V_{\mathrm{int}}(r)`$ potential for two ground state atoms approaches $`C_6/r^6`$ at large internuclear separation $`r`$. The associated van der Waals length scale is $`x_0=(2\mu C_6/\mathrm{}^2)^{1/4}/2`$ Gribakin93 ; Weiner99 ; Williams99 . It gives the approximate size of the potential, that is, the wavefunction takes on its asymptotic scattering form for $`rx_0`$. We have previously approximated the exact Born-Oppenheimer potential by the energy-dependent Fermi pseudopotential Bolda02a ; Blume02 ; Huang57

$$\widehat{V}_{\mathrm{eff}}(r;E)=\frac{4\pi \mathrm{}^2a_{\mathrm{eff}}(E)}{m}\delta (𝐫)\frac{}{r}r.$$
(9)
The dependence on collision energy $`E`$ is due to the energy-dependent scattering length, defined as

$$a_{\mathrm{eff}}(E)=\frac{\mathrm{tan}\delta _0(E)}{k},$$
(10)
where $`E=\mathrm{}^2k^2/2\mu `$ and $`\delta _0(E)`$ is the $`s`$-wave collisional phase shift. $`S`$-wave scattering predominates for ultracold collisions except for the case of identical fermions.

The pseudopotential approximation is valid provided the van der Waals length scale $`x_0`$ is less than the smallest harmonic oscillator width, $`x_0\text{min}\{d_{},d_z\}`$ Bolda02a . We assume that the energy shifts due to higher partial waves are negligible in comparision to those of the $`s`$-wave. This tends to be true for ultracold collisions because the centrifugal barrier heights are large compared to the collision energy. We also neglect inelastic losses, that is, the imaginary part of the scattering length is much smaller than the real part Tiesinga00a .

Since spherical symmetry is broken by the anisotropic potential, $`\widehat{L}^2`$ does not commute with the Hamiltonian. Consequently, partial waves with the same parity with respect to $`l`$ are coupled. While the projection $`\widehat{L}_z`$ of angular momentum on the $`z`$-axis does commute with the Hamiltonian, only its $`m_l=0`$ eigenstates are affected by $`s`$-wave scattering. Thus we only compute the energies of even partial wave, $`m_l=0`$ states. We do not consider odd partial waves, should they be present for distinguishable bosons, because they have negligible energy shifts in the limit of very low collision energy.

## III Method of eigenvalue solution and self-consistent energies

We need to solve the eigenvalue problem for the Hamiltonian Eq. (4) self-consistently, because of the energy-dependent scattering length in the pseudopotential term. As in Ref. Bolda02a , this is done in two steps. We first obtain the eigenvalues $`E_i(a/d,A)/\mathrm{}\omega `$ in scaled trap energy units for fixed values of $`A`$ and the energy-independent scaled scattering length $`a/d`$. The self-consistent energy eigenvalues for an actual system with $`a_{\mathrm{eff}}(E)`$ from Eq. (10) are then found graphically, for a given $`A`$, by superposing a plot of $`E_i(a/d,A)`$ as a function of $`a/d`$ and a plot of $`a_{\mathrm{eff}}(E)/d`$, with $`E`$ as the ordinate and $`a_{\mathrm{eff}}/d`$ as the abscissa (cf. Fig. 1 and Fig 2 of Bolda02a ). The points where the curves intersect determine the self-consistent energies.

In the rest of this section we focus on obtaining $`E_i(a/d,A)`$, since $`a_{\mathrm{eff}}(E)`$ can be found from a standard free-space scattering calculation. We use the partial wave expansion of the wavefunction,and expand each partial wave, with the exception of $`s`$-waves, in isotropic harmonic oscillator eigenfunctions of frequency $`\omega `$. For the $`s`$-wave part, we use the analytic eigenfunctions of the isotropic harmonic oscillator with a Fermi pseudopotential proportional to scattering length $`a`$ Busch98a . These automatically incorporate the singular nature of the wavefunction at the origin. Since $`m_l=0`$ we can set the spherical coordinate $`\varphi =0`$, and consequently we write

$`\psi (r,\theta ;a)`$ $`=`$ $`{\displaystyle \underset{n}{}}c_{n0}Q_n(r;a)Y_{00}(\theta ,0)`$ (11)
$`+{\displaystyle \underset{n}{}}{\displaystyle \underset{l>0}{}}c_{nl}R_{nl}(r)Y_{l0}(\theta ,0)`$
$$Q_n(r;a)=\frac{2a}{\sqrt{\pi }d^2}\sqrt{\frac{\nu _n}{a}}\mathrm{\Gamma }(\nu _n)U(\nu _n,\frac{3}{2};\frac{r^2}{d^2})e^{\frac{r^2}{2d^2}}$$
(12)
$$R_{nl}(r)=\left[\frac{2(n!)}{\mathrm{\Gamma }(n+l+\frac{3}{2})}\right]^{\frac{1}{2}}\left(\frac{r}{d}\right)^{\frac{l}{2}+\frac{3}{4}}L_n^{(l+\frac{1}{2})}\left(\frac{r^2}{d^2}\right)e^{\frac{r^2}{2d^2}}$$
(13)
where the summations are over all nonnegative integers $`n`$ and even $`l>0`$. Here $`U`$ is the (second) Kummer confluent hypergeometric function, $`L_n^{(l+\frac{1}{2})}`$ are the Laguerre polynomials, and $`\mathrm{\Gamma }`$ is the Gamma function Abramowitz72 . The nonintegral $`s`$-wave quantum numbers $`\nu _n`$ for an isotropic trap are determined by

$$\frac{a}{d}=\frac{1}{2}\mathrm{tan}\pi \nu _n\frac{\mathrm{\Gamma }(\nu _n+1)}{\mathrm{\Gamma }(\nu _n+\frac{3}{2})}.$$
(14)
This equation also is used in the calculation of the derivative in Eq. (12).

We need the Hamiltonian matrix elements in the basis used in Eq. (11). The isotropic part of the Hamiltonian Eq. (4) gives only diagonal matrix elements:

$$H_{nl;nl}^{(0)}=\left(2n+l+\frac{3}{2}\right)\mathrm{}\omega $$
(15)
for $`l>0`$. Only the $`s`$-wave diagonal matrix elements,

$$H_{n0;n0}^{(0)}=\left(2\nu _n+\frac{3}{2}\right)\mathrm{}\omega ,$$
(16)
are affected at low energy by atom-atom interactions proportional to the scattering length. The anisotropic part of the Hamiltonian, Eq. (4) contributes both diagonal and off-diagonal matrix elements:

$$H_{nl;n^{}l^{}}^{(1)}=\sqrt{\frac{4\pi }{5}}\mu \omega ^2\mathrm{\Lambda }n,l|r^2Y_{20}(\theta ,0)|n^{},l^{}$$
(17)
for all principal quantum numbers $`n,n^{}`$ and even partial waves $`l,l^{}`$. The derivation of these matrix elements is given in the Appendix. The diagonalization of the Hamiltonian matrix is straightforward with sparse matrix eigenvalue routines. For the most extreme anisotropies considered in this paper (as small as $`0.01`$ and as large as $`100`$), we required a maximum $`l=600`$ and a maximum $`n=600`$ to compute the lowest few eigenvalues for all values of scattering length. We have checked that the correct solutions are approached as $`A1`$.

## IV Quasi-one-dimensional trap

We have found solutions in the quasi-one-dimensional regime $`A1`$. Figure 1 shows $`E_i(a/d,A=0.01)`$ versus $`a/d`$. The figure also shows the eigenvalues $`E_i^{1D}(a/d)`$ for a purely one-dimensional model, corresponding to interaction via a delta function in $`z`$ and trapping along $`z`$ only. Following Ref. Busch98a , the eigenvalues are

$$E_i^{1D}=\left(\frac{3}{2}\right)^{\frac{1}{2}}(1+\frac{A^2}{2})^{\frac{1}{2}}\left[1+A\left(2\nu _i^{1D}+\frac{1}{2}\right)\right]\mathrm{}\omega $$
(18)
where the $`\nu _i^{1D}`$ satisfy

$$\mathrm{tan}\pi \nu _i^{1D}\frac{\mathrm{\Gamma }(\nu _i^{1D}+1)}{\mathrm{\Gamma }(\nu _i^{1D}+\frac{1}{2})}=g^{1D}.$$
(19)
The one-dimensional interaction, $`\mathrm{}\omega g^{1\mathrm{D}}\delta (z/d)`$, is related to the three-dimensional scattering length through Olshanii98

$$g^{1D}=\frac{\left(\frac{3}{2}\right)^{\frac{1}{4}}A^{\frac{1}{2}}(1+\frac{A^2}{2})^{\frac{1}{4}}\frac{a}{d}}{11.4603\left(\frac{3}{2}\right)^{\frac{1}{4}}(1+\frac{A^2}{2})^{\frac{1}{4}}\frac{a}{d}}.$$
(20)
Our eigenvalues agree well with those from the one-dimensional model for $`E>0.8\mathrm{}\omega `$. To get this agreement, it is crucial to include the renormalization in the denominator of Eq. (20). The levels for $`E>1.23\mathrm{}\omega `$ represent trap levels aligned along the weak trapping direction with spacing about $`2\mathrm{}\omega _z=0.024\mathrm{}\omega `$. In Fig. 1, there is no difference within our numerical accuracy between our lowest three trap levels and those from the one-dimensional model maximum energy difference, while the difference for the highest trap level shown is $`0.002\mathrm{}\omega `$. The eigenvalues for large positive or negative scatttering length approach the same asymptotic values. The lowest energy state is not predicted accurately by the one-dimensional model; we discuss this special state further in Sec. VI.

The wavefunction $`r\psi (x=r\mathrm{sin}\theta ,y=0,z=r\mathrm{cos}\theta )`$ , corresponding to the second lowest energy, with $`A=0.01`$ and $`a/d=25`$, is plotted in Fig. 2. In Fig. 2(a) the one dimensional nature of the wavefunction is apparent on a scale large compared with $`d`$. The variation in $`x`$ is approximately Gaussian, while scattering results in the dip along the line $`z=0`$. A close-up of the wavefunction in Fig. 2(b) reveals how the scattering crosses over to a isotropic three-dimensional character at short interatomic distance. The function $`r\psi `$ is finite at the origin as a consequence of the pseudopotential scattering.

We have compared some of our results to a recent diffusion quantum Monte-Carlo study on the ground state of interacting bosons in elongated traps Blume02b . Our lowest trap state agrees with the ground state of that method to within one percent for anisotropies in the range $`1A0.01`$ for a fixed positive value of the scattering length, even though that study assumed a hard-core potential of size $`a`$.

## V Quasi-two-dimensional trap

We have also computed solutions in the quasi-two-dimensional regime, $`A1`$. Figure 3 shows $`E_i(a/d,A=100)`$ versus $`a/d`$. The figure also shows the eigenvalues $`E_i^{2D}(a/d)`$ for a purely two-dimensional model, corresponding to interaction via zero-range two-dimensional scattering and trapping in the $`x,y`$-plane only. Again following Busch98a , the eigenvalues of a two-dimensional trap are

$$E_i^{2D}=\left(\frac{3}{2}\right)^{\frac{1}{2}}(1+\frac{A^2}{2})^{\frac{1}{2}}\left[\frac{A}{2}+2\nu _i^{2D}+1\right]\mathrm{}\omega ,$$
(21)
where the $`\nu _i^{2D}`$ satisfy

$$\mathrm{\Psi }(\nu _i^{2D})=\frac{1}{g^{2D}}$$
(22)
and $`\mathrm{\Psi }`$ is the digamma function Abramowitz72 . The two-dimensional scattering is mediated by an interaction strength related to the three-dimensional scattering length through

$$g^{2D}=\frac{\left(\frac{3}{2\pi ^2}\right)^{\frac{1}{4}}A^{\frac{1}{2}}\left(1+\frac{A^2}{2}\right)^{\frac{1}{4}}\frac{a}{d}}{1+\left(\frac{3}{2\pi ^2}\right)^{\frac{1}{4}}\mathrm{ln}\left(\frac{0.915A}{4\pi }\right)A^{\frac{1}{2}}\left(1+\frac{A^2}{2}\right)^{\frac{1}{4}}\frac{a}{d}}.$$
(23)
This expression for $`g^{2D}`$ is derived by simple algebra from equations in Ref. Petrov00b .

Our eigenvalues agree well with this model for $`E>0.8\mathrm{}\omega `$. To obtain this agreement it is crucial to include the renormalization in the denominator of Eq. (23). The trap levels above $`E0.9\mathrm{}\omega `$ are spaced by about $`2\mathrm{}\omega _{}=0.03462\mathrm{}\omega _{}`$. In Fig. 3, the difference between our lowest trap level and that from the two-dimensional model is $`0.0015\mathrm{}\omega `$ at $`a/d=\pm 20`$; the difference for the highest trap level shown is $`0.002\mathrm{}\omega `$ at the same values of scattering length. The eigenvalues for large positive or negative scattering length approach the same asymptotic values. Again we note the exception that the lowest energy state is not predicted well by the renormalized two-dimensional model (see Sec. VI). A feature of the two-dimensional physics seen in Fig. 3 is that for $`|a|d`$, all eigenvalue curves except the lowest have nearly the same slope.

A sample wavefunction $`r\psi (x=r\mathrm{sin}\theta ,y=0,z=r\mathrm{cos}\theta )`$ , corresponding to the second lowest energy, at $`A=100`$ and $`a/d=25`$, is plotted in Fig. 4. (Note that the $`x`$ and $`z`$ axes are interchanged from Fig. 2.)

The two-dimensional nature of the physics is immediately apparent on a scale large compared with $`d`$. The variation in $`z`$ is approximately Gaussian, with a dip along $`x=0`$. The close-up Fig. 4(b) reveals how the scattering crosses over to a isotropic three-dimensional character at short interatomic distance, as in the quasi-one-dimensional case.

## VI Trap-induced bound state

Curiously, for large magnitudes of the scattering length, the lowest energy eigenvalue appears near $`\mathrm{}\omega /2`$ for all values of trap anisotropy. This behavior can be explained by examining the wavefunctions of these states plotted in Figs. 5 and 6, at $`a=25`$, for $`A=0.01`$ and $`100`$ respectively.

Both wavefunctions are more nearly isotropic than those corresponding to higher levels at the same parameters (see Figs. 2 and 4), and their extent is roughly the mean trap length scale $`d`$. Similar wavefunctions are obtained for large positive scattering length.

The energy of the lowest state is plotted in Fig. 7 for $`a/d=\pm 100`$ as a function of anisotropy.

The shape of this curve can be understood from a perturbative picture in the anisotropic interaction $`H^{(1)}`$, defined in Eq. (7). To zeroth-order in this picture the state is the lowest level of the isotropic trap of frequency $`\omega `$ and a pseudopotential with scattering length $`a`$. (Recall that for an isotropic trap with $`a\mathrm{}`$ the lowest energy is $`\mathrm{}\omega /2`$.) There is no first-order correction in $`H^{(1)}`$, while to second-order the energy gets a small correction proportional to $`\mathrm{\Lambda }^2`$. Indeed, we find that a quadratic fit in $`\mathrm{\Lambda }`$ at $`a/d=100`$ results in $`E_0/\mathrm{}\omega =0.50540.0562\mathrm{\Lambda }^2`$; compare with the zeroth-order energy $`0.5056\mathrm{}\omega `$ for the lowest state in an isotropic trap at the same scattering length. Similarly, at $`a/d=100`$ we find $`E_0/\mathrm{}\omega =0.49410.0549\mathrm{\Lambda }^2`$ compared with the zeroth-order energy $`0.4943\mathrm{}\omega `$ in the isotropic case.

When neglecting the pseudopotential, the energy of the lowest state is $`E(a=0)/\mathrm{}\omega =(1+A/2)/\sqrt{2/3+A^2/3}>\sqrt{3}/2`$, so we should properly consider a state to be bound when its energy is lower than this. Thus we denote such a state the trap-induced bound state. One can also think of it as an artificial molecule with an extent given roughly by the size of the tightest trap dimension. As the scattering length approaches zero from the positive side, we recover the usual molecular bound state with an energy below zero.

## VII Feshbach resonance in quasi-one-dimensional trap

One possible way of varying the atomic interaction strengths experimentally is through the use of a tunable Feshbach resonance state. Consequently we describe a quasi-one-dimensional magnetically tunable Feshbach resonance by using the self-consistent energy method with the eigenvalues of Fig. 1. We consider two Na atoms in their lowest hyperfine levels, for which an $`s`$-wave Feshbach resonance occurs near 90.9 mT Inouye98 ; vanAbeelen99 ; Mies00 . The scattering length is highly energy- and magnetic-field- dependent near the resonance. We use the effective scattering length from a close-coupling calculation, as described in Ref. Bolda02a ; Mies00 . The trap frequencies are taken to be $`\omega _{}/2\pi =612`$ kHz, $`\omega _z/2\pi =6.12`$ kHz, so the trap anisotropy is $`A=0.01`$ and $`\omega /2\pi =500`$ kHz. Using the procedure outlined at the beginning of Sec. III, we predict the eigenvalues as a function of applied magnetic field near the resonance in Fig. 8. As the magnetic field is tuned through resonance, the lowest state goes continuously from a molecular state with $`E<0`$, to the trap-induced bound state $`E\mathrm{}\omega /2`$, to the lowest quasi-one-dimensional trap state at $`E\sqrt{3/2}\mathrm{}\omega 1.23\mathrm{}\omega `$. For the trap frequencies chosen here one can change from a molecular bound state to a trap state by tuning the magnetic field 0.01 mT. The trap states are smoothly shifted up by $`2\mathrm{}\omega _z`$ as the magnetic field is increased.

## VIII Conclusion

We have argued that an energy-dependent pseudopotential approach may be used to calculate eigenvalues of two ultracold atoms colliding in an axisymmetric harmonic trap. Furthermore, we have numerically solved for the eigenvalues of the axisymmetric trap with an $`s`$-wave pseudopotential interaction proportional to a scattering length. These results can be considered a generalization of the isotropic trap case previously solved Busch98a ; Bolda02a . Our results show that one- and two- dimensional interaction regimes can be reached, but that the interactions become renormalized by the tight trapping potential when the magnitude of the effective scattering length is large compared with a mean trap length. Remarkably, in the case of scattering length of large magnitude, we find a nearly isotropic state with energy near $`\mathrm{}\omega /2`$ for all values of trap anisotropy. This is a trap-induced bound state. The size of state is controllable by the tightest trap frequency. We show by an example that this state can be reached with the current techniques of magnetically tunable Feshbach resonance applied to atoms in an optical lattice.

The numerical techniques used in this article may also be useful when the interatomic interaction becomes sufficiently long-range compared with the tightest trap direction or is anisotropic such that the $`s`$-wave pseudopotential approach becomes insufficient. This is of importance, for example, for dipole-dipole interactions Brennen99 ; Goral02 ; Derev03 . The possibility of manipulating shape resonances (such as the $`d`$-wave resonance in Na collisions) with the trap should not be overlooked. The effective-scattering-length model and partial wave expansion in our numerical technique should also solve the eigenvalue problem for atoms colliding in separated traps, as can occur in a state-dependent optical lattice. This is particularly important for the proposals on quantum computing with neutral atoms. We are continuing work in this area.

## Acknowledgments

Discussions with M. Olshanii and B. Gao were helpful. We thank D. Blume for making available her results on two-body ground state energies. ELB was supported from the National Research Council. ET and PSJ acknowledge support from the Office of Naval Research.

## Appendix: Matrix elements of anisotropic potential

In this appendix we evaluate the matrix elements of the anisotropic potential term $`H_{nl;n^{}l^{}}^{(1)}`$ in the basis of partial waves, isotropic harmonic oscillator functions and the irregular $`s`$-wave oscillator eigenfunctions.

The partial wave expansion applied to the angle-dependent factor of the Hamiltonian is evaluated with the three-spherical-harmonic formula,

$`I_{ll^{}}`$ $`=`$ $`\sqrt{{\displaystyle \frac{4\pi }{5}}}{\displaystyle 𝑑\mathrm{\Omega }Y_{l^{}0}^{}(\mathrm{\Omega })Y_{20}(\mathrm{\Omega })Y_{l0}(\mathrm{\Omega })}`$ (24)
$`=`$ $`\sqrt{{\displaystyle \frac{2l+1}{2l^{}+1}}}2l;00|2l;l^{}0^2.`$
Evaluation of the Clebsch-Gordan coefficient $`2l;00|2l;l^{}0`$ shows that

$`I_{ll}`$ $`=`$ $`{\displaystyle \frac{l(l+1)}{(2l3)(2l+3)}}`$ (25)
$`I_{l,l+2}`$ $`=`$ $`I_{l+2,l}={\displaystyle \frac{3(l+1)(l+2)}{2(2l+3)\sqrt{(2l+1)(2l+5)}}}`$ (26)
while all other angular matrix elements are zero.

This leaves the radial factor of the matrix element to be computed. Throughout the remainder of the appendix kets refer to the radial part of the basis functions only, so that for (cf. Eqs. (12) and (13))

$`r|nl`$ $`=`$ $`R_{nl}(r)`$ (27)
$`r|n0`$ $`=`$ $`Q_n(r;a).`$ (28)
In this notation

$$H_{nl;n^{}l^{}}^{(1)}=\mu \omega ^2\mathrm{\Lambda }I_{ll^{}}nl|r^2|n^{}l^{}.$$
(29)
The analytic evaluation of the radial matrix elements $`nl|r^2|n^{}l`$ is most conveniently carried out using $`n`$\- and $`l`$-ladder operators. From the factorization method applied to the radial Schrödinger equation Infeld51 , the $`n`$-raising and $`n`$-lowering operators are

$$\widehat{b}_{nl}^\pm =\pm r\frac{}{r}\pm \frac{1}{2}r^2+2n+l+\frac{1}{2}$$
(30)
$`\widehat{b}_{nl}^{}|nl`$ $`=`$ $`\sqrt{2n(2n+2l+1)}|n1,l`$ (31)
$`\widehat{b}_{nl}^+|n1,l`$ $`=`$ $`\sqrt{2n(2n+2l+1)}|nl,`$ (32)
while the $`l`$-raising operator is

$$\widehat{}_l^+=\frac{}{r}+r\frac{l}{r}$$
(33)
such that

$$\widehat{}_l^+|nl=2\sqrt{n}|n1,l+1.$$
(34)
(Note that these operators are applied to the normalized radial eigenfunctions.)

For $`l=l^{}>0`$, we use the fact that

$$r^2=2n+l+\frac{3}{2}\frac{1}{2}\widehat{b}_{n+1,l}^+\frac{1}{2}\widehat{b}_{n,l}^{}$$
(35)
and orthonormality to find

$`nl|r^2|nl`$ $`=`$ $`2n+l+{\displaystyle \frac{3}{2}}`$ (36)
$`n+1,l|r^2|nl`$ $`=`$ $`nl|r^2|n+1,l`$ (37)
$`=`$ $`{\displaystyle \frac{1}{2}}\sqrt{2(n+1)(2n+2l+3)},`$
with all other equal-$`l`$ matrix elements vanishing. (We recognize Eq. (36) as a consequence of the quantum virial theorem for the isotropic harmonic oscillator.)

For $`l^{}=l+2`$ (but $`l0`$), we begin with

$$nl|r^2|n^{},l+2=\frac{nl|r^2\widehat{}_{l+1}^+\widehat{}_l^+|n^{}+2,l}{4\sqrt{(n^{}+1)(n^{}+2)}}$$
(38)
and use the operator identity

$`r^2\widehat{}_{l+1}^+\widehat{}_l^+`$ $`=`$ $`[\widehat{b}_{n+1,l}^{}+2(n+1)+2]\times `$ (39)
$`\left[\widehat{b}_{n+2,l}^{}+2(n+2)\right]`$
to obtain the non-vanishing matrix elements

$`nl|r^2|n,l+2`$ $`=`$ $`{\displaystyle \frac{1}{2}}\sqrt{(2n+2l+3)(2n+2l+5)}`$ (40)
$`n+1,l|r^2|n,l`$ $`+`$ $`2=\sqrt{2(n+1)(2n+2l+5)}`$ (41)
$`n+2,l|r^2|n,l`$ $`+`$ $`2=\sqrt{(n+1)(n+2)}.`$ (42)
For the special case of $`l=0`$, we use the expansion of the irregular solutions in terms of $`|nl`$ Busch98a  and apply the above matrix elements. This results in

$`n0|r^2|n^{}2=\sqrt{{\displaystyle \frac{2\mathrm{\Gamma }(n^{}+\frac{7}{2})}{\pi \mathrm{\Gamma }(n^{}+1)}}}a\sqrt{{\displaystyle \frac{\nu _n}{a}}}`$ $`\times `$
$`\left({\displaystyle \frac{1}{n^{}\nu _n}}{\displaystyle \frac{2}{n^{}+1\nu _n}}+{\displaystyle \frac{1}{n^{}+2\nu _n}}\right)`$ . (43)
(The matrix elements $`n0|r^2|n^{}0`$ are not needed since $`I_{00}=0`$.)