# Spin monopoles with Bose–Einstein condenstates

## Abstract

We study the feasibility of preparing a Bose–Einstein condensed sample of atoms in a 2D spin monopole. In this state, the atomic internal spins lie in the $`x`$$`y`$ plane, and point in the radial direction.

After the successful generation of Bose–Einstein condensation of alkali atoms, the creation and manipulation of certain macroscopic quantum states remains as one of the fundamental goals in the field of Atomic Physics . During the last year, several ways to create vortices and solitons have been proposed . Under certain circumstances these states are stable , which has motivated several experimental groups to try corresponding experiments. In this letter we study a new kind of macroscopic quantum state for these atomic samples, what we call a 2D spin monopole. It is a state in which the atomic spin (i.e., the magnetization vector) points in the radial direction in the $`x`$$`y`$ plane \[Fig. 1(a)\]. We show that this state is stable under realistic conditions, and analyze a method to generate it which only requires current experimental technology. We will first study a one dimensional situation, in which the condensate is confined in a ring, and then we will generalize it to the 3D case.

We consider a Bose–Einstein condensed sample of $`N`$ atoms. The atoms have two internal levels, $`|`$ and $`|`$, which are relevant for the dynamics . We assume that half of the atoms is in each internal level. The atomic motion is confined by an external trap with a ring shape \[see Fig. 1(b)\] . In the limit where the motion along the radial ($`r`$) and axial ($`z`$) direction is frozen, the dynamics of the motional state only depends on the polar angle ($`\theta `$). We can write the wavefunction of the condensate as

$$|\mathrm{\Psi }(\theta ,\tau )=\varphi _1(\theta ,\tau )|+\varphi _2(\theta ,\tau )|,$$
(1)   
where $`\varphi _{1,2}`$ are the motional wavefunctions corresponding to the internal states $`|,`$, respectively and fulfill the coupled Gross–Pitaevskii Equations:

$`i{\displaystyle \frac{d}{d\tau }}\varphi _1`$ $`=`$ $`\left[{\displaystyle \frac{d^2}{d\theta ^2}}+\delta _1+u_{11}|\varphi _1|^2+u_{21}|\varphi _2|^2\right]\varphi _1,`$ (3)   
$`i{\displaystyle \frac{d}{d\tau }}\varphi _2`$ $`=`$ $`\left[{\displaystyle \frac{d^2}{d\theta ^2}}+\delta _2+u_{12}|\varphi _1|^2+u_{22}|\varphi _2|^2\right]\varphi _2,`$ (4)   
with normalization

$$\frac{1}{2\pi }_0^{2\pi }𝑑\theta |\varphi _{1,2}(\theta ,\tau )|^2=1.$$
(5)   
Here, $`\delta _{1,2}`$ denotes the energy of the two internal states and $`u_{ij}=u_{ji}`$ ($`i,j=1,2`$) describe their mutual interactions. All the quantities in Eqs. (1) have been rescaled so that they are dimensionless. In particular, the $`u_{ij}`$ are proportional to the number of atoms $`N`$, the corresponding scattering length, and the square of the ring radius \[see Eq. (22) below\].

The state (1) with

$$\varphi _1^{\mathrm{mp}}(\theta ,\tau )=e^{i\mu _1\tau },\varphi _2^{\mathrm{mp}}(\theta ,\tau )=e^{i\theta }e^{i\mu _2\tau },$$
(6)   
is a stationary solution of Eqs. (1) with $`\mu _j=j1+u_{1j}+u_{2j}+\delta _j`$ ($`j=1,2`$). Furthermore, defining the Pauli operator $`\stackrel{}{\sigma }=(\sigma _x,\sigma _y,\sigma _z)`$ as usual, we have that the state (1) with (6) is an eigenstate of $`\stackrel{}{\sigma }\stackrel{}{n}`$ with $`\stackrel{}{n}=[\mathrm{cos}(\theta \mu \tau ),\mathrm{sin}(\theta \mu \tau ),0]`$ and $`\mu \mu _2\mu _1`$. This means that for $`\mu =0`$ the spins are always pointing outwards in the $`x`$$`y`$ plane, whereas for $`\mu 0`$ they oscillate as shown in Fig. 1(a). Thus, the stationary state (6) can be considered as an oscillating spin monopole in two dimensions.

