# Stability of dark solitons in a harmonic trap

This question is based upon a research article which I am trying to reproduce. One of the main result of this paper is the condition on transverse confinement of the Bose-Einstein Condensate(BEC) to make the black soliton solution stable. The equation governing BEC is the Gross-Pitaevskii(GP) equation, given by

$i\hbar \frac{\partial\psi}{\partial t}=-\frac{\hbar^{2}}{2m}\frac{\partial^{2}\psi}{\partial x^{2}}+g\psi|\psi|^{2}+V_{ext}\: \psi$

Here, $|\psi|^{2}$ gives the density of the condensate. We can see when $V_{ext}=0$, the above equation has a solution, in one dimension(say z), of the form $\psi(z,t)=\tanh{cz} \:e^{-i\mu t}$, where $c$ accounts for the constants.

Suppose now that we are working in a cylindrical geometry, such that the $V_{ext}=\omega_{z}z^{2}+\omega_{r}r^{2}$ and $\omega_{z}<\omega_{r}$, meaning the radial confinement is stronger than the axial confinement. In such a case, one can obtain a solitonic solution with a nodal plane perpendicular to the axial direction. This can be done by using the split operator method and imaginary time evolution.

Now, comes the question of stability of the solitonic solution. One can perturb $\psi$ with a perturbation of the form $\psi\rightarrow\psi+\delta\psi$, where $\delta\psi = u(z)e^{iq.r-i\epsilon t}+v(z)e^{-iq.r+i\epsilon t}$. So, basically, we are looking for small amplitude oscillations where the soliton, whose nodal plane is perpendicular to the axial(z) direction, gives out energy in the radial direction. Putting the form of $\delta\psi$ in the GP equation, we get the following set of equations for $f_{\pm}(z)=u(z)\pm v(z)$, where $\psi_{0}(z)$ is the density profile of the soliton in the z-direction(axial direction) in the presence of the trap($V_{ext}$)

$\epsilon f_{-}(z)=\Big[-\frac{\hbar^{2}}{2m}\big(\frac{\partial^{2}}{\partial z^{2}}-q^{2}\big)-\mu+V_{ext}+3g\psi_{0}^{2}(z) \Big]f_{+}(z)$

$\epsilon f_{+}(z)=\Big[-\frac{\hbar^{2}}{2m}\big(\frac{\partial^{2}}{\partial z^{2}}-q^{2}\big)-\mu+V_{ext}+g\psi_{0}^{2}(z) \Big]f_{-}(z)$

To obtain a stability condition, we need a dispersion relation between $\epsilon$ and $q$. However, as you can see in the above set of equations, there are a lot of terms with $z$ dependences, including a derivative in the $z$ direction. The authors of the paper say that they have numerically solved these set of equations to obtain a dispersion relation. How does one do that?

• You are, basically, asking that somebody explain that paper to you. Furthermore, that paper does not seem to be a computational one, it is more related with solving a differential equation (Eq. 11). If you have particular questions about it, then this might be the right place to ask, otherwise I think that you are not in the right track here. – nicoguaro May 8 '17 at 20:21
• @nicoguaro If you read my question, I have explained the paper in one line already. So it should be clear that I am not asking somebody to explain the paper to me. Since you have misunderstood my question, I shall take your advice and write the equations here itself instead of giving a link. – Abhijit May 9 '17 at 5:48
• What are the boundary conditions and domain for the variable $z$? – nicoguaro May 9 '17 at 16:16
• @nicoguaro That depends on the length scale given by GP equation. This length is called the correlation length $\xi$. Roughly, the idea is to extend the domain of the $z$ such that the base function($\psi_{0}(z)$), on top which perturbations are introduced, is more or less constant even if you extend $z$ more. The boundary condition on $\psi_{0}(z)$ is that it goes to $0$ as $z$ goes to infinity in presence of trap($V_{ext}$). The seed usually used to find $\psi_{0}(z)$ using split-operator method is $\tanh{z}$. Since $f_{\pm}(z)$ are fluctuations, same boundary and domain conditions apply. – Abhijit May 10 '17 at 5:46

As short answer, you need to:

1. Sweep over your wavenumber $q\in [q_\min, q_\max]$;

