Two dipolar atoms in a harmonic trap

We observe a remarkable progress in experiments with ultracold quantum gases. Many are performed with large number of atoms in a single trap. However, progress is also made at a level of a few atoms in a trap. These experiments are performed with cold atoms distributed between the wells of an optical lattice. This way, with a help of tunable parameters of interaction, using the Feshbach resonances [1,2], and the properties of the lattice itself, one can access in a controlled way vital models of condensed-matter physics (for a recent review see [3]). Many experiments in optical lattices are performed with the Mott insulator phase [4–6] where a well-defined, small number of atoms is confined in every well. Another set of a few atoms in a trap experiments is offered by the setting available in Selim Jochim's lab [7]. Detailed properties of such systems crucially depend on the properties of atom-atom interaction. This interaction is best tested if exactly two atoms are present. Early analytic predictions for contact interacting atoms [8] where positively verified in precise spectroscopic experiments [9]. New twist to the problem is introduced by the long-range dipole-dipole (DD) interactions [10,11]. Although negligible in the early days of quantum gases experiments, dipole-dipole interactions are getting more and more relevant with the condensation of chromium [12,13], erbium [14–16] and recently dysprosium [17–20]. The dipolar interaction couples the spin degree of freedom with the orbital angular momentum. This leads to the well-known Einstein-de Haas effect [21]. To observe this effect with chromium atoms, where DD interaction is just a perturbation, properly resonant magnetic field strength must be used [22]. Of course a direct coupling to the orbital angular momentum is possible for sufficiently strong DD interactions. For the large systems it has been noted using a conventional mean-field approach [23]. It is the purpose of this letter to present are exact analysis of the role of DD interactions for two atoms trapped in a spherically symmetric harmonic potential. Initially we introduce our theoretical model. The simplicity of the harmonic potential allows to separate the center-of-mass degree of freedom. The relative-motion Hamiltonian remains spherically symmetric. Utilizing this symmetry we may construct the energy eigenstates using the angular-momentum algebra. What remains is the set of coupled radial Schrödinger equations linking components of the wave function corresponding to orbital angular momenta differing by two units. We model the radial component of the interaction by the dipolar expression modified by the infinite repulsive sphere at short distances. Finally we present our results. We note the anti-crossings of the energy levels as a function of the dipolar coupling constant. This dependence may be tuned by the change of the trapping frequency. The most striking feature is the possibility of adiabatically pumping our system from the s-wave to the d-wave relative motion.

Let us consider two identical dipolar atoms (fermions or bosons) of a spin (a total angular momentum of an atom) $f_{1}=f_{2}$ moving in an isotropic harmonic trap. The Hamiltonian of such a system can be written as

$$\begin{equation} H=-\frac{1}{2}\nabla_{1}^{2}-\frac{1}{2}\nabla_{2}^{2}+\frac{1}{2}r_{1}^{2}+\frac{1}{2} r_{2}^{2}+V(\mathbf{r}_{1}-\mathbf{r}_{2}), \end{equation} \tag{ 1 }$$

where $\mathbf{r}_{1}$ and $\mathbf{r}_{2}$ are the position vectors of the two atoms. We are using harmonic-oscillator units, in which $\hbar \omega$ is a unit of energy and the characteristic size of the ground state of the trap $\sqrt{\frac{\hbar}{m\omega}}$ is a length unit. An interaction potential $V(\mathbf{r}_{1}-\mathbf{r}_{2})$ is a sum of a short-range (SR) and a long-range magnetic dipole-dipole interaction potentials. We model SR potential as a spherically symmetric barrier described later in the next section (for different SR models see, e.g. [8,10,11,24–26]). The magnetic dipole-dipole interaction potential $V_{dd}(\mathbf{r}_{1}-\mathbf{r}_{2})$ can be expressed in the following form:

$$\begin{eqnarray} V_{dd}&= \frac{\mu_{0}(\mu_{B}g_{j})^{2}}{4\pi\left | \mathbf{r}_{1} - \mathbf{r}_{2} \right |^{3}} \Bigl[\mathbf{F}_{1}\cdot \mathbf{F}_{2}-3\left(\mathbf{F}_{1} \cdot \mathbf{n} \right) \left( \mathbf{F}_{2} \cdot \mathbf{n} \right) \Bigr], \end{eqnarray} \tag{ 2 }$$