In order to analyze the stability of the spin monopole we have carried out a perturbative linear analysis. We consider a small perturbation around the monopole solution, such that $`\varphi _{1,2}(\theta ,\tau )=\varphi _{1,2}^{\mathrm{mp}}(\theta ,\tau )+ϵ\alpha _{1,2}(\theta ,\tau )`$. We then expand Eq. (1) up to first order in $`ϵ`$ obtaining a linear set of coupled equations for $`\alpha _{1,2}`$ and $`\alpha _{1,2}^{}`$. Expanding $`\alpha _{1,2}(\theta ,\tau )=_{n=\mathrm{}}^{\mathrm{}}\alpha _{1,2}^{(n)}(\tau )e^{in\theta }`$ and substituting in these equations, we obtain

$$i\frac{d}{d\tau }\stackrel{}{\alpha }^{(n)}=E_n\stackrel{}{\alpha }^{(n)},$$
(7)   
where

$$\stackrel{}{\alpha }^{(n)}=[\alpha _1^{(n)},\alpha _1^{(n)},\alpha _2^{(n+1)},\alpha _2^{(n+1)}]^\mathrm{T},$$
(8)   
is a column vector. Here, $`_nK_n+H^{\mathrm{int}}`$, with

$`E`$ $`=`$ $`\mathrm{diag}(1,1,1,1),`$ (10)   
$`K_n`$ $`=`$ $`\mathrm{diag}[n^2,n^2,(n+1)^21,(n1)^21],`$ (11)   
$`H^{\mathrm{int}}`$ $`=`$ $`\left(\begin{array}{cc}u_{11}& u_{21}\\ u_{12}& u_{22}\end{array}\right)\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right),`$ (16)   
are $`4\times 4`$ matrices and $``$ denotes tensor product. The stability analysis can be carried out by diagonalizing the matrices $`E_n`$. Complex eigenvalues $`\pm i\lambda `$ correspond to exponentially growing solutions \[$`\alpha (\tau )=\alpha (0)e^{\lambda \tau }`$\], i.e. the stationary solution is dynamically unstable. Real eigenvalues $`\lambda `$ lead to small oscillatory solutions \[$`\alpha (\tau )=\alpha (0)e^{\pm i\lambda \tau }`$\], which may correspond to higher or lower values of the Gross–Pitaevskii energy functional with respect to the stationary solution. In the first case, the stationary solution is stable, whereas in the latter one the system is energetically unstable; this means that if energy can be taken out from the system at a given rate, then the stationary solution is unstable on that time scale. Another interesting case occurs when one of the matrices is not simple (i.e., non diagonalizable), since in that case one has solutions which only grow polynomically with time \[for example, for a $`2\times 2`$ Jordan matrix corresponding to a degenerate eigenvalue $`\lambda `$, $`\alpha (\tau )=\alpha (0)(1+\lambda \tau )`$\]. On the other hand, if one is only interested in checking whether the solution is stable or not, (without distinguishing among the different types of instability) one can simply study the positivity of $`_n`$. If all matrices $`_n0`$ and the projector operator $`P_0^n`$ on the kernel of $`_n`$ commutes with $`E`$ for all $`n`$, then the stationary solution is stable. One can proof this statement as follows: first, the positivity of $`_n`$ ensures the positivity of the Gross–Pitaevskii energy functional and therefore there cannot be energetically instabilities; secondly, all eigenvalues of $`E_n`$ are also eigenvalues of $`_n^{1/2}E_n^{1/2}`$, which is hermitian (has real eigenvalues) and therefore there cannot be dynamical instabilities; finally, if $`[P_0^n,E]=0`$ one can easily prove that if a vector belongs to the kernel (range) of $`_n`$, then it also belongs to the kernel (range) of $`E_n`$, and therefore this last matrix is simple.

