Stability of dark solitons in a harmonic trap

Question

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?

Answer 1

As short answer, you need to:

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

2. Approximate your differential equation

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

to obtain a discrete problem

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

3. 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] $$