Deflation-based Identification of Nonlinear Excitations of the 3D Gross–Pitaevskii equation

The nonlinear Schrödinger (NLS) or Gross–Pitaevskii (GP) equation ablowitz; ablowitz1; sulem; pethick; stringari; siambook is a fundamental partial differential equation that combines dispersion and nonlinearity. It has been central to a variety of areas of mathematical physics for several decades. The NLS/GP model has facilitated a universal description of a wide range of phenomena, including electric fields in optical fibers hasegawa, Langmuir waves in plasmas zakh1, freak waves in the ocean slunaev, and Bose-Einstein condensates (BECs). In the past 25 years since the experimental realization of atomic BECs, the NLS/GP model has enabled the theoretical identification and experimental observation of a wide range of coherent structures, including (but not limited to) dark djf and bright tomio solitary waves, two-dimensional vortical patterns and lattices fetter1; fetter2 as well as vortex lines and rings komineas.

The examination of three-dimensional (3D) systems has been a key frontier of recent studies in BECs. Recent theoretical advances have enabled the capturing of a number of such states mateo2014chladni. Some, especially topological ones such as skyrmions, monopoles and Alice rings ruost1; ruost2 have been of particular interest since the early exploration of BECs, while others such as knots irvine have been studied more recently. In this manuscript, we apply a powerful numerical technique called deflation farrell2015deflation; charalampidis2018computing; charalampidis2019bifurcation to identify multiple solutions of the 3D NLS/GP equation.

Many of the solutions obtained by this process are identifiable as nonlinear extensions of solutions of the linear limit of the problem, or as bifurcations therefrom. Yet other solutions are highly unexpected and are not previously known, to the best of our knowledge. Without deflation, it would be very difficult to identify these complex (literally and figuratively) topological stationary points of the infinite-dimensional energy landscape. In fact, as we increase the atom number of the system, we observe this complexity to be substantially enhanced and lead to states which, while stationary, are not straightforwardly decomposable into simpler linear or nonlinear building blocks. To further investigate the nature of the identified solutions, we compute their spectral linearization (so-called Bogoliubov-de Gennes) modes, and conduct transient simulations of prototypical unstable states to explore the dynamical behavior of their instabilities. The present work showcases, in our view, the utility and potential impact of the deflation method to complex 3D physical problems well beyond atomic BECs.

The structure of this paper is as follows. In Section II, we present the model and computational techniques employed in this work. Our numerical results on the existence, stability and selected transient simulations of nonlinear excitations are demonstrated in Section III. Section IV summarizes our findings and presents directions for future study.

The 3D NLS/GP model of interest is of the form pethick; stringari; siambook:

$\displaystyle i\psi_{t}=-\frac{1}{2}\nabla^{2}\psi+|\psi|^{2}\psi+V(\mathbf{r})\psi,$

(1)

subject to homogeneous Dirichlet conditions on the boundary of the domain $D=[-6,6]^{3}$. Here, $\psi=\psi(\mathbf{r},t)$ plays the role of the suitably normalized (see siambook for details) wavefunction, while $V$ is the external confining potential of the form $V(\mathbf{r})=\frac{1}{2}\Omega^{2}|\mathbf{r}|^{2}$, a spherically symmetric trap of strength $\Omega$, which we fix to $\Omega=1$. The boundary conditions do not affect the solutions for this choice of the trap strength since the domain is chosen large enough so that the solutions vanish well before reaching the boundary. Using the standard standing wave decomposition $\psi(\mathbf{r},t)=e^{-i\mu t}\phi(\mathbf{r})$ (where $\mu>0$ is the chemical potential), we obtain the stationary NLS/GP elliptic problem of the form:

$$F(\phi)\doteq-\frac{1}{2}\nabla^{2}\phi+|\phi|^{2}\phi+V(\mathbf{r})\phi-\mu\phi=0.$$

(2)

This equation is discretized using piecewise cubic Lagrange finite elements on a structured hexahedral grid using the Firedrake finite element library rathgeber2016. Multiple solutions to the discretized problem are sought using deflation, which we briefly describe here.

Suppose that Newton’s method has discovered an isolated root $\phi_{1}$ of $F$. Deflation constructs a new problem $G$ via

$$G(\phi)\doteq\left(\frac{1}{\|\phi-\phi_{1}\|^{2}}+1\right)F(\phi),$$

(3)

where $\|\cdot\|$ is a suitable norm, in this case the $H^{1}$ norm. The essential idea is that $\|\phi-\phi_{1}\|^{2}$ approaches $0$ faster than $F(\phi)$ does as $\phi\to\phi_{1}$, hence avoiding the convergence to $\phi_{1}$ of a Newton iteration applied to $G$. The addition of $1$ ensures that $G(\phi)\approx F(\phi)$ far from $\phi_{1}$. By applying Newton’s method to $G$, an additional root $\phi_{2}\neq\phi_{1}$ can be found, and the process repeated (by premultiplying with additional factors) until Newton’s method fails to converge from the available initial guesses.