We start out by analyzing the positivity of $`_n`$. We consider two cases: (a) $`u_{11}u_{22}<u_{12}^2`$. The matrix $`H^{\mathrm{int}}`$ has a negative eigenvalue; for $`n=0`$ we have $`K_0=0`$ and therefore $`_0`$ has a negative eigenvalue and the solution is not stable. (b) $`u_{11}u_{22}>u_{12}^2`$. Now $`H^{\mathrm{int}}>0`$. Given that $`K_n0`$ for $`|n|>2`$ we can restrict ourselves to the case $`|n|=1`$. Since the first three diagonal minors of $`_n`$ are positive, we just have to impose that the determinant be positive. We obtain that the monopole (6) is stable if $`u_{11}u_{22}u_{12}^2+u_{22}/23u_{11}/23/4>0`$. Now, we investigate which kind of instabilities occur in the opposite case. We can again restrict ourselves to the case $`|n|=1`$, and therefore we have to find the eigenvalues $`\lambda `$ of the matrix $`E_1`$. In the relevant case for which $`u_{11}=u_{22}u`$ (see below) this can be done analytically. We obtain:

$$\lambda =1\pm \left[2\left(1+u\pm \sqrt{1+2u+u_{12}^2}\right)\right]^{1/2}.$$
(17)   
Thus, we have the following scenario:

$`\begin{array}{cc}\hfill u_{12}^2u^2u3/4& \text{Stable}\\ \hfill u^2u3/4<u_{12}^2<u^2& \text{Energetically unstable}\\ \hfill u^2<u_{12}^2& \text{Dynamically unstable.}\end{array}`$

The eigenvalues as well as the stability diagram is shown in Fig. 2. For $`u>u_{12}`$ and when the interaction energy becomes more important than the (rotational) kinetic energy ($`u1`$), the solution is stable as long as $`uu_{12}1`$. Since both $`u`$ and $`u_{12}`$ are proportional to the number of atoms $`N`$, we have that by increasing $`N`$ one can completely stabilize the spin monopole. An interesting situation occurs for $`u=u_{12}`$; in that case the matrix $`E_1`$ is not simple, which implies that the perturbation only grows linearly with time. On the other hand, the dynamical instability occurring for $`u<u_{12}`$ also appears for the homogeneous stationary solution and therefore it simply corresponds to a phase separation of the two components (internal states), as could be expected.

In order to establish if the monopole continues being stable in a realistic 3–dimensional situation, we consider a trapping potential of the form \[Fig. 1(c)\]

$$V(r,z)=\frac{1}{2}m\omega _z^2z^2+\frac{1}{2}m\omega _r^2r^2+V_0e^{r^2/(2\sigma ^2)}.$$
(18)   
This corresponds to a standard dipole trap with an off resonant Gaussian laser beam propagating along the $`z`$ direction . Here $`V_0`$ gives the ac–Stark shift at the center of the trap. The equilibrium point $`R`$ along the radial direction is given by $`R=2\sigma ^2\mathrm{ln}[V_0/(m\omega _r^2\sigma ^2)]`$. In order to have a ring trap, we choose the parameters such that $`R\mathrm{\Delta }r`$, where $`\mathrm{\Delta }r`$ is the typical size of the atomic cloud along the $`z`$ and $`r`$ directions. Writing $`|\mathrm{\Psi }(r,z,\theta )=\psi _1(r,z,\theta )|+\psi _2(r,z,\theta )|`$, the Gross–Pitaevskii Equation governing the evolution is

$`i\mathrm{}{\displaystyle \frac{d}{dt}}\psi _1`$ $`=`$ $`\left[{\displaystyle \frac{\mathrm{}^2^2}{2m}}+V+\mathrm{}\stackrel{~}{\delta }_1+\stackrel{~}{u}_{11}|\psi _1|^2+\stackrel{~}{u}_{12}|\psi _2|^2\right]\psi _1,`$ (19)   
$`i\mathrm{}{\displaystyle \frac{d}{dt}}\psi _2`$ $`=`$ $`\left[{\displaystyle \frac{\mathrm{}^2^2}{2m}}+V+\mathrm{}\stackrel{~}{\delta }_2+\stackrel{~}{u}_{21}|\psi _1|^2+\stackrel{~}{u}_{22}|\psi _2|^2\right]\psi _2.`$ (20)   
Here $`\stackrel{~}{u}_{ij}=4\pi \mathrm{}^2a_{ij}/m`$ with $`a_{ij}=a_{ji}`$ being the s–wave scattering lengths corresponding to the different collisions, and $`\stackrel{~}{\delta }_{1,2}`$ two constant offsets. The functions $`\psi _{1,2}`$ are normalized to the number of atoms $`N_{1,2}`$ in each internal state.