where $\mathbf{n}=\frac{\mathbf{r}_{1}-\mathbf{r}_{2}}{\left| \mathbf{r}_{1} - \mathbf{r} _{2}\right|}$, $\mu_{0}$ stands for the vacuum magnetic permeability, $\mu_{B}$ indicates the Bohr magneton, g_j is the Landé g-factor and F is the total angular momentum of an atom (spin vector). Thus for the atomic spin quantum number f half integer we have fermions and for f integer we have bosons.

Foregoing Hamiltonian may be divided into two parts, a center-of-mass part and a relative part, i.e. $H=H_{CM}+H_{\textit{Rel}}$ with

$$\begin{equation} \begin{array}{rcl} H_{CM}&=&-\displaystyle \frac{1}{2}\nabla_{R}+\frac{1}{2}\mathbf{R}^{2},\\ H_{\textit{rel}}&=&-\displaystyle \frac{1}{2}\nabla_{r}+\frac{1}{2}\mathbf{r}^{2}+ V_{VdW}(r) \\ &&+\,\displaystyle \frac{g_{dd}}{r^{3}}\left[ \mathbf{F}_{1}\cdot \mathbf{F}_{2}-3 \left( \mathbf{F}_{1} \cdot \mathbf{n} \right) \left( \mathbf{F}_{2} \cdot \mathbf{n} \right) \right] \end{array} \end{equation} \tag{ 3 }$$

which can be diagonalized separately. Here $\mathbf{R}=\frac{1}{\sqrt{2}}\left( \mathbf{r} _{1}+\mathbf{r}_{2} \right)$ is the center-of-mass coordinate and $\mathbf{r}=\frac{1} {\sqrt{2}}\left( \mathbf{r}_{1}-\mathbf{r}_{2} \right)$ stands for the relative-motion coordinate¹. The strength of the dipole-dipole interaction is characterized by the $g_{dd}= \frac{\mu_{0}(\mu_{B}g_{j})^{2}}{8 \sqrt{2} \pi}$. The eigenvalues of the H_CM are simply that of the harmonic oscillator. In order to investigate the relative motion of the two atoms we observe that the total angular momentum is conserved:

$$\begin{equation} \left[\mathbf{F}+\mathbf{L},H_{rel} \right]\equiv\left[ \mathbf{J},H_{\textit{rel}} \right]= 0, \end{equation} \tag{ 4 }$$

where J stands for the total angular momentum of the system which is a sum of the total spin operator $\mathbf{F}=\mathbf{F}_{1}+\mathbf{F}_{2}$ and the orbital-momentum operator of the relative motion of the atoms L. The spherical symmetry of the system means that it is convenient to solve the relative-motion problem in a total-angular-momentum basis. The eigenfunction of the system in the chosen basis can be written as

$$\begin{eqnarray} \Psi^{jm_{j}}_{n}(\mathbf{r})&= \sum_{l,f}a^{jm_{j}lf}_{n}\psi^ {jm_{j} lf}_{n}\left( \mathbf{r} \right)\nonumber\\ &= \sum_{l,f}a^{jm_{j}lf}_{n}\phi_{n}^{jlf}(r)\left|jm_{j}lf\right\rangle\nonumber\\ &= \sum_{l,f}a^{jm_{j}lf}_{n}\phi_{n}^{jlf}(r)\sum_{\substack {m_{l},m_{f}\\&&m_{l}+m_{f}=m_{j}}}\nonumber\\ & \quad \times\, C^{jm_{j}}_{lm_{l}fm_{f}} \left|lfm_{l}m_{f}\right\rangle. \end{eqnarray} \tag{ 5 }$$

Here j denotes the total-angular-momentum quantum number and m_j the magnetic total-angular-momentum number, l and m_l stand for the orbital-momentum and the magnetic orbital-momentum quantum numbers, respectively. The total spin and its projection values are indicated by f and m_f and $C^{jm_{j}}_{lm_{l}fm_{f}}$ denotes Clebsch-Gordan coefficients [27]. Eigenfunctions are enumerated by the $n=0,1,\ldots$ number and $a^{jm_{j}lf}_{n}$ indicate constant coefficients.

Our goal is to derive the radial Shrödinger equations for $\phi^{jlf}_{n}$ with given j, l, f. We are now interested in the result of acting with the V_dd operator on the single state $\psi^{jm_{j}lf}_{n}(\mathbf{r})$. Using spin operators properties it can be shown that

