# Numerical problem with Hartree-Fock equations for dilute Bose gas

I have to solve the following set of equations self-consistently: \begin{align} n_c(\mathbf{r}) & = \frac{1}{g}\left[\mu - V_{\rm ext}(\mathbf{r}) - 2 g n_{T}(\mathbf{r}) \right] \\[3mm] n_{T}(\mathbf{r})& = \frac{1}{\lambda_{T}^{3}} Li_{3/2}\left(e^{\frac{\mu - V_{\rm eff}(\mathbf{r})}{k_B T}} \right) \\[3mm] V_{\rm eff}(\mathbf{r})& = V_{\rm ext}(\mathbf{r}) + 2 g [n_c(\mathbf{r}) + n_{T}(\mathbf{r})] \\[3mm] \mu &= g\left[n_c(0) + 2n_{T}(0) \right] \end{align} where $g$ and $\lambda_{T}$ are constants (apart from temperature $T$ and Boltzmann constant $k_{B}$). $Li_{3/2}(x)$ is the polylogarithm function. This is basically the Hartree-Fock equations for Bose-Einstein condensate $n_c$ and thermal cloud of excitations $n_T$. I tried to solve it iteratively but without success.

# UPDATE

My method is the following:

I start with thermal density $n_T$ set exactly to zero and adjust $\mu$ in such a way that integral over coordinate space with $n_c$ gives $1$. Then I start the iteration loop:

1) Having the chemical potential and effective potential I update the thermal density using Eq. no. 2 2) I find chemical potential and condensate density in such a way that the integral over coordinate space with $n_c$ given by Eq. no. 1 is $1$.

Repeat the process until final convergence is reached (I compare condensate density - integral over coordinate space before and after the iteration step). I set the precision to $1.0E-11$, but sometimes it stop at $1.0E-05$ and never decreases.

Here is my quick python code: http://rextester.com/BNZM39456

Do you know of any source where I can find a code for this problem and compare with mine?

• You say it does not converge. What do you get specifically after many iterations? – Lewis Miller Dec 24 '15 at 15:59
• @LewisMiller I did some update and provided my Python code for this problem. My simulations for small temperatures looked like ideal gas case. – WoofDoggy Dec 25 '15 at 19:28
• It looks like you get convergence to 5 significant digits but were hoping for 11. I took a look at your code and noticed that you are using 64 grid points in each direction for your numerical integration. Have you tried increasing that number? It is possible that 5 digits is all you can expect for that number of grid points. – Lewis Miller Dec 27 '15 at 21:32