Let us first make the connection between the 3–D model and the ring. We consider the simple situation in which the motions along the radial and $`z`$ direction are frozen. The conditions of validity of this situation will be discussed below. If the number of atoms in each internal level is the same $`N_1=N_2N/2`$, we can reduce the full three dimensional problem to the ring case studied above by writing $`\psi _{1,2}(r,z,\theta ,t)=[N/(4\pi )]^{1/2}f_{1,2}(r,z)\varphi _{1,2}(\theta ,t)`$, multiplying the coupled Gross–Pitaevskii equations by $`f_i^{}`$ and integrating in $`r`$ and $`z`$. We obtain that $`\varphi _{1,2}`$ satisfy the Eqs. (1) with

$`u_{ij}`$ $`=`$ $`2NR^2a_{ij}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑z{\displaystyle _0^{\mathrm{}}}r𝑑r|f_i(r,z)|^2|f_j(r,z)|^2,`$ (22)   
$`\delta _i`$ $`=`$ $`{\displaystyle \frac{2mR^2}{\mathrm{}^2}}(\mathrm{}\stackrel{~}{\delta }_i+ϵ_i),\tau ={\displaystyle \frac{\mathrm{}}{2mR^2}}t,`$ (23)   
where $`ϵ_i`$ is the expectation value of the kinetic plus potential energy with the wavefunction $`f_i`$.

In general, we look for stationary solutions of the form

$`\psi _1(r,z,\theta ,t)`$ $`=`$ $`\sqrt{N/2}f_1(r,z)e^{i\stackrel{~}{\mu }_1t},`$ (25)   
$`\psi _2(r,z,\theta ,t)`$ $`=`$ $`\sqrt{N/2}f_2(r,z)e^{i\theta }e^{i\stackrel{~}{\mu }_2t},`$ (26)   
with $`f_{1,2}`$ normalized real functions satisfying $`\stackrel{~}{L}_{1,2}f_{1,2}=0`$ where