Previous applications of deflation to the study of BECs in 2D interleaved with continuation in $\mu$, capturing solutions as they bifurcate from known ones charalampidis2018computing; charalampidis2019bifurcation. This strategy is too expensive in 3D and so a different approach is taken here. We fix $\mu=6$ and exploit the linear (low-density, i.e. $|\phi|^{2}\to 0$) limit states to furnish a large number of initial guesses for Newton’s method. The algorithm proceeds as follows. Given an initial guess, the inner loop applies Newton’s method and deflation until no more solutions are found. The outer loop iterates over the available initial guesses, and terminates when no guess yields any solutions. We emphasize that at each application of Newton’s method, all previously computed solutions are deflated, to avoid their rediscovery.

The initial guesses used were the eigenstates of the linear limit in Cartesian, cylindrical and spherical coordinates. The Cartesian eigenstates are given by

$$\ket{k,m,n}\doteq H_{k}(\sqrt{\Omega}x)H_{m}(\sqrt{\Omega}y)H_{n}(\sqrt{\Omega}z)e^{-\Omega r^{2}/2},$$

(4)

with associated energy (i.e. chemical potential) $E_{k,m,n}\doteq(k+m+n+3/2)\Omega$. The $H_{k,m,n}$ in (4) stand for the Hermite polynomials and $k,m$ and $n$ are nonnegative integers. The cylindrical eigenstates are given by

$$\ket{K,l,n}_{\textrm{cyl}}\doteq q_{K,l}(R)e^{il\theta}H_{n}(\sqrt{\Omega}z)e^{-\Omega(R^{2}+z^{2})/2},$$

(5)

with $E_{K,l,n}\doteq(2K+|l|+n+3/2)\Omega$ where $K$, $n$ are nonnegative integers, and $l=0,\pm 1,\pm 2,\dotsc$. The radial profile $q_{K,l}$ in (5) is given by $q_{K,l}\sim r^{l}L_{K}^{l}(\Omega R^{2})e^{-\Omega R^{2}/2}\Omega$ where $L_{K}^{l}$ are the associated Laguerre polynomials in $R=\sqrt{x^{2}+y^{2}}$.

Finally, the spherical eigenstates are given by $\ket{K,l,m}_{\textrm{sph}}$, where the radial part is similar but now in the spherical variable $r=\sqrt{x^{2}+y^{2}+z^{2}}$, and the angular part is described by the spherical harmonics $Y_{lm}(\theta,\phi)$ with $E_{K,l,m}=(2K+l+3/2)\Omega$. The quantum numbers $K$ and $l$ are nonnegative integers and $m=0,\pm 1,\dotsc,\pm l$. All these states with $E\leq\mu=6$ were used in the process described above.

Once a solution has been discovered, the next step is the consideration of the spectral stability of the solutions via the well-established pethick; stringari; siambook Bogoliubov-de Gennes (BdG) analysis. More specifically, we assume the following perturbation ansatz around a stationary solution $\phi^{0}$:

$$\tilde{\psi}(\mathbf{r},t)=e^{-i\mu t}\left\{\phi^{0}(\mathbf{r})+\epsilon[a(\mathbf{r})e^{i\omega t}+b^{\ast}(\mathbf{r})e^{-i\omega^{\ast}t}]\right\},$$

(6)

where $\epsilon$ is a (formal) small perturbation parameter, $\omega$ is the eigenfrequency, and $(a,b)^{\top}$ the corresponding eigenvector. After substituting Eq. (6) into the time-dependent NLS equation [cf. Eq. (1)] we obtain the following complex eigenvalue problem

$$\begin{pmatrix}A_{11}&A_{12}\\ -A_{12}^{*}&-A_{11}\end{pmatrix}\begin{pmatrix}a\\ b\end{pmatrix}=\rho\begin{pmatrix}a\\ b\end{pmatrix},$$

(7)

where $\rho=-\omega$ and the matrix elements are given by

	$\displaystyle A_{11}$	$\displaystyle=-\frac{1}{2}\nabla^{2}+2|\phi^{0}|^{2}+V(\mathbf{r})-\mu,$		(8a)
	$\displaystyle A_{12}$	$\displaystyle=\left(\phi^{0}\right)^{2}.$		(8b)

We solve the above eigenvalue problem for the eigenfrequencies $\omega$ and eigenvectors $(a,b)^{\top}$ by using a Krylov–Schur algorithm with a shift-and-invert spectral transformation stewart2002, implemented in the SLEPc library hernandez2005slepc (details about the decomposition of Eq. (7) into real and imaginary parts are presented in Appendix LABEL:app_spec). Upon convergence of the eigenvalue solver, we draw conclusions on the stability characteristics of the stationary state $\phi^{0}$ i.e., real $\omega$ implies stability (vibrations), while complex $\omega$ is associated with instability.