$$\mathcal{L}_q \mathbf{f} = \epsilon \mathbf{f}$$

to obtain a discrete problem

$$[H(k)]\lbrace U \rbrace = E[S]\lbrace U \rbrace\, ;$$

1. Solve the generalized eigenvalue problem to obtain the admissible frequencies $\epsilon$ for that particular $q$.

Regarding point 2 and 3, there are several methods. Some of the most popular ones are discussed in this answer [1]:

• Finite Difference Methods;
• Finite Element Methods; and
• Direct Variational methods.

Let's rewrite the equation as above, $\mathcal{L}_q \mathbf{f} = \epsilon \mathbf{f}$, with

$$\mathcal{L}_q = \begin{bmatrix} 0 &-\frac{1}{2}\nabla^2 + \hat{V}(q) + 3\psi^2_0\\ -\frac{1}{2}\nabla^2 + \hat{V}(q) + \psi^2_0 &0\ \end{bmatrix}\, ,\quad \mathbf{f} = \begin{Bmatrix} f_{-}\\ f_{+}\end{Bmatrix}\, ;$$

and $\hat{V(q)} = V_\text{ext} - \mu - q^2$, where I re-scaled some variables for convenience. The operator $\mathcal{L}$ does not seem to be self-adjoint, if that's the case, the resulting matrices might not be Hermitian (depending on the scheme).

If one uses a weighted residual approach, the resulting functional is of the form (I would double-check, it's easy to make silly mistakes)

$$\Pi_q[\mathbf{f}, \mathbf{w}] = \frac{1}{2}\int_\mathbb{R} \nabla w_1 \nabla f_{+}\, dz + \frac{1}{2}\int_\mathbb{R} \nabla w_2 \nabla f_{-}\, dz +\\ \int_\mathbb{R} |\psi_0|^2[3 w_1 f_{+} + w_2 f_{-}]\, dz + \int_\mathbb{R} \hat{V}(q)[w_1 f_{+} + w_2 f_{-}]\, dz - \epsilon\left[\int_\mathbb{R} w_1 f_{-}\, dz + \int_\mathbb{R} w_2 f_{+}\, dz\right]$$

• Ok. I had tried the finite difference method, but I could not reproduce the results. There is a small difference between what the link [1] which you have cited has done and what I had done. Suppose you take $\frac{d^{2}\psi}{dz^{2}}=\frac{\psi(z+\delta z)-2\psi(z)+\psi(z-\delta z)}{(\delta z)^{2}}$. Now, if $(-z_{0},z_{0})$ is my domain, then the definition for second derivative isn't valid for the boundary points, right? We will have to bring the $\psi(z\pm\delta z)$ in our $z$ domain, right? However, in the link provided by you, the same definition has been used for the boundary points. – Abhijit May 10 '17 at 20:57
• Yes, the definition is not valid for the boundary points, but you are already changing the domain from $(-\infty, \infty)$ to $[-z_0, z_0]$. You could use backward and forward differences in the two extreme nodes and then impose Dirichlet boundary conditions by removing the corresponding rows and columns associates with the first and last degree of freedom. Maybe some of the comments in this other answer can help. Knowing how $V$, $\mu$ and $\psi_0$ could also help answering. – nicoguaro May 10 '17 at 21:33
• Suppose as a start we set $V=0$, $\mu$ as some constant and $\psi_{0}=\tanh{z}$, I am not sure I fully understand the last part of what you said. I agree for the extreme nodes, we have to take forward and backward differences. But then, I do not understand the part where you say "impose Dirichlet boundary conditions by removing the corresponding rows and columns associates with the first and last degree of freedom". Maybe I am missing something obvious? Can you elucidate this part? – Abhijit May 11 '17 at 15:04
• @Abhijit, let's say that you discretize your domain in $N + 2$ points, then, we need boundary conditions, namely $\mathbf{f}(x_0) = \mathbf{f}(x_{N+1}) =0$. And you really need to solve for the interior points ($x_1, x_2, \cdots, x_N$), and for those you can form the centered difference, based on the information on neighbors. You can form your system neglecting the $x_0$ and $x_(N+1)$, solve the eigenvalue problem and then place black the removed zero values at the end ... that is more or less the idea. – nicoguaro May 13 '17 at 14:11