# BEC in a rotating disc

Goodmorning everybody, I have to run a numerical simulation of a Bose-Einstein condensate on a rotating disc. Now, my problem is that I became suspicious about the equation I'm using, since the final result I obtain is not what I expected. First of all, the Gross-Pitaevskii equation (in a z-axis symmetrical trap) for the amplitude of the wave function of a condensate with a vortex is $$\mu f(\rho,z) = -\frac{\hbar^2}{2m}\left( \frac{1}{\rho}\frac{\partial}{\partial \rho}\left( \rho \frac{\partial f(\rho,z)}{\partial \rho} \right) - \frac{f(\rho)l^2}{\rho^2}+\frac{\partial^2 f(\rho,z)}{\partial z^2} \right) + V(\rho,z)f(\rho,z) + U_0f^2(\rho,z)\cdot f(\rho,z)$$ where $l$ is the quantum of circulation, $\mu$ the chemical potential and $U_0 = 4\pi\hbar^2a/m$ ($a$ is the scattering length). Now, I want to consider the equation on a disk, so I put $f=f(\rho)$. To solve it numerically, I need to perform the substitution $$f(\rho) = \frac{u(\rho)}{\sqrt \rho }$$ (while in a spherical confinment I would have put $f(\rho) = u(\rho)/\rho)$). The kinetic term is therefore given by $$\frac{1}{\rho}\frac{\partial}{\partial \rho} \left( \rho \frac{\partial f(\rho)}{\partial \rho} \right) = \frac{1}{\rho}\frac{\partial f(\rho)}{\partial \rho}+\frac{\partial^2f(\rho)}{\partial\rho^2} = \frac{u''(\rho)}{\sqrt \rho}+\frac{u(\rho)}{4\rho^{5/2}}$$ so that the first derivative is absent (I need that because I'm solving it using the Numerov's algorithm). Putting all together $$\mu u(\rho) = -\frac{\hbar^2}{2m}\frac{d^2u(\rho)}{d \rho^2} + \frac{\hbar^2}{2m\rho^2}\left( l^2-\frac{1}{4} \right)u(\rho) + V(\rho,z)u(\rho) + \frac{U_0}{\rho}u^2(\rho)\cdot u(\rho)$$ where, as external potential I use $V(\rho) = 1/2 m \omega^2\rho^2$. It seemed a good equation, until I noticed that if I try to test it in the case of the non-rotating disc ($l=0$, in this case the equation above must become exactly the GP equation in cylindrical coordinates) and non-interacting gas ($a=0$, that means $U_0=0$) the equation mantains a centrifugal term, and so the ground state solution vanishes as $\rho$ approches 0, that I don't think it's correct. What do you think? What did I do wrong? I assume that my program works, because the simulation in the spherical trap ended well (and I paid attention to modify the normalization condition and other stuff), but if you don't find any errors in my calculations, then I will go check again the algorithms! Thank you!!

• Your centrifugal term goes as $(l^2-\frac14)$ in the equation for $u$, so it changes sign from a centrifugal barrier to a well at $l=0$. Are you sure that this still impedes the formation of a ground state? Aug 9, 2013 at 14:09
• More generally, I would start by factoring out the angular momentum dependence as $f(\rho)=\rho^lg(\rho)$ to do away with the bad $1/\rho^2$ term before actually solving the equation. That should at least tell you the behaviour to expect from your solution for $u$ at $\rho\rightarrow0$ in the cases $l=0$ and $l>0$, even if you do need to use $\rho$ for numerical purposes. Aug 9, 2013 at 14:12
• The derivation seems correct. But are you really sure that you have added back the factor $\sqrt{\rho}$ at the end of your calculation/algorithm? Aug 10, 2013 at 14:06
• Also, when $l=0$, there is a $-\infty$ singularity at origin, how do you solve this numerically? Aug 10, 2013 at 14:14
• First of all, thank you for your answers!! @hwlau: there's no problem in solving it numerically because you impose that the solution $u$ is 0 in the origin, so that the divergence does not occur (is the same trick you use to solve the Schrodinger equation for example in spherical coordinates for the exited states, and it works believe me..). Reguarding the $\sqrt(\rho)$, I'm sure I factorized my solution correctly, because I created a special function to normalize the wave function and the division by $\sqrt \rho$ is done there! Aug 10, 2013 at 16:33