Finally, we explore the dynamical evolution of unstable solutions via transient numerical simulations of Eq. (1). To that end, let $\phi^{0}$ be an unstable stationary solution discovered by deflation, and $(a,b)^{\top}$ be its most unstable eigendirection normalized according to

$$\int_{D}\left(\,|a|^{2}+|b|^{2}\,\right)\,dx=1.$$

(9)

We integrate Eq. (1) forward in time until $t=50$ using the following perturbed solution as initial state

$$\psi(x,y,z,t=0)=\phi^{0}+\epsilon[a+b^{*}],$$

(10)

thus perturbing $\phi^{0}$ along its most unstable eigendirection with perturbation parameter $\epsilon$ chosen to be $0.1$. Next, let $\Delta t$ be the time step-size ($\Delta t=5\times 10^{-2}$ in this work) such that $t_{n}=n\Delta t$ and $\psi^{(n)}\doteq\psi(\mathbf{r},t_{n})$ with $n\geq 0$. Then, for given $\psi^{(n)}$ at the $n$th time step, $\psi^{(n+1)}$ is obtained (implicitly) by a modified Crank-Nicolson method delfour1981finite:

		$\displaystyle i\frac{\psi^{(n+1)}-\psi^{(n)}}{\Delta t}=$		(11)
		$\displaystyle\left(-\frac{1}{2}\nabla^{2}+V(\mathbf{r})+\frac{1}{2}(|\psi^{(n+1)}|^{2}+|\psi^{(n)}|^{2})\right)\frac{\psi^{(n+1)}+\psi^{(n)}}{2},$

where cubic finite elements are employed for the spatial discretization as before. At each time step $n$, a nonlinear problem is solved by using Newton’s method. It should be pointed out in passing that the time-marching scheme employed in this work [cf. Eq. (11)] preserves both the squared $L^{2}$ norm (i.e., atom number)

$$N(\psi)\doteq\int_{D}|\psi|^{2}\,dx,$$

(12)

and the energy of the solutions

$$\mathcal{E}(\psi)=\int_{D}\bigg{\{}{\frac{1}{4}|\nabla\psi|^{2}+\frac{1}{2}V(\mathbf{r})|\psi|^{2}+\frac{1}{4}|\psi|^{4}\bigg{\}}}dx,$$

(13)

to machine precision. We now turn to discussing the solutions obtained through the application of these numerical methods for the 3D NLS/GP problem.

We briefly describe the physical meaning of the quantum numbers for the Cartesian, cylindrical and spherical states, as they are useful in what follows. In the case of the Cartesian eigenfunctions, the quantum numbers $k,m$, and $n$ are simply the numbers of cuts along the $x$-, $y$-, and $z$-axes respectively. For instance, in Fig. 1, panel (a) represents a $\ket{0,0,2}$ Cartesian state with 2 cuts along the $z$-axis (and $\pi$ phase differences across them), while panel (b) is $\ket{1,1,0}$, bearing one planar cut along the $x$-axis, and one along the $y$-axis. Combinations of states are also possible, such as the one in panel (c) of $\ket{2,0,0}+r\ket{0,2,0}+\ket{0,0,2}$ (in the particular example of this panel $r\approx 3.39$), which forms a 2D ring along the $y$- and $z$-axes embedded in 3D space.

Vortical structures and rings can be identified in the cylindrical system of coordinates. Here, $K$ denotes the number of cylindrical (nodal) surfaces, $l$ the topological charge of the configuration and $n$ the number of planar cuts along the $z$-axis. For example, panel (d) of Fig. (1) is a so-called ring dark soliton state (extended in 3D) $\ket{1,0,0}_{\rm cyl}$ that has been recently considered in wenl, and panel (e) is the $\ket{0,2,0}_{\rm cyl}$ state, i.e., a vortex line, piercing through the BEC with topological charge $l=2$.

Finally, in the spherical representation, $K$ denotes the number of spherical (nodal) shells within the solution, $l-m$ denotes the number of planar cuts along the $z$-axis, and $m$ denotes the topological charge of vortical lines. Panel (f) is the $\ket{1,0,0}_{\rm sph}$ state corresponding to a spherical shell dark solitary wave, which is also connected with recent work wenl1.

The ground state of the system (starting at $\mu=3/2$) is known to always be spectrally and nonlinearly stable pethick; stringari. The case of the 1st excited states (e.g. dipolar states and single vortex lines) emanating from $\mu=5/2$ is interesting but reasonably well understood on the basis of corresponding 2D studies siambook, since no fundamentally novel states appear to emerge in 3D. Indicatively, LABEL:fig_2_5(a) shows a $\ket{1,0,0}$ Cartesian state with one cut along the $x$-axis whereas LABEL:fig_2_5(b) presents the $\ket{0,1,0}_{\textrm{cyl}}$ state corresponding to a single vortex line with topological charge $l=1$. The rotations of these solutions along the $x$, $y$, and $z$ axes such as the $\ket{0,1,0}$ and $\ket{0,0,1}$ Cartesian states are also obtained by deflation but are not reported. A typical example of a relevant solutio