$`\stackrel{~}{L}_n`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}\left[{\displaystyle \frac{d^2}{dr^2}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{d}{dr}}+{\displaystyle \frac{d^2}{dz^2}}{\displaystyle \frac{(n1)^2}{r^2}}\right]+V(r,z)`$ (28)   
$`\mathrm{}\stackrel{~}{\delta }_n\mathrm{}\mu _n+\stackrel{~}{u}_{1n}f_1^2+\stackrel{~}{u}_{2n}f_2^2.`$   
We will concentrate in the case $`\stackrel{~}{u}_{11}=\stackrel{~}{u}_{22}\stackrel{~}{u}`$ (equivalently, $`a_{11}=a_{22}a`$) which is the relevant experimental situation described below. In the limit $`R\mathrm{\Delta }r`$ the centrifugal term in $`\stackrel{~}{L}_2`$ can be approximated by a constant $`\mathrm{}^2/(2mR^2)`$. If there is no phase separation ($`a>a_{12}`$) we find $`f_1(r,z)=f_2(r,z)f(r,z)`$ and $`\stackrel{~}{\mu }\stackrel{~}{\mu }_2\stackrel{~}{\mu }_1=\mathrm{}/(2mR^2)+\stackrel{~}{\delta }_2\stackrel{~}{\delta }_1`$. Solution (2) is again an eigenstate of $`\stackrel{}{n}\stackrel{}{\sigma }`$ with $`\stackrel{}{n}=[\mathrm{cos}(\theta \stackrel{~}{\mu }t),\mathrm{sin}(\theta \stackrel{~}{\mu }t),0]`$, being the eigenvalue proportional to $`|f(r,z)|^2`$. Thus it represents a state where the atoms have their spin in the $`x`$-$`y`$ plane forming a 2D monopole. The length of the spin at each point depends on the corresponding local density. We analyze now the stability of such a solution. As before, we linearize around the solution (2) by adding a small quantity $`\alpha _{1,2}`$ and expanding it in powers of $`e^{i\theta }`$. We obtain (7) where now $`_n=\stackrel{~}{L}+\stackrel{~}{K}_n+\stackrel{~}{H}^{\mathrm{int}}`$ a $`4\times 4`$ matrix operator with $`\stackrel{~}{L}=\mathrm{diag}(\stackrel{~}{L}_1,\stackrel{~}{L}_1,\stackrel{~}{L}_2,\stackrel{~}{L}_2)`$, $`\stackrel{~}{K}_n=\mathrm{}^2/(2mr^2)K_n`$, and $`\stackrel{~}{H}^{\mathrm{int}}=f^2H^{\mathrm{int}}`$. Using again $`R\mathrm{\Delta }r`$ we have $`L_1=L_2`$ and therefore $`\stackrel{~}{L}=L_1`$ times the $`4\times 4`$ identity matrix. Since we are only interested in the stability, we analyze the positivity of $`_n`$. As before, we just have to study the cases $`|n|<2`$. Let us distinguish two situations:

Weak interactions: in the limit $`N(a+a_{12})/R1`$ we have that the interaction energy $`N\stackrel{~}{u}/[2\pi R(\mathrm{\Delta }r)^2]`$ is much smaller than the harmonic oscillator quantum $`\mathrm{}\omega `$. In that case, $`\mathrm{\Delta }ra_0`$ where $`a_0=[\mathrm{}/(m\omega )]^{1/2}`$ is the size of harmonic potential ground state. The spectrum of $`_{0,1}`$ is dominated by $`\stackrel{~}{L}`$. The radial and $`z`$ dependence give rise to excitation energies $`k\mathrm{}\omega `$ ($`k`$ integer). The lowest excitations $`k=0`$ correspond to $`\alpha (r,z)f(r,z)`$, and therefore we obtain $`\stackrel{~}{\lambda }=\mathrm{}^2/(2mR^2)\lambda `$ ($`\mathrm{}\omega `$ in absolute value), where $`\lambda `$ is given in (17) with $`u`$ and $`u_{12}`$ given in (22). Thus, for excitation energies lower than $`\mathrm{}\omega `$ the problem fully reduces to the ring case.

Strong interactions: In the opposite limit, we are in the Thomas–Fermi regime, where $`\mathrm{\Delta }r=a_0[32N(a+a_{12})/R]^{1/4}`$. Now, one cannot simply separate radial and $`z`$ excitations from ring excitations. The excitation spectrum of $`_{0,1}`$ is dominated by $`\stackrel{~}{L}+\stackrel{~}{H}^{\mathrm{int}}`$. It is convenient to diagonalize $`\stackrel{~}{H}^{\mathrm{int}}`$, and consider the eigenfunctions separately. (a) Consider $`\stackrel{}{\alpha }^{(n)}=(g_1,g_2)(1,1)`$: in this case $`\stackrel{~}{H}^{\mathrm{int}}`$ is zero, and therefore the excitation frequencies correspond to those of $`\stackrel{~}{L}`$, which are of the order of $`k\mathrm{}^2/[2m(\mathrm{\Delta }r)^2]\mathrm{}^2/(2mR^2)`$. (b) For $`\stackrel{}{\alpha }^{(n)}=(g,g)(1,1)`$, $`H^{\mathrm{int}}`$ gives $`(\stackrel{~}{u}+\stackrel{~}{u}_{12})f^2`$, whereas for $`\stackrel{}{\alpha }^{(n)}=(g,g)(1,1)`$ it gives $`(\stackrel{~}{u}\stackrel{~}{u}_{12})f^2`$: the lowest energy will be of the order of $`N(\stackrel{~}{u}\stackrel{~}{u}_{12})/[2\pi R(\mathrm{\Delta }r)^2]`$. As long as this energy is larger than $`\mathrm{}^2/(2mR^2)`$ we can consider separately the cases (a) and (b) treating $`K_n`$ as a perturbation; in both the correction is positive, i.e. the monopole is stable. In the opposite case, one has to be more careful in the perturbation analysis, since one cannot separate the cases (a) and (b); the excitation energies may become negative. Thus, we obtain a necessary condition for stability $`N(\stackrel{~}{u}\stackrel{~}{u}_{12})/[2\pi R(\mathrm{\Delta }r)^2]\mathrm{}^2/(2mR^2)`$. Using (22) we can write this condition as $`uu_{12}>1`$, which coincides with the basic stability condition derived for the ring. It means that the interactions have to be sufficiently strong to stabilize the monopole.

In order to be specific, we will propose now a particular configuration to create the spin monopole. We consider an alkali atom in a ground $`F=1`$ hyperfine state. We will assume that the energy of the $`m_F=0`$ level is made higher (by using an off–resonant laser or radio-frequency field), so that it is not involved in the dynamics. In this case we can identify $`|=|F=1,m_F=1`$ and $`|=|F=1,m_F=1`$, the collisions do not change spin, and $`\stackrel{~}{u}_{11}=\stackrel{~}{u}_{22}\stackrel{~}{u}`$. Although the current status does not allow to make a clear statement about the sign of $`\stackrel{~}{u}\stackrel{~}{u}_{12}`$ it seems that for $`{}_{}{}^{23}Na`$ it is positive, which is also agrees with recent experimental results concerning miscibility . In order to generate the 2–D spin monopole (2) we propose to use an off-resonant Raman beam. The atoms are initially condensed in the internal $`|`$ state. A Raman laser that connects the states $`|`$ and $`|`$ is then switched on. It should have the appropriate spatial dependence so that the angular momentum in the $`z`$ direction is changed by one unit . Denoting by $`\mathrm{\Omega }(z,r,\theta )`$ the effective Rabi frequency, the evolution equations are the above Gross–Pitaevskii Equations but with a coupling term between $`\psi _1`$ and $`\psi _2`$ proportional to $`\mathrm{\Omega }`$; the Raman detuning is incorporated to the definitions of $`\stackrel{~}{\delta }_{1,2}`$. Initially, one takes $`\stackrel{~}{\delta }_1\stackrel{~}{\delta }_2`$, so that the laser does not affect the internal atomic state, since it is effectively out of resonance. Then, $`\stackrel{~}{\delta }_2\stackrel{~}{\delta }_1`$ is changed adiabatically until $`\stackrel{~}{\delta }_2\stackrel{~}{\delta }_10`$. The method is more robust than the one used to generate vortices since the spatial wavefunctions $`f_{1,2}(r,z)`$ remain practically constant with our setup. Actually, in this case one can simply use a $`\pi /2`$ laser pulse, taking $`\stackrel{~}{\delta }_1=\stackrel{~}{\delta }_2`$. This allows to generate the monopole state in a much faster time scale, which will be of the order of several inverse trap frequencies.

In order to evaluate our proposal, we have performed a 3–dimensional numerical simulation of the Gross–Pitaevskii equations in the presence of the laser for the creation of the spin monopole. We have used an optimized three-dimensional collocation Fourier method with typically 80$`\times `$80$`\times `$40 collocation points and integrating in time with a symmetrized split-step operator technique. The results are shown in Fig. 3. In Fig. 3(a) we have plotted the final density profile along the $`x`$ axis at the end of the generation process, whereas in Fig. 3(b) we have plotted the evolution of the population of the internal levels until the monopole is generated. After this process, we switch of the laser and apply an internal energy shift so that the spin start precessing, as is shown in Fig. 3(c). For this figure we have taken realistic parameters. For example, taking Na with $`\omega _r=100`$Hz, $`a_0=2\mu `$m, and $`a=52a_B`$ we have that the number of atoms is of the order of $`2\times 10^5`$. In order to ensure the stability, we have evolved the formed state in imaginary time (renormalizing the state after each evolution step) for about 20 trap oscillation times, without noticing any instability.

In summary, we have shown that when the interactions are sufficiently strong, a 2D spin monopole becomes stable in a ring trap. We have performed the stability analysis, both in a ring situation as well as in the full 3D case. We have shown a way to prepare such a state, and verified with a full 3D numerical simulation that one can prepare it with current experimental parameters.

J.J.G-R and V.M.P-G. have been partially supported by CYCIT grants PB96-0534 and PB95-0839. The work in Innsbruck has been supported by the Austrian Science Foundation.