$$\begin{equation} V_{dd}\psi^{jm_{j}lf}_{n}(\mathbf{r})=\frac{g_{dd}}{r^{3}}\sum_{l',f'}\alpha_{ll'ff'}\psi ^{jm_{j}l'f'}_{n}(\mathbf{r}) \end{equation} \tag{ 6 }$$

with the following selection rules:

$$\begin{equation} \begin{array}{rcl} l'=l+\Delta l, \qquad\Delta l=0,\pm 2,\\ f'=f+\Delta f, \quad\,\Delta f=0,\pm 2. \end{array} \end{equation} \tag{ 7 }$$

A value of the scalar coefficient $\alpha_{ll'ss'}$ is expressed by a product of Clebsh-Gordan coefficients determined by the angular-momentum algebra.

The preceding result might be understood by the fact that the dipole-dipole interaction operator is symmetric with respect to the exchange of the two particles. Thus it does not change a symmetry of the given $\psi^{jm_{j}lf}_{n}(\mathbf{r})$. Knowing (6) we are able to find the radial Schrödinger equation for the $\chi_{n}^{jlf}(r)\equiv r\phi_{n}^{jlf}(r)$ by the direct calculation

$$\begin{eqnarray} & -\,\frac{1}{2}\frac{\text{d}^{2}}{\text{d}r^{2}}\chi_{n}^{jlf}(r)+\frac{1}{2}r^{2}\chi_{n}^{jlf}(r)+ \frac{l(l+1)}{2r^{2}}\chi_{n}^{jlf}(r) \nonumber\\ & +\,\frac{g_{dd}}{r^{3}}\sum_{l',f'}\alpha_{ll'ff'}\chi^{jl'f'}_{n}(r)=E^{j}_{n}\chi_{n} ^{jlf}(r), \end{eqnarray} \tag{ 8 }$$

where $E^{j}_{n}$ is an eigenvalue and the short-range potential $V_{VdW}(r)$ is incorporated in the boundary conditions (see the next section). Here $\sqrt{\frac{\hbar ^5}{m^3\omega}}$ is a unit of the g_dd in the harmonic-oscillator units.

As can be seen in (8) in order to find a $\chi_{n}^{jlf}(r)$ one has to solve a system of the radial Schrödinger equations for a fixed total angular momentum number j. Note that the number of equations in the system is determined by the maximum value of the total spin: $f_{max}=f_{1}+f_{2}$.

We are interested in solving the system of the radial Schrödinger equations introduced in the previous section, in particular for the total angular momentum j = 0. In order to accomplish this task we remind that the interaction potential $V(\mathbf{r}_{1} -\mathbf{r}_{2})$ introduced earlier consists of the short-range potential also. We use a simple model of such a potential in the following form [26]:

$$\begin{equation} V_{SR}=\begin{cases} 0, &\text{for }r > b=100 \,r_{0}, \\ \infty,& \text{for } r \leq b=100 \,r_{0}, \end{cases} \end{equation} \tag{ 9 }$$

where r₀ is the Bohr radius. We motivate our choice by the fact that the scattering length a₀ for a scattering process of a single particle on an infinite spherically symmetric potential barrier is equal to the radius of barrier, i.e. $b=a_{0}$. A value of b is determined by the numerical calculations for the dysprosium atoms [28].

