A finite element toolbox for the Bogoliubov-de Gennes stability analysis of Bose-Einstein condensates

We present a finite element toolbox for the computation of Bogoliubov-de Gennes modes used to assess the linear stability of stationary solutions of the Gross-Pitaevskii (GP) equation. Applications concern one (single GP equation) or two-component (a system of coupled GP equations) Bose-Einstein condensates in one, two and three dimensions of space. An implementation using the free software FreeFem++ is distributed with this paper. For the computation of the GP stationary (complex or real) solutions we use a Newton algorithm coupled with a continuation method exploring the parameter space (the chemical potential or the interaction constant). Bogoliubov-de Gennes equations are then solved using dedicated libraries for the associated eigenvalue problem. Mesh adaptivity is proved to considerably reduce the computational time for cases implying complex vortex states. Programs are validated through comparisons with known theoretical results for simple cases and numerical results reported in the literature. Program summary Program Title: FFEM_BdG_toolbox.zip CPC Library link to program files: https://doi .org /10 .17632 /dgypyc34gb .1 Licensing provisions: Apache 2.0 Programming language: FreeFem++ (v 4.12) free software (www .freefem .org) Nature of problem: The software computes Bogoliubov-de Gennes (BdG) complex modes of Bose-Einstein condensates described by the Gross-Pitaevskii (GP) equation. BdG equations are obtained by linearizing the GP equation (or the system of coupled GP equations) around a stationary solution. Obtained BdG modes are used to assess on the stability of stationary states. Solution method: Stationary states of the GP equation are obtained by a Newton algorithm. Parameter space is explored using a continuation on the chemical potential. Once the stationary (complex or real) state is captured accurately, BdG modes are computed by solving the associated eigenvalue problem with the ARPACK library. Complex eigenvalues and eigenvectors are computed and stored. The wave function is discretized by P2 (piece-wise quadratic) Galerkin triangular (in 2D) or tetrahedral (in 3D) finite elements. Mesh adaptation is implemented to reduce the computational time. Examples are given for stationary states in one-and two-component Bose-Einstein condensates.