From the angular-momentum algebra we also deduced that for eigenstate with j = 0 the total spin number is equal to the orbital quantum number, i.e. l = f. Thus for such states the corresponding coefficient matrix $\alpha_{ll'ff'}$ reduces to the $\alpha _{ll'}$ matrix. We calculate them for the various atomic spin values, i.e. $f_{1}=f_{2}=\frac{1}{2},1,\frac{3}{2} \text{ and } \frac{21}{2}$. Our results can be found in the supplementary material supplementary.pdf.

Knowledge of the $\alpha_{ll'ff'}$ coefficients allows us to solve numerically the system of the radial Schrödinger equations of the form presented in (8). We used the multi-parameter shooting method. We set the $b=0.04$ in the harmonic-oscillator units which corresponds for the dysprosium-like atoms at the trap frequency $\omega \approx 2\pi\ 2.8\ \text{kHz}$. For such a trap frequency the $g_{dd}=0.0006$ in the harmonic-oscillator units. Our system admits two control parameters that may be changed by experimenters. Note that the g_dd in the harmonic-oscillator units depends on the trap frequency as $\sqrt{\omega}$, so it is tunable. One may also change the scattering length a₀ by the optical Feshbach resonances [29–32], so that the b value in the harmonic-oscillator units may be kept constant while one changes the trap frequency.

In fig. 1 we present the eigenvalues $E^{0}_{n}$ with $n=0,1,2$ as a function of g_dd for atoms with different spins. For atoms with the spin $f_{1}=f_{2}=1$, $\frac{3}{2}$ and $\frac{21}{2}$ we consider only solutions for the even orbital-angular-momentum quantum number l. In the case of odd l results are qualitatively the same.

Zoom In Zoom Out Reset image size

For spin $\frac{1}{2}$ atoms the energy values rise very slowly as g_dd rises. The radial part of $\psi^{0011}_{n}(\mathbf{r})$ is simply the $\phi_{n}^{011}(r)$, so the expected value of the orbital-angular-momentum operator $\langle L^{2}\rangle$ is constant and equal $\langle L^{2}\rangle =2$ for all n.

For the higher spin values we observe more complex behaviour. First of all, the energy values for n = 0, 1 and 2 are highly dependent on the value of g_dd. For low values of g_dd eigenvalues vary slightly, then for higher values they decrease rapidly. We observe also the presence of anti-crossings between consecutive lines $E^{0}_{n}(g_{dd})$ accompanied by changes of the $\langle L^{2}\rangle$. This is due to changes in the structure of the radial part of eigenstates. From (8) we notice that the radial part of the eigenstate is a linear combination of $\phi_{n}^{0ll}(r)$ where in this case $l \in \{0,2,\ldots,2\cdot f_{1} \}$. As the g_dd rises the weight of each $\phi_{n}^{0ll}(r)$ function varies, i.e. the values of the $a^{00ll}_{n}$ coefficient vary. For instance, we see that for low g_dd the ground state consists of almost only the s-state $(\phi_{0}^{000}(r))$, whereas as we increase the trap frequency, contributions of the functions with higher l grow. The ground state starts to "rotate". This feature resembles the Einstein-de Haas effect [21], although it is caused only by the internal spin-spin interactions between two atoms without any influence of external fields.

Moreover, as values of $a^{00ll}_{n}$ for l > 0 grow and $a^{0000}_{n}$ decreases also mutual orientation between the atoms starts to favour attractive regions of the DDI over repulsive regions, which can be observed in the analysis of the angular part of the eigenstates (see also fig. IIa in [33]). For the spin $f_{1}=f_{2}= \frac{1}{2}$ atoms, such a behaviour is impossible as all three angular parts of the eigenstates are multiplied by the same $a^{0000}_{n}\phi_{n}^{000}(r)$ expression which is almost being unchanged as g_dd rises. This fact explains qualitative difference between the dependence of the eigenvalues on g_dd for spin $\frac{1}{2}$ and higher spin values.

Figure 1 also illustrates that the bigger atomic spin is, the lower trap frequency is needed to observe above effects. In addition, the effect of changes in the expected value of orbital angular momentum is stronger for larger atomic spin values. It seems that at least it is possible to check our model experimentally using the system of the dysprosium atoms with the $\frac{21}{2}$ spin.

The nature of anti-crossings in fig. 1 can be explained by Landau-Zener theory [34,35] as depicted in fig. 2. As an example we used spin-$\frac{3}{2}$ atoms. A composition of the eigenstate corresponding to the eigenvalue $E_{n}(g_{dd})$ is not conserved along a given energy line, but it propagates along straight lines upward or downward. This type of effect was already observed by Kanjilal et al., although for a simpler system consisting of two aligned dipoles [10].

Zoom In Zoom Out Reset image size

Motivated by experiments under development [18,19,36] we based our calculations on dysprosium parameters. Our model of the dipole-dipole interactions between two atoms reveals a non-trivial dependence of two atoms in a harmonic-trap system on the trap frequency. We showed that increasing ω the system undergoes an analog of Einstein-de Hass effect. Such a behaviour is a result of spin-spin interaction and its coupling to the orbital angular momentum. We have also found the Landau-Zener anti-crossings in the energy levels of the system. Our results may be checked experimentally for the dysprosium atoms. Of course, proposed model is oversimplified in this case as dysprosium atoms are not exactly spherically symmetric [28].

The authors acknowledge fruitful conversations with Mariusz Gajda. This work was supported by the (Polish) National Science Center Grant No. DEC-2012/04/A/ST2/00090.