domain: math phys

text:
Georges Sadaka

Victor Kalt

Ionut Danaila
email: ionut.danaila@univ-rouen.fr

Frédéric Hecht

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

Keywords: Super-computer, Linux, Windows

We present a finite element toolbox for the computation of Bogoliubovde 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 1

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.

Introduction

Since their first observation [START_REF] Anderson | Observation of Bose-Einstein condensation in a dilute atomic vapor[END_REF][START_REF] Davis | Bose-Einstein condensation in a gas of sodium atoms[END_REF], Bose-Einstein condensates have become a powerful experimental framework for the study of waves and excitations in superfluids and nonlinear systems. The study of wave related structures (solitons, vortices) and their stability is an active area of research, and many efforts have been devoted to the developments of new experimental techniques for the creation and the study of new excited states. We can mention the use of rotation [START_REF] Madison | Vortices in a stirred Bose-Einstein condensate[END_REF][START_REF] Haljan | Driving Bose-Einstein condensate vorticity with a rotating normal cloud[END_REF], imprinting techniques manipulating the phase of the wave function [START_REF] Leanhardt | Imprinting vortices in a Bose-Einstein condensate using topological phases[END_REF][START_REF] Becker | Oscillations and interactions of dark and dark-bright solitons in Bose-Einstein condensates[END_REF], counterflows techniques [START_REF] Yan | Multiple darkbright solitons in atomic Bose-Einstein condensates[END_REF], the use of anisotropic potentials [START_REF] Theocharis | Multiple atomic dark solitons in cigar-shaped Bose-Einstein condensates[END_REF] or multicomponent BEC settings (Wang and Kevrekidis, 2017). A large variety of wave or vortex-related states could be thus obtained. Basic examples refer to single vortex lines (with I-, U-or S-shape) in rotating BEC [START_REF] Aftalion | Three-dimensional vortex configurations in a rotating Bose-Einstein condensate[END_REF], vortex rings and one-dimensional solitons. More complex states with multiple vortex rings, vortex stars, hopfions and solitons can be created [START_REF] Crasovan | Three-dimensional parallel vortex rings in Bose-Einstein condensates[END_REF][START_REF] Bisset | Robust vortex lines, vortex rings, and hopfions in three-dimensional Bose-Einstein condensates[END_REF]Wang et al., 2017). In multicomponent BECs, dark-bright [START_REF] Charalampidis | Bifurcation analysis of stationary solutions of two-dimensional coupled Gross-Pitaevskii equations using deflated continuation[END_REF] and dark-antidark [START_REF] Danaila | Vector dark-antidark solitary waves in multicomponent Bose-Einstein condensates[END_REF] states can be also obtained.

The main application of the programs presented in this paper is the study of the stability of such solitary waves or vortex states that are theoretically or numerically found as stationary solutions to the Gross-Pitaevskii (GP) equation [START_REF] Gross | Structure of a quantized vortex in boson systems[END_REF][START_REF] Pitaevskii | Vortex lines in an imperfect Bose gas[END_REF]. The linearization of the GP equation around a given stationary solution results in the Bogoliubov-de Gennes (BdG) system of equations [START_REF] Pitaevskii | Bose-Einstein condensation[END_REF]. Solving the BdG eigenvalue problem provides linear modes, for which their stability could be studied. The present toolbox thus contains two distinct parts: the computation of stationary solution of the GP equation (or a system of coupled GP equations) and the computation of complex eigenvalues and modes for the associated BdG system.

Concerning the computation of stationary states of the GP equation, a large variety of discretization methods were suggested in the open literature: spectral methods [START_REF] Dion | Ground state of the time-independent Gross-Pitaevskii equation[END_REF][START_REF] Caliari | GSGPEs: A Matlab code for computing the ground state of systems of Gross-Pitaevskii equations[END_REF][START_REF] Antoine | a Matlab toolbox to solve Gross-Pitaevskii equations I: Computation of stationary solutions[END_REF], finite-elements [START_REF] Marojević | ATUS-PRO: A FEM-based solver for the time-dependent and stationary Gross-Pitaevskii equation[END_REF][START_REF] Vergez | A finite-element toolbox for the stationary Gross-Pitaevskii equation with rotation[END_REF] or finite-differences [START_REF] Muruganandam | Fortran programs for the time-dependent Gross-Pitaevskii equation in a fully anisotropic trap[END_REF][START_REF] Vudragović | C programs for solving the time-dependent Gross-Pitaevskii equation in a fully anisotropic trap[END_REF][START_REF] Caplan | NLSEmagic: Nonlinear Schrödinger equation multi-dimensional Matlab-based GPU-accelerated integrators using compact high-order schemes[END_REF][START_REF] Kong | LOD-ms for Gross-Pitaevskii equation in Bose-Einstein condensates[END_REF][START_REF] Hohenester | OCTBEC a Matlab toolbox for optimal quantum control of Bose-Einstein condensates[END_REF][START_REF] Kumar | C and Fortran OpenMP programs for rotating Bose-Einstein condensates[END_REF]. Programs written in Fortran [START_REF] Dion | Ground state of the time-independent Gross-Pitaevskii equation[END_REF][START_REF] Muruganandam | Fortran programs for the time-dependent Gross-Pitaevskii equation in a fully anisotropic trap[END_REF], C [START_REF] Vudragović | C programs for solving the time-dependent Gross-Pitaevskii equation in a fully anisotropic trap[END_REF][START_REF] Caplan | NLSEmagic: Nonlinear Schrödinger equation multi-dimensional Matlab-based GPU-accelerated integrators using compact high-order schemes[END_REF], Matlab [START_REF] Caliari | GSGPEs: A Matlab code for computing the ground state of systems of Gross-Pitaevskii equations[END_REF][START_REF] Antoine | a Matlab toolbox to solve Gross-Pitaevskii equations I: Computation of stationary solutions[END_REF][START_REF] Caplan | NLSEmagic: Nonlinear Schrödinger equation multi-dimensional Matlab-based GPU-accelerated integrators using compact high-order schemes[END_REF][START_REF] Hohenester | OCTBEC a Matlab toolbox for optimal quantum control of Bose-Einstein condensates[END_REF], FreeFem++ [START_REF] Vergez | A finite-element toolbox for the stationary Gross-Pitaevskii equation with rotation[END_REF] or C and Fortran with OpenMP (Kishor [START_REF] Kumar | C and Fortran OpenMP programs for rotating Bose-Einstein condensates[END_REF] have been shared. All these works use a common approach to the problem, which is to find the stationary states of the GP equation under the constraint of the total mass conservation (i. e. the L2-norm of the wave function is fixed). When studying the stability of excitations in the BdG framework, another approach is commonly adopted: the chemical potential is used as a convenient parameter to explore all possible states (and bifurcations) and thus the total mass varies from one state to another. This second approach has been already applied with finite elements [START_REF] Charalampidis | Bifurcation analysis of stationary solutions of two-dimensional coupled Gross-Pitaevskii equations using deflated continuation[END_REF][START_REF] Danaila | Vector dark-antidark solitary waves in multicomponent Bose-Einstein condensates[END_REF][START_REF] Boullé | Deflationbased identification of nonlinear excitations of the three-dimensional Gross-Pitaevskii equation[END_REF][START_REF] Carretero-González | Vortex nucleation in a dissipative variant of the nonlinear Schrödinger equation under rotation[END_REF], finite differences [START_REF] Bisset | Robust vortex lines, vortex rings, and hopfions in three-dimensional Bose-Einstein condensates[END_REF]Wang et al., 2017;[START_REF] Charalampidis | Bifurcation analysis of stationary solutions of two-dimensional coupled Gross-Pitaevskii equations using deflated continuation[END_REF] or spectral methods (Wang et al., 2017) and will be used in this work. To solve the BdG eigenvalue problem, specialized libraries are generally used: ARPACK [START_REF] Lehoucq | ARPACK users' guide: solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods[END_REF], SLEPc [START_REF] Hernandez | SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems[END_REF] or FEAST [START_REF] Polizzi | Density-matrix-based algorithm for solving eigenvalue problems[END_REF]. Such libraries offer the flexibility to be easily interfaced with different types of discretization, since only final matrices of the eigenvalue problem are needed. A mathematical study of the BdG equation with numerical comparisons between spectral and finite difference discretizations has been recently reported in [START_REF] Gao | Numerical methods for Bogoliubov-de Gennes excitations of Bose-Einstein condensates[END_REF].

To the best of our knowledge, FACt (Fluctuations in Atomic Condensates, [START_REF] Roy | FACt: FORTRAN toolbox for calculating fluctuations in atomic condensates[END_REF] is the only publicly accessible code for the BdG problem. It considers thermal excitations of the BECs at non-zero temperatures by solving the BdG equation in two component BECs with a pseudo-spectral method. The present FreeFem++ toolbox uses a different (PDE-oriented) formulation of the problem and has the advantage to hide all technicalities related to the implementation of the finite-element method and the interface with eigenvalue libraries (such as ARPACK) [START_REF] Hecht | New developments in Freefem++[END_REF]. The user can thus focus on the physical and mathematical model, and eventually on the numerical algorithm solving the problem. The high level programming language offered in FreeFem++ and the syntax close to mathematical formulations make the implementation of numerical algorithms very easy. Finite elements algorithms were successfully used to solve the GP equation [START_REF] Vergez | A finite-element toolbox for the stationary Gross-Pitaevskii equation with rotation[END_REF][START_REF] Danaila | A finite element method with mesh adaptivity for computing vortex states in fast-rotating Bose-Einstein condensates[END_REF] or the BdG problem [START_REF] Danaila | Vector dark-antidark solitary waves in multicomponent Bose-Einstein condensates[END_REF] and recently to identify vortices in a quantum field [START_REF] Kalt | Identification of vortices in quantum fluids: Finite element algorithms and programs[END_REF]. Another advantage of the present toolbox is to use mesh adaptivity to reduce the mesh size and the computational time. Solving the BdG problem for complex cases becomes thus possible using personal computers.

The structure of the paper is as follows. In Sect. 2, we introduce the GP and BdG models. Sections 3 and 4 describe the numerical methods used for the computation the stationary states and BdG modes. We present various benchmarks used for the validation of our codes in Sects. 5 and 6. The architecture of the programs and a description of parameter and output files are given in Sect. 7. Finally, we summarize the main features of the toolbox and present some possible extensions in Sect. 8.

The Gross-Pitaevskii model and Bogoliubovde Gennes equations

The Gross-Pitaevskii equation

In the zero temperature limit, the Gross-Pitaevskii equation describes the time-evolution of the complex-valued macroscopic wave function ψ :

D → C, with D ∈ R d the domain of the d-dimensional condensate (d = 1, 2, 3): i ∂ψ ∂t = - 2 2m ∇ 2 ψ + V trap ψ + g|ψ| 2 ψ, ( 1
)
where V trap (x) is the external trapping potential, the reduced Planck constant and m the atomic mass. The nonlinear term models the interaction between atoms, and for d = 3 (3D condensate) g = 4π 2 as m , with a s the scattering length. For d = 1 or 2 (1D or 2D condensates), the nonlinear interaction constant g is specified according to the dimension reduction [START_REF] Frantzeskakis | Dark solitons in atomic Bose-Einstein condensates: from theory to experiments[END_REF][START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF]. We consider the case of harmonic trapping potentials:

V trap (x) = m 2 ω 2 x x 2 + ω 2 y y 2 + ω 2 z z 2 , ( 2
)
where ω x , ω x , ω x are trapping frequencies. The atomic density n(x) = |ψ(x)| 2 vanishes outside the condensate due to the trapping, which implies that homogeneous Dirichlet boundary conditions should be imposed for the wave function (ψ = 0 on ∂D). The corresponding GP energy is:

E(ψ) = D 2 2m |∇ψ(x, t)| 2 + V trap (x)|ψ(x, t)| 2 + g 2 |ψ(x, t)| 4 dx, (3)
and the total number of atoms:

N (ψ) = D ψψ dx = D |ψ| 2 dx, ( 4
)
where ψ denotes the complex conjugate.

Stationary solutions to the GP equation ( 1) are obtained by imposing the form ψ(x, t) = φ(x)e -i µt , ( 5)

with µ the chemical potential. The stationary wave function φ is then solution of the stationary GP equation:

-

2 2m ∇ 2 φ + V trap φ + g|φ| 2 φ = µφ, ( 6
)
which is a nonlinear eigenvalue problem. Note that from (5) we infer that |ψ| 2 = |φ| 2 and thus N (φ) = N (ψ) and E(ψ) = E(φ). The chemical potential is related to the number of atoms by the relation:

µ = 1 N (φ) E(φ) + g 2 D |φ| 4 dx . (7)
In this work, we compute stationary solutions for fixed values of the chemical potential µ. Branches of solutions are followed by continuation on µ. To catch such branches, two limits associated to the value of µ can be considered as initial condition. In the case of low density (corresponding to a small number of particles), the nonlinear term in (1) can be neglected to obtain the linear GP equation. In the case of an harmonic potential (2), ψ can then be described as an eigenstate of the quantum harmonic oscillator. Using a separation of variables, these eigenstates can be written as a product of Hermite and Laguerre polynomials or spherical harmonics, depending on the dimension and the coordinate system. As an example, for the 2D BEC with trapping potential (2) with

ω x = ω y = ω ⊥ , the eigenstates formulated in Cartesian coordinates are |k, l ∝ H k ( √ ω ⊥ x)H l ( √ ω ⊥ y)e -1 2 ω ⊥ (x 2 +y 2 )
, where H k , H l are Hermite polynomials; k and l index the eigenstates and correspond to the number of cuts in the condensate along the x and y axes respectively. Solutions in the linear limit corresponding to various exited states of the condensate have been analysed in many studies (e. g. [START_REF] Crasovan | Three-dimensional parallel vortex rings in Bose-Einstein condensates[END_REF][START_REF] Boullé | Deflationbased identification of nonlinear excitations of the three-dimensional Gross-Pitaevskii equation[END_REF].

The other limit is the Thomas-Fermi limit, associated to large values of µ. In this case, the kinetic energy becomes negligible when compared to the nonlinear term. The stationary GP equation ( 6) reduces to:

µφ = V trap φ + g|φ| 2 φ, ( 8
)
which gives an analytical expression for the atomic density:

n TF = |ψ TF | 2 = 1 g (µ -V trap ) + . (9)

The Bogoliubov-de Gennes equation

The Bogoliubov-de Gennes model is based on the linearisation of (1) assuming that:

ψ(x, t) = (φ(x) + δφ(x, t))e -i µt , ( 10
)
where φ(x) is a stationary state satisfying Eq. ( 6) and δφ a small perturbation. Inserting (10) in (1), we obtain, after neglecting second order terms in δφ, an evolution equation for the perturbation δφ:

i ∂δφ ∂t = Hδφ -µδφ + 2g|φ| 2 δφ + gφ 2 δφ, ( 11
)
where we denoted by H ≡ -2 2m ∇ 2 + V trap the linear part of the Hamiltonian. Considering perturbations of the form

δφ(x, t) = A(x)e -iωt + B(x)e iωt , (12)
we obtain, after separating terms in e -iωt and e iωt , the Bogoliubov-de Gennes (BdG) system of equations [START_REF] Pitaevskii | Bose-Einstein condensation[END_REF][START_REF] Castin | Bose-einstein condensates in atomic gases: Simple theoretical results[END_REF]:

H -µ + 2g|φ| 2 gφ 2 -gφ 2 -(H -µ + 2g|φ| 2 ) A B = ω A B . ( 13
)
Note that the BdG equation ( 13) is a linear eigenvalue problem, since φ is fixed and H is a linear real operator. The present toolbox computes, for a given complex stationary state φ ∈ C, solutions (ω, A, B) to the BdG equation (13), with ω denoting complex eigenvalues and (A, B) complex eigenvectors.

The following properties of the BdG eigenvalue problem can be obtained by elementary manipulations and will be useful to check the accuracy of calculations:

1. If (ω, A, B) is solution to (13), then (-ω, B, A) is also a solution. This property is obtained by taking the conjugate of (13).

2. ω = 0 is always an eigenvalue (the zero-energy mode). It can be checked that the full solution is (0, αφ, -αφ), α ∈ C [START_REF] Pitaevskii | Bose-Einstein condensation[END_REF] and represents following (12) a time invariant (small) excitation δφ = (α -α)φ. This corresponds in (10) to a gauge transformation and, consequently, it does not add any physical excitation to the system.

3. ω is also an eigenvalue. If φ ∈ R, it is easy to see from ( 13) that (ω, A, B) is also a solution. For the general case of φ ∈ C this property also holds and it can be proved using the Hamiltonian nature of the problem [START_REF] Castin | Bose-einstein condensates in atomic gases: Simple theoretical results[END_REF].

4. If we multiply the first equation of the system (13) by A and the second by B, integrate over the domain D and then sum the two equations, we obtain that:

δE = ω D |A| 2 -|B| 2 dx ∈ R, (14)
which is generally presented in the literature in the equivalent form [START_REF] Pitaevskii | Bose-Einstein condensation[END_REF][START_REF] Gao | Numerical methods for Bogoliubov-de Gennes excitations of Bose-Einstein condensates[END_REF]Kevrekidis et al., 2008):

(ω -ω) D |A| 2 -|B| 2 dx = 0. ( 15
)
From ( 14) or (15) we can draw two main observations that are important to interpret the results of the BdG analysis:

1. If the BdG modes are normalized such that D (|A| 2 -|B| 2 ) dx = 0, then we infer from (15) that only real eigenvalues are possible. These correspond to elementary excitations. A mathematical study of the properties of the BdG modes when ω and φ are real is offered in [START_REF] Gao | Numerical methods for Bogoliubov-de Gennes excitations of Bose-Einstein condensates[END_REF]. Moreover, the quantity δE in ( 14) represents the energy difference between the stationary ψ and the perturbed state ψ + δψ [START_REF] Pitaevskii | Bose-Einstein condensation[END_REF]. The Krein signature K was introduced as the sign of the energy difference, K = sign(δE) [START_REF] Chernyavsky | Krein signature in Hamiltonian and PT-symmetric systems[END_REF]. If K > 0 for all modes, then φ is the global minimum of the energy, i. e. the ground state. On the contrary, if there exists a mode with K < 0, then the excitation reduces the energy of the system and the stationary state is thus energetically unstable, i. e. excited state (or local minimum of the energy).

2. If D (|A| 2 -|B| 2 ) dx = 0, complex eigenvalues ω = ω r +iω i are possible.

If ω i = 0, then the BdG mode is dynamically unstable.

Two-component BECs

Mixtures of BECs have been experimentally created, either by considering different atomic species or by using hyperfine states of a single isotope. We consider two-component BECs described by the following system of two coupled GP equations for wave functions ψ 1 and ψ 2 :

           i ∂ψ 1 ∂t = - 2 2m ∇ 2 + V trap + g 11 |ψ 1 | 2 + g 12 |ψ 2 | 2 ψ 1 , i ∂ψ 2 ∂t = - 2 2m ∇ 2 + V trap + g 21 |ψ 1 | 2 + g 22 |ψ 2 | 2 ψ 2 . ( 16
)
Coefficients g 11 and g 22 represent interactions between atoms of the same species while g 12 and g 21 describe interactions between different species. The total energy is the sum of the GP energy of each component:

E(ψ 1 , ψ 2 ) = D 2 i=1   2 2m |∇ψ i | 2 + V trap |ψ i | 2 + 1 2 2 j=1 g ij |ψ i | 2 |ψ j | 2   dx. ( 17
)
Similarly to the one component case, stationary states are sought as ψ 1 = φ 1 e -i µ 1 t and ψ 2 = φ 2 e -i µ 2 t , with chemical potentials µ 1 and µ 2 . We obtain the following system of equations:

           µ 1 φ 1 = - 2 2m ∇ 2 + V trap + g 11 |φ 1 | 2 + g 12 |φ 2 | 2 φ 1 , µ 2 φ 2 = - 2 2m ∇ 2 + V trap + g 21 |φ 1 | 2 + g 22 |φ 2 | 2 φ 2 . ( 18
)
For the linear stability analysis we consider perturbations of the form:

   δφ 1 (x) = A(x)e -iωt + B(x)e iωt , δφ 2 (x) = C(x)e -iωt + D(x)e iωt , ( 19
)
and obtain the BdG equations for the two-component case:

M      A B C D      = ω      A B C D      , ( 20
)
where the matrix M can be presented in the form:

M =       M 11 g 11 φ 2 1 g 12 φ 1 φ 2 g 12 φ 1 φ 2 -g 11 φ 1 2 M 22 -g 12 φ 1 φ 2 -g 12 φ 1 φ 2 g 21 φ 1 φ 2 g 21 φ 1 φ 2 M 33 g 22 φ 2 2 -g 21 φ 1 φ 2 -g 21 φ 1 φ 2 -g 22 φ 2 2 M 44       , (21)
with:

             M 11 = H -µ 1 + 2g 11 |φ 1 | 2 + g 12 |φ 2 | 2 , M 22 = -M 11 , M 33 = H -µ 2 + g 21 |φ 1 | 2 + 2g 22 |φ 2 | 2 , M 44 = -M 33 . ( 22
)

Scaling

Various forms of scaling are used in the literature (see [START_REF] Danaila | From Classical to Quantum Turbulence: Basic Concepts and Models[END_REF] for a unified form of the GP scaling). To allow one to switch between different forms, we first introduce a reference (trapping) frequency ω s which will define a time scale t s and a length scale x s (the corresponding harmonic oscillator length):

t s = 1 ω s , x s = mω s . ( 23
)
We then introduce a reference value ψ s for the wave function and scale variables as:

x → x x s , t → t t s , ψ → ψ ψ s . ( 24
)
The same scaling is used for the stationary state φ. The dimensionless form of the time-dependent stationary GP equation (1) becomes:

i ∂ψ ∂t = - 1 2 ∇ 2 ψ + C trap ψ + β|ψ| 2 ψ, ( 25
)
with

C trap (x) = 1 ω s V trap (x), β = gψ 2 s ω s . ( 26
)
Note that the coefficient 1/2 in front of the Laplacian in Eq. ( 25) comes from the choice (23) for the reference length, since all terms were divided by the quantity (of the dimension of an energy):

ω s = mx 2 s ω 2 s = 2 mx 2 s . ( 27
)
From ( 2) we infer that the non-dimensional trapping potential takes the form:

C trap (x, y, z) = 1 2 ω 2 x x 2 + ω 2 y y 2 + ω 2 z z 2 , where ω x,y,z → ω x,y,z ω s . ( 28
)
Similarly, the stationary GP equation ( 6) becomes:

- 1 2 ∇ 2 φ + C trap φ + β|φ| 2 = µφ, where µ → µ ω s . ( 29
)
Finally, the BdG system of equations takes the non-dimensional form:

H -µ + 2β|φ| 2 gφ 2 -gφ 2 -(H -µ + 2β|φ| 2 ) A B = ω A B , (30)
where

H ≡ -1 2 ∇ 2 + C trap is dimensionless, ω → ω/ω s and A → A/ψ s , B → B/ψ s .
In the two component case, the wave functions are scaled as

ψ 1 → ψ 1 ψ s,1
and ψ 2 → ψ 2 ψ s,2 , and the system (18) for the stationary state becomes:

         µ 1 φ 1 = - 1 2 ∇ 2 + C trap + β 11 |φ 1 | 2 + β 12 |φ 2 | 2 φ 1 , µ 2 φ 2 = - 1 2 ∇ 2 + C trap + β 21 |φ 1 | 2 + β 22 |φ 2 | 2 φ 2 , (31)
where

β ij → g ij ωs ψ 2 s,j , µ i → µ i ωs .
It follows that the non-dimensional form of the BdG system for the two-component BEC is obtained from ( 20)-( 22), by replacing coefficients g ij with β ij and using the non-dimensional form of H.

Computing stationary solutions of the GP equation

Newton method for a single component BEC

Stationary solutions of Eq. ( 29) are computed using a Newton method.

Considering that φ = φ r + iφ i , we obtain the following system of equations after separating real and imaginary parts:

       - 1 2 ∇ 2 φ r + C trap φ r + βf (φ r , φ i )φ r -µφ r = 0, - 1 2 ∇ 2 φ i + C trap φ i + βf (φ r , φ i )φ i -µφ i = 0, (32)
where f (φ r , φ i ) implements the non-linear (interaction) term. Here

f (φ r , φ i ) = |φ| 2 = φ 2 r + φ 2 i
, but the method is described (and programmed) for a general expression of f that can be easily changed in the toolbox.

We consider homogeneous Dirichlet boundary conditions for φ r and φ i , i. e. φ r = φ i = 0 on ∂D, and set the classical Sobolev spaces V = H 1 0 (D) for φ r and φ i . The weak formulation of ( 32) can be thus written as follows: find

(φ r , φ i ) ∈ V × V = V 2 , such that for all test functions (v r , v i ) ∈ V 2 :        F r (φ r , φ i , v r ) = D (C trap -µ)φ r v r + D 1 2 ∇φ r • ∇v r + D βf (φ r , φ i )φ r v r = 0, F i (φ r , φ i , v i ) = D (C trap -µ)φ i v i + D 1 2 ∇φ i • ∇v i + D βf (φ r , φ i )φ i v i = 0.
(33) Starting from an initial guess (φ 0 r , φ 0 i ), solution increments

q = φ k r -φ k+1 r , s = φ k i -φ k+1 i , k ≥ 0, ( 34
)
are computed using the Newton algorithm:

       ∂F r ∂φ r φr=φ k r ,φ i =φ k i ∂F r ∂φ i φr=φ k r ,φ i =φ k i ∂F i ∂φ r φr=φ k r ,φ i =φ k i ∂F i ∂φ i φr=φ k r ,φ i =φ k i        q s = F r (φ k r , φ k i , v r ) F i (φ k r , φ k i , v i ) , (35)
with the corresponding weak formulation:

                           D (C trap -µ)qv r + D 1 2 ∇q • ∇v r + D β ∂f ∂φ r (φ k r , φ k i ) φ k r q + ∂f ∂φ i (φ k r , φ k i ) φ k r s + f (φ k r , φ k i )q v r = D (C trap -µ)φ k r v r + D 1 2 ∇φ k r • ∇v r + D βf (φ k r , φ k i )φ k r v r , D (C trap -µ)sv i + D 1 2 ∇s • ∇v i + D β ∂f ∂φ r (φ k r , φ k i ) φ k i q + ∂f ∂φ i (φ k r , φ k i ) φ k i s + f (φ k r , φ k i )s v i = D (C trap -µ)φ k i v i + D 1 2 ∇φ k i • ∇v i + D βf (φ k r , φ k i )φ k i v i . ( 36
)
Note that the metalanguage used in FreeFem++ enables the implementation of Eqs. (36) in a form very similar to mathematical formulae, which is appreciable to rapidly build bug-free numerical codes.

Newton method for a two-component BEC

For a two-component BEC, we solve the following system, obtained from (31) after separating real and imaginary parts for φ 1 = φ 1r +iφ 1i and φ 2 = φ 2r +iφ 2i :

                           - 1 2 ∇ 2 φ 1r + (C trap -µ 1 )φ 1r + β 11 f (φ 1r , φ 1i )φ 1r + β 12 f (φ 2r , φ 2i )φ 1r = 0, - 1 2 ∇ 2 φ 1i + (C trap -µ 1 )φ 1i + β 11 f (φ 1r , φ 1i )φ 1i + β 12 f (φ 2r , φ 2i )φ 1i = 0, - 1 2 ∇ 2 φ 2r + (C trap -µ 2 )φ 2r + β 21 f (φ 1r , φ 1i )φ 2r + β 22 f (φ 2r , φ 2i )φ 2r = 0, - 1 2 ∇ 2 φ 2i + (C trap -µ 2 )φ 2i + β 21 f (φ 1r , φ 1i )φ 2i + β 22 f (φ 2r , φ 2i )φ 2i = 0.
(37) We consider again homogeneous Dirichlet boundary conditions, i. e. φ 1r = φ 1i = φ 2r = φ 2i = 0 on ∂D. The weak formulation of (37) can be written as follows: find (φ 1r , φ 1i , φ 2r , φ 2i ) ∈ V 4 , such that for all test functions

(v 1r , v 1i , v 2r , v 2i ) ∈ V 4 :                          F 1r = D (C trap -µ 1 )φ 1r v 1r + D 1 2 ∇φ 1r • ∇v 1r + D (β 11 f (φ 1r , φ 1i ) + β 12 f (φ 2r , φ 2i ))φ 1r v 1r = 0, F 1i = D (C trap -µ 1 )φ 1i v 1i + D 1 2 ∇φ 1i • ∇v 1i + D (β 11 f (φ 1r , φ 1i ) + β 12 f (φ 2r , φ 2i ))φ 1i v 1i = 0, F 2r = D (C trap -µ 2 )φ 2r v 2r + D 1 2 ∇φ 2r • ∇v 2r + D (β 21 f (φ 1r , φ 1i ) + β 22 f (φ 2r , φ 2i ))φ 2r v 2r = 0, F 2i = D (C trap -µ 2 )φ 2i v 2i + D 1 2 ∇φ 2i • ∇v 2i + D (β 21 f (φ 1r , φ 1i ) + β 22 f (φ 2r , φ 2i ))φ 2i v 2i = 0. ( 38
)
The Newton step for increments

q 1 = φ k 1r -φ k+1 1r , s 1 = φ k 1i -φ k+1 1i , q 2 = φ k 2r -φ k+1 2r , s 2 = φ k 2i -φ k+1 2i , ( 39
)
13 consists then in solving the following four equations:

D (C trap -µ 1 )q 1 v 1r + D 1 2 ∇q 1 • ∇v 1r + D (β 11 f (φ k 1r , φ k 1i ) + β 12 f (φ k 2r , φ k 2i ))q 1 v 1r + D β 11 ∂f ∂φ r (φ k 1r , φ k 1i )φ k 1r q 1 + ∂f ∂φ i (φ k 1r , φ k 1i )φ k 1r s 1 v 1r + D β 12 ∂f ∂φ r (φ k 2r , φ k 2i )φ k 1r q 2 + ∂f ∂φ i (φ k 2r , φ k 2i )φ k 1r s 2 v 1r = D (C trap -µ 1 )φ k 1r v 1r + D 1 2 ∇φ k 1r • ∇v 1r + D (β 11 f (φ k 1r , φ k 1i ) + β 12 f (φ k 2r , φ k 2i ))φ k 1r v 1r , (40) D (C trap -µ 1 )s 1 v 1i + D 1 2 ∇s 1 • ∇v 1i + D (β 11 f (φ k 1r , φ k 1i ) + β 12 f (φ k 2r , φ k 2i ))s 1 v 1i + D β 11 ∂f ∂φ r (φ k 1r , φ k 1i )φ k 1i q 1 + ∂f ∂φ i (φ k 1r , φ k 1i )φ k 1i s 1 v 1i + D β 12 ∂f ∂φ r (φ k 2r , φ k 2i )φ k 1i q 2 + ∂f ∂φ i (φ k 2r , φ k 2i )φ k 1i s 2 v 1i = D (C trap -µ 1 )φ k 1i v 1i + D 1 2 ∇φ k 1i • ∇v 1i + D (β 11 f (φ k 1r , φ k 1i ) + β 12 f (φ k 2r , φ k 2i ))φ k 1i v 1i , (41) D (C trap -µ 2 )q 2 v 2r + D 1 2 ∇q 2 • ∇v 2r + D (β 22 f (φ k 2r , φ k 2i ) + β 21 f (φ k 1r , φ k 1i ))q 2 v 2r + D β 21 ∂f ∂φ r (φ k 1r , φ k 1i )φ k 2r q 1 + ∂f ∂φ i (φ k 1r , φ k 1i )φ k 2r s 1 v 2r + D β 22 ∂f ∂φ r (φ k 2r , φ k 2i )φ k 2r q 2 + ∂f ∂φ i (φ k 2r , φ k 2i )φ k 2r s 2 v 2r = D (C trap -µ 2 )φ k 2r v 2r + D 1 2 ∇φ k 2r • ∇v 2r + D (β 21 f (φ k 1r , φ k 1i ) + β 22 f (φ k 2r , φ k 2i ))φ k 2r v 2r , (42) D (C trap -µ 2 )s 2 v 2i + D 1 2 ∇s 2 • ∇v 2i + D (β 22 f (φ k 2r , φ k 2i ) + β 21 f (φ k 1r , φ k 1i ))s 2 v 2i + D β 21 ∂f ∂φ r (φ k 1r , φ k 1i )φ k 2i q 1 + ∂f ∂φ i (φ k 1r , φ k 1i )φ k 2i s 1 v 2i + D β 22 ∂f ∂φ r (φ k 2r , φ k 2i )φ k 2i q 2 + ∂f ∂φ i (φ k 2r , φ k 2i )φ k 2i s 2 v 2i = D (C trap -µ 2 )φ k 2i v 2i + D 1 2 ∇φ k 2i • ∇v 2i + D (β 21 f (φ k 1r , φ k 1i ) + β 22 f (φ k 2r , φ k 2i ))φ k 2i v 2i . (43)
Again, the implementation of Eqs. ( 40)-( 43) with FreeFem++ is very similar to the mathematical formulation.

Finite element implementation with FreeFem++

The algorithms presented below are implemented using the free software FreeFem++ [START_REF] Hecht | New developments in Freefem++[END_REF]. We illustrate in this section the main principles of programming used in building the toolbox and the numerical settings for the BdG problem.

One of the main advantages offered by FreeFem++ is to program cumbersome formulae in a compact form, close to the mathematical formulation. For example, the system (36) is implemented in a Macro (precisely BdG 1comp/A macro/Macro problem.edp) in which integral terms are easy to identify: NewMacro problemGP macro f ( ur , ui ) ( ur ˆ2 + ui ˆ2) // macro dfdur ( ur , ui ) (2.* ur ) // macro dfdui ( ur , ui ) (2.* ui ) // varf vGP ([ q , s ] , [ vr , vi ]) = intN ( Th , qforder = ord ) (( Ctrap -mu ) * q * vr + .5* grad ( q ) '* grad ( vr ) + ( Ctrap -mu ) * s * vi + .5* grad ( s ) '* grad ( vi ) + beta * ( f ( phir , phii ) * q * vr + f ( phir , phii ) * s * vi ) + beta * phir * vr *( dfdur ( phir , phii ) * q + dfdui ( phir , phii ) * s ) + beta * phii * vi *( dfdur ( phir , phii ) * q + dfdui ( phir , phii ) * s ) ) + intN ( Th , qforder = ord ) (( Ctrap -mu ) * phir * vr + .5* grad ( phir ) '* grad ( vr ) + ( Ctrap -mu ) * phii * vi + .5* grad ( phii ) '* grad ( vi ) + beta * f ( phir , phii ) * ( phir * vr + phii * vi ) ) BCGP ; EndMacro Another advantage of this formulation is that it can be used for any dimension (d = 1, 2 or 3) and any available type of finite elements, by simply declaring these values in the files defining the computational case. For example, for the 1D dark-soliton test case (file BdG 1comp/INIT/1D DS.inc):

macro dimension 1 // macro FEchoice P2 //
These choices are then used in the main programs to define the finiteelement spaces. For example, in FFEM GP 1c 1D 2D 3D.edp: func Pk = [ FEchoice , FEchoice ]; ... meshN Th ; // Local mesh fespace Wh ( Th , FEchoice ) ; fespace Whk ( Th , Pk ) ; ... Wh < complex > phi ; // Wavefunction Whk [q , s ] , [ phir , phii ];

For all programs in this toolbox, we use P 2 (piece-wise quadratic) finite elements.

FreeFem++ also offers a fast mesh generator for 1D, 2D or 3D configurations. The mesh (generically identified as Th) is made of segments in 1D, triangles in 2D and tetrahedrons in 3D. The initial solution is built specifically for each case as an approximation of the state we want to study (see Sects. 5 and 6). Newton iterations are stopped when one of two following criteria is satisfied:

q s ∞ < q , F r F i 2 < F . ( 44
)
In practice, we use q = 10 -8 and F = 10 -16 and in all considered cases both criteria are satisfied simultaneously. To achieve convergence in the Newton algorithm, the choice of the solver for the linear system resulting from (32) or ( 37) is very important. For 1D and 2D problems we solve the system with a direct LU method using the library MUMPS. For 3D problems, we use a GMRES method, preconditioned by an incomplete LU factorization. Branches of stationary solutions are followed by a continuation method on the parameter µ 0 ≤ µ ≤ µ f . In practice, we start from a value µ 0 for which the initial condition is sufficiently close to the stationary state and use this converged state as an initial guess for the Newton method with chemical potential µ 0 + δµ. The process is repeated until µ f is reached. This is especially useful when following states from the linear limit to high values of µ. For the two-component case, the continuation is done first on µ 1 and µ 2 and then on the inter-component interactions β 12 and β 21 .

An important tool in FreeFem++ is mesh adaptation, that considerably helps in reducing the computational time while keeping a high degree of accuracy. The mesh is adapted in 2D using the standard adaptmesh command of FreeFem++ which creates a new mesh adapted to the Hessian of the solution. In 3D, the adaptation is done through the libraries mshmet and mmg [START_REF] Dapogny | Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems[END_REF] which are directly linked to FreeFem++ . When using continuation, we adapt the mesh for different values of µ. Mesh adaptation is mandatory for the complicated test cases, especially in 3D or for the two-component cases: using a refined mesh for the entire domain would lead to a large memory consumption and an excessively long computational time.

Solving the BdG equations

The BdG problem (30) is solved using the ARPACK library [START_REF] Lehoucq | ARPACK users' guide: solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods[END_REF]. It is directly interfaced with FreeFem++ and uses an Arnoldi method to compute the eigenvalues and eigenvectors of a given matrix. We use the following weak formulation corresponding to (30):

       D 1 2 ∇A • ∇v 1 + D (C trap -µ)Av 1 + D 2β|φ| 2 Av 1 + D βφ 2 Bv 1 = ω D Av 1 , - D 1 2 ∇B • ∇v 2 - D (C trap + µ)Bv 2 - D 2β|φ| 2 Bv 2 - D βφ 2 Av 2 = ω D Bv 2 . (45)
The bilinear terms in the left hand side of this equation form the finite element matrix M that is sent to ARPACK. To check the accuracy of the eigenvalue computation, we compute the residual:

M A B -ω A B ∞ . ( 46
)
Numerical tests showed that using a shift leads to an increased accuracy: the residual decreases to 10 -7 and eigenvalues are closer to the expected values when compared to known results. We use a shift σ = 10 -4 , which is implemented by adding the following term to the matrix:

- D σ(Av 1 + Bv 2 ). ( 47
)
For the two-component case, the numerical method is similar and based on the following weak formulation corresponding to (20)-( 22):

                                                               D 1 2 ∇A • ∇v 1 + D (C trap -µ 1 )Av 1 + D 2β 11 |φ 1 | 2 + β 12 |φ 2 | 2 Av 1 + D β 11 φ 2 1 Bv 1 + D β 12 φ 1 φ 2 Cv 1 + D β 12 φ 1 φ 2 Dv 1 = ω D Av 1 , - D 1 2 ∇B • ∇v 2 - D (C trap -µ)Bv 2 - D 2β 11 |φ 1 | 2 + β 12 |φ 2 | 2 Bv 2 - D β 11 φ 1 2 Av 2 - D β 12 φ 1 φ 2 Cv 2 - D β 12 φ 1 φ 2 Dv 2 = ω D Bv 2 , D 1 2 ∇C • ∇v 3 + D (C trap -µ)Cv 3 + D 2β 22 |φ 2 | 2 + β 21 |φ 1 | 2 Cv 3 + D β 21 φ 1 φ 2 Av 3 + D β 21 φ 1 φ 2 Bv 3 + D β 22 φ 2 2 Dv 3 = ω D Cv 3 , - D 1 2 ∇D • ∇v 4 - D (C trap -µ)Dv 4 - D 2β 22 |φ 2 | 2 + β 21 |φ 1 | 2 Dv 4 - D β 21 φ 1 φ 2 Av 4 - D β 21 φ 1 φ 2 Bv 4 - D β 22 φ 2 2 Cv 4 = ω D Dv 4 . (48)

Validation test cases for the one-component BEC

We start by validating the codes for the one-component BEC against well known benchmarks. A summary of the considered cases, together with typical computational times and mesh sizes, is provided in Tab. 1. When mesh adaptation is used, we indicate the size of the mesh for the last step of the continuation procedure. Note that we considered for all cases the non-dimensional equations in the setting for which β = 1.

1D case: ground state

The first test case is the computation of eigenvalues of the ground state of a one-dimensional BEC with trapping potential V trap = 1 2 mω 2 z z 2 . In the Thomas-Fermi limit, the explicit expressions for eigenvalues are known (Kevrekidis

ω TF n = ω z n(n + 1) 2 , n ∈ N. ( 49
)
We compare in Tab. 2 numerical and theoretical values of eigenvalues ω for µ = 6 and ω z = 0.025. The Thomas-Fermi solution (9) was used to initialize the Newton algorithm. We could check from Tab. 2 that the computed eigenvalues verify the following expected properties (see also Sect.

2.2):

(i) all eigenvalues are real (the stationary state is dynamically stable) and form pairs (+ω, -ω), (ii) the first eigenvalue is ω = 0, (iii) the other eigenvalues correspond to theoretical predictions (49), (iv) all Krein signatures are positive (the stationary state is energetically stable).

1D case: dark soliton

We analyse for the second 1D test case an excited state, obtained by adding a dark solition to the Thomas-Fermi density previously computed. The initial condition for the Newton algorithm is thus built as:

φ init DS = √ n TF tanh( √ µz). ( 50
)
We plot in Fig. 1(a) the initial condition and the converged stationary state. Eigenvalues are displayed in Tab. 3. As expected, all eigenvalues are real, as the dark soliton is dynamically stable in 1D. A complete characterization of the BdG modes is offered in [START_REF] Frantzeskakis | Dark solitons in atomic Bose-Einstein condensates: from theory to experiments[END_REF].

Re(ω)

Im(ω) K ω TF n from ( 49) ω 1 -2.89857e-15 2.16087e-07 1 ω TF 0 = 0 ω 2 6.18933e-15 -2.16087e-07 1 ω 3 -0.025 -8.80682e-11 1 ω TF 1 = ω z = 0.025 • The mode with ω 4 ≈ ωz √ 2 ≈ 0.017677669 is the anomalous mode; it is the only mode with a negative Krein signature. It is represented in Fig. 1(b) and we retrieve the profile obtained in [START_REF] Law | Quantum fluctuations of coupled dark solitons in a trapped Bose-Einstein condensate[END_REF].

ω
• The dipole or Kohn mode at ω 6 ≈ ω z = 0.025 corresponds to oscillations of the center of mass of the condensate.

• The quadrupole mode (or the breathing mode) is obtained for ω 8 ≈ ω z √ 3 ≈ 0.04330127. This mode is particular to the one-dimensionality of the system.

Re(ω)

Im(ω) K ω 1 -2.57971e-07 -3.29591e-15 1 ω 2 2.57971e-07 3.25101e-15 1 ω 3 -0.0178197 -1.09577e-12 -1 ω 4 0.0178197 1.00791e-12 -1 ω 5 -0.025 -5.42312e-12 1 ω 6 0.025 6.69980e-12 1 ω 7 -0.0435553 -8.24994e-12 1 ω 8 0.0435553 9.60204e-12 1 ω 9 -0.0616151 4.28088e-13 1 ω 10 0.0616151 2.16934e-13 1

2D case: ground state

We switch now to 2D BEC configurations with trapping potential V trap = 1 2 mω 2 ⊥ r 2 , where r 2 = x 2 + y 2 . For this case, the eigenvalues in the Thomas-Fermi limit are also known [START_REF] Kevrekidis | Distribution of eigenfrequencies for oscillations of the ground state in the Thomas-Fermi limit[END_REF]:

ω TF m,k = ω ⊥ m + 2k 2 + 2k(1 + m), ( 51
)
where m, k ≥ 0 are integers. We present in Tab. 4 the first 20 eigenvalues computed for µ = 6 and ω ⊥ = 0.2, with and without mesh adaptation. We find all Krein signatures to be 1, which is the expected result in the absence of topological excitations. This is a perfect case to check that computations using mesh adaptation provide the same results as computations with a refined fixed mesh. Results in Tab. 4 show that this is indeed the case for our BdG solvers.

2D case: dark soliton

Following the same procedure as in the 1D case, we add to the previously computed 2D ground state a centered dark soliton (Fig. 2). This is an

No mesh adaptation

With mesh adaptation Re(ω) Im(ω) K Re(ω) Im(ω) K ω m,k from ( 51) ω 1 -2.07687e-06 5.60174e-16 1 -6.24135e-15 1.37474e-07 1 ω TF 0,0 = 0 ω 2 2.07687e-06 -5.08061e-16 1 6.23144e-15 -1.37474e-07 1 ω 3 -0.2 -4.04752e-11 1 -0.2 1.16809e-11 1 interesting case to test the continuation procedure in following a branch of stationary solutions. The initial condition is given by the |1, 0 state in the linear limit:

ω TF 1,0 = 0.2 ω 4 0.2 1.28499e-
φ DS = ω ⊥ 2π H 0 ( √ ω ⊥ x)H 1 ( √ ω ⊥ y)e -1 2 ω ⊥ (x 2 +y 2 ) , ( 52
)
where H n are Hermite polynomials. We set ω ⊥ = 0.2 and follow this solution up to µ = 3. Real and imaginary parts of eigenvalues are presented in Fig. 2(a,b) and are identical to the results published in Middelkamp et al. (2010a). This state does not have an azimuthal symmetry. Due to the space discretization, there exists a preferred direction along which the soliton will tend to align itself. When adapting the mesh, this direction changes and the wave function will then rotate. To avoid this phenomenon, we only adapt the mesh every 5 iterations during the continuation procedure. This permits to optimally adapt the size of the mesh while reducing the effects of the rotation.

Table 1 shows that this is an efficient approach to reduce the computational time. The final adapted mesh is presented in Fig. 2(c), with the corresponding atomic density |φ| 2 in Fig. 2(d).

2D case: central vortex

We compute another solution studied in Middelkamp et al. (2010a,b). It consists of a disk-shaped BEC with a centered vortex. The initial condition is given by the |0, 1 state in cylindrical coordinates (r, θ):

φ V S ∝ rL 1 0 (ω ⊥ r 2 )e iθ e -1 2 ω ⊥ r 2 , ( 53
)
where L 1 0 is the Laguerre polynomial. We set, as in the previous case, ω ⊥ = 0.2. Eigenvalues computed with and without mesh adaptation are displayed in Fig. 3(a) and (b). We checked that both methods give the same results as those obtained in Middelkamp et al. (2010a,b). The atomic density is presented in Fig. 3(c) and (d) for two values of µ. The first BdG modes (A, B) for this state are displayed in Fig. 4 by plotting their modulus coloured by the phase. We can distinguish:

• The zero-energy mode ω = 0 (Fig. 4(a)), associated to the phase invariance of the GP equation.

• The anomalous mode (Fig. 4(b)) corresponds to the following approximation of its eigenvalue in the Thomas-Fermi limit (Middelkamp et al., 2010a):

ω = ω 2 ⊥ 2µ ln(A µ ω ⊥ ) ≈ 0.03261667238, A ≈ 2 √ 2π ≈ 8.886. ( 54
)
This value is represented by a dashed green line in Figs. 3 (a) and (b).

• The dipole or Kohn mode (Fig. 4(c)) corresponds to ω = ω ⊥ and is independent of µ.

• The 4-th mode (Fig. 4(d)) corresponding to ω = µ -2ω ⊥ in the linear limit (Middelkamp et al., 2010b). This value is represented by a continuous green line in Figs.

3D case: ground state

Computing the BdG modes for a 3D BEC is a challenging numerical problem.

Even with mesh adaptation, the number of degrees of freedom is high and increases with the size of condensate (i. e. with µ). For this test case, we computed the spectrum of the ground state of a spherical BEC with trapping potential V trap = 1 2 mω 2 ⊥ r 2 , where r 2 = x 2 + y 2 + z 2 . We set ω ⊥ = 1. The eigenvalues presented in Fig. 5 are in very good agreement with numerical results obtained in [START_REF] Bisset | Robust vortex lines, vortex rings, and hopfions in three-dimensional Bose-Einstein condensates[END_REF]. This case shows that our finite-element toolbox can be used to study simple 3D configurations. For more complicated states, the use of parallelization is mandatory to reduce the computational time and memory requirements.

Validation test cases for the two-component BEC

For the two-component BEC, we compute BdG modes for the dark-antidark solitary 1D or 2D waves studied in [START_REF] Danaila | Vector dark-antidark solitary waves in multicomponent Bose-Einstein condensates[END_REF]. Antidark solitary waves are bright solitary waves on top of a finite background. Such states appear in a two-component system with inter-component repulsion: a dark soliton or a vortex (or a ring) in one component will induce an effective potential, through the inter-component nonlinearity, on the second component.

The result is that atoms of the second component are attracted into the dip of the first one. We consider the system (31) in the case of repulsive intercomponent interactions with miscibility condition 0 ≤ β 12 < √ β 11 β 22 that ensures that the two components co-exist outside the dark-antidark state. To simplify the case study, since only the ratio between non-linear interaction constants matters, we set

β 11 = β 22 = β = 1, β 12 = β 21 and 0 < β 12 < β.
The considered two-component cases are summarized in Tab. 5 displaying the necessary computational times and mesh sizes. Note that the codes consider independent values for coefficients β ij , 1 ≤ i, j ≤ 2, and thus can be used to study configurations different from those analysed in [START_REF] Danaila | Vector dark-antidark solitary waves in multicomponent Bose-Einstein condensates[END_REF].

1D case: dark-antidark soliton

The first state is a dark-antidark solitary wave in 1D. We set a soliton solution (constructed as in Eq. ( 50)) in the first component, while the second component is in the Thomas-Fermi ground state. Obtained eigenvalues are shown in Fig. 6(a) and (b) and correspond to the results of [START_REF] Danaila | Vector dark-antidark solitary waves in multicomponent Bose-Einstein condensates[END_REF]. The small imaginary instability around β 12 = β 21 = 0.8 is well resolved.

Profiles of the atomic density for different values of the interaction coefficient are presented in Fig. 6(c).

2D two-component case: ring-antidark-ring state

With the toolbox, we provide two test cases for 2D configurations: the vortexantidark and ring-antidark solitary waves. We show here only the case with ring-antidark solitary waves. The first component contains a ring soliton and the second is in the ground state. Results are shown in Fig. 7(a) and (b) for the real and imaginary parts of the eigenvalues. Figure 7(c) shows for the atomic density profiles which correspond to the figures presented in [START_REF] Danaila | Vector dark-antidark solitary waves in multicomponent Bose-Einstein condensates[END_REF].

Description of the programs

In this section, we first describe the architecture of the programs and the organisation of the provided files. We then present the input parameters and the structure of the output files.

Program architecture

Codes and data files forming the BdG toolbox are stored in the FFEM BdG toolbox directory, which is organized around two main subdirectories: BdG 1comp and BdG 2comp, corresponding to the one-and two-component codes. Each subdirectory contains two main files: FFEM GP $case.edp, which is the main FreeFem++ script file for the computation of the stationary state, and FFEM BdG $case.edp which is the main FreeFem++ script file for the computation of the BdG eigenvalues. To run the computation of the Gross-Pitaevskii stationary state, the user can use the command FreeFem++ FFEM GP $case.edp. BdG eigenvalues can then be computed with the command FreeFem++ FFEM BdG $case.edp. Parameter files for the examples presented in this paper are stored in the INIT folder.

The obtained solutions are saved in the dircase directory. Depending on the output format selected by the user, data files are generated in specific folders for visualization with Tecplot, Paraview or Gnuplot. We also provide in the folder Figures ready-made layouts for Tecplot. The user can thus obtain the figures from this paper using newly generated data. More details about the output structure are given in Sect. 7.4.

The complete architecture of the BdG 1comp directory is the following (the architecture of the BdG 2comp directory is almost identical):

1. FFEM GP $case.edp: the main script for the computation of the GP stationary states.

2. FFEM BdG $case.edp: the main script for computing eigenvalues.

3. param num common.inc: a parameter file for the main numerical parameters. For some examples, it is necessary to run the case with different parameters (e.g. with and without mesh adaptation) before opening the layout.

6. A macro: directory containing macros used in the main scripts.

Macros and functions

The different macros and functions used in the toolbox for the sequential code are stored in the A macro folders:

• Macro BdGsolve.edp: macro computing the BdG eigenvalues corresponding to matrices in Eqs. ( 45) and ( 48).

• Macro createdir.edp: macro creating the file structure of the dircase folder.

• Macro GPsolve.edp: macro computing the GP stationary state with a Newton method (see Eqs. ( 36) and ( 40)-( 43)).

• Macro meshAdapt.edp: macro adapting the mesh to the wave function.

• Macro operator.edp: definitions of useful macros and functions: gradients, energy (3), chemical potential (7), Hermite polynomials, etc. Also contains a macro creating a spherical mesh for 3D problems.

• Macro output.edp: macros used to save data in Tecplot and Paraview formats.

• Macro plotEigenvector.edp: macro plotting the real and imaginary parts of a BdG eigenvector.

• Macro plotphi.edp: macro plotting the complex wave function. The user can press "k" to alternate between plots of the density, phase and real and imaginary parts of the wave function.

• Macro problem.edp: definitions of the weak formulations for the GP (Eqs. ( 36) or ( 48)) and BdG problems (Eqs. ( 45) or ( 40)-( 43)).

• Macro restart.edp: macros used to save and load the wave function to or from FreeFem++ files.

• Macro saveData.edp: macro saving the stationary wave function.

• Macro saveEigenvalues.edp: macro saving the BdG eigenvalues and eigenvectors.

Input parameters

Parameters are separated in two files. Numerical parameters used in all computations are in param num common.inc. Files in the INIT directory specify physical parameters describing the state that will be studied during a computation and numerical parameters specific to this problem. The files distributed with the toolbox provide a variety of examples that can be used as a starting point when selecting parameters for the study of new states.

(1) In the file param num common.inc, the parameters are:

• displayplot: controls the output information to plot. Possible values range from 0 (no plots), to 2 (plots data at all iterations of the Newton method and all eigenvectors computed by the BdG code).

• iwait: a Boolean indicating if the code must wait for user input when a plot is shown (true) or it can continue (false) with the next plot.

• cutXY, cutXZ, cutYZ: (only for 3D cases) Booleans indicating whether to plot cuts of the wave function along the different axis at x = 0, y = 0 or z = 0.

• Tecplot: a Boolean indicating whether to save data in the Tecplot format.

• Paraview: a Boolean indicating whether to save data in the Paraview format (only in 2D and 3D).

• adaptinit: if true, the initial solution is recomputed after the first mesh adaptation.

• adaptmeshFF: determines if mesh adaptation is used (true) or not (false).

• useShift: a Boolean indicating whether to use a shift when computing the BdG eigenvalues (see Eq. ( 47)).

• Nadapt: if mesh adaptation is used, the mesh is adaptated every Nadapt iterations during the continuation.

• Nplot: the wave function is plotted every Nplot iterations during the continuation.

• Nsave: the wave function is saved for Paraview or Tecplot every Nsave iterations during the continuation.

• Nrst: the wave function is saved for the BdG computation every Nrst iterations during the continuation.

• tolerrF: the tolerance value of F in Eq. ( 44).

• tolNewton: the tolerance value of q in Eq. ( 44).

• shift: the value of the shift σ used when computing eigenvalues.

• newtonMax: the maximum number of Newton iterations.

(2) In the file $case.inc, stored in the INIT directory, the parameters are:

• General parameters for the case:

• dimension: the dimension of the problem (1, 2 or 3).

• FEchoice: the type of finite element used. Usually P 2.

• nev: the number of eigenvalues computed by the BdG code.

• Parameters used to restart a computation:

• restart: a boolean indicating if the initial solution is a restart from a previous computation. If true, the solution and mesh stored in fcaserestart for the value of µ given by murestart will be used as initial solution.

• murestart: the initial value of µ in the case of a restart.

• fcaserestart: the folder where the initial solution is stored in the case of a restart.

• Parameters of the continuation:

• kpol, lpol, mpol: integers defining the initial state in the linear limit.

• startmu: the initial value of µ.

• endmu: the final value of µ.

• dmu: the increment in µ during the continuation.

• facmu: when using the linear limit, the initial value of µ is given by facmu • µ |klm .

• Coefficients of the GP equation:

• beta: the nonlinear coefficient (we set β = 1 in all test cases).

• ax, ay, az: the frequencies of the trapping potential along the three axes.

• Ctrap: a function defining the trapping potential.

• Parameters for the mesh creation:

• Dx: the distance between points on the mesh border.

• scaledom: a coefficient used to control the size of the domain: the mesh radius is given by Rdom = scaledom r TF where r TF is the Thomas-Fermi radius.

• createMesh: a macro creating the initial mesh Th.

• Parameters for the mesh adaptation:

• errU: the interpolation error level.

• hmin: the minimum length of a mesh element edge in the new mesh.

• hmax: the maximum length of a mesh element edge in the new mesh.

• adaptratio: the ratio for a prescribed smoothing of the metric. No smoothing is done if the value is less than 1.1.

• Parameters for the initial solution:

• initname: the name given to the initial solution.

• initcond: a macro defining the initial solution for the phi variable.

• Definitions of the boundary conditions:

• BCGP: the boundary conditions used in the GP code for Eqs. ( 36) and ( 40)-( 43).

• BCBdG: the boundary conditions used in the BdG code for Eqs. ( 45) and (48).

• fcase: the name given to the current computation.

• dircase: the directory where the results are stored.

1. The Gnuplot folder contains two files:

• Informations about the stationary states are stored in the GP results.dat file. The columns are in order: the non-linear coefficient β (β 12 and β 21 ), the imposed chemical potential µ (µ 1 and µ 2 ), the number of Newton iterations used for this value of µ, the errors F and q (44), the computed value of the chemical potential ( 7) (computed values of µ 1 and µ 2 ), the number of atoms (4) (the number of atoms in the two components), the GP energy (3), the mesh size, the number of degrees of freedom and the CPU time to compute the stationary state.

• BdG eigenvalues are stored in the BdG results.dat file. The columns are in order: the non-linear coefficient β (β 12 and β 21 ), the imposed chemical potential µ (µ 1 and µ 2 ), the eigenvalue number between 0 and nev, the real and imaginary part of the eigenvalues, the Krein signature and its sign (the Krein signature and its sign for the two components), the residual ( 46) and the CPU time to compute the eigenvalues.

2. The Paraview folder contains the wave functions stored as .vtk or .vtu and .pvd files:

• phi init.vtu and phi final.vtu are the initial and final solutions.

• phi mu $mu.vtu contains the stationary wave function for a given value of µ.

• phi mu1 $mu1 mu2 $mu2.vtu contains the stationary wave function for given values of µ 1 and µ 2 in the first continuation.

• phi beta12 $beta12 beta21 $beta21.vtu contains the stationary wave function for given values of β 12 and β 21 in the second continuation.

3. The Paraview Eigenvectors folder contains the eigenvectors stored as:

• eVec mu $mu $nev.vtu in the one component code.

• eVec beta12 $beta12 beta21 $beta21 mu1 $mu1 mu2 $nev.vtu in the two component code.

4. The RST folder contains the stationary states stored as FreeFem++ files.

The names are:

• RST-$mu.rst or RST-$mu1-$mu2-$beta12-$beta21.rst for the data.

• RSTTh-$mu or RSTTh-$mu1-$mu2-$beta12-$beta21 for the mesh files. The file extensions are .mesh (in 1D), .msh (in 2D) or .meshb (in 3D).

5. The Tecplot folder contains the wave functions stored as .dat Tecplot files:

• phi init.dat and phi final.dat are the initial and final solutions.

• phi mu $mu.dat contains the stationary wave function for a given value of µ.

• phi mu1 $mu1 mu2 $mu2.dat contains the stationary wave function for given values of µ 1 and µ 2 in the first continuation.

• phi beta12 $beta12 beta21 $beta21.dat contains the stationary wave function for given values of β 12 and β 21 in the second continuation.

6. The Tecplot Eigenvectors folder contains the eigenvectors stored in the Tecplot format:

• eVec mu $mu $nev.dat in the one component code.

• eVec beta12 $beta12 beta21 $beta21 mu1 $mu1 mu2 $nev.dat in the two component code.

7. The Tecplot Eigenvalues folder contains the eigenvalues stored in the Tecplot format. Filenames are BdG results $i.dat. Each file contains the i-th smallest eigenvalue for each value of µ (or β 12 in the two-component code).

Summary and conclusions

The aim of the toolbox presented in this paper is the computation of stationary states and BdG modes of one-and two-component BECs in 1D and 2D. The use of mesh adaptation enables an efficient computation of stationary states by adapting the spatial discretization to the topology of the considered state. This makes possible the study of various 1D and 2D problems and even a simple 3D configuration without parallelization. The toolbox was created with FreeFem++ , a free and open-source finite element software for the study of partial differential equations. The method consists of two steps: (i) a Newton method, combined with a continuation on the chemical potential µ or the inter-component interaction parameters β 12 and β 21 , is used to obtain the stationary state of the GP equation, (ii) the BdG modes are obtained by solving an eigenvalue problem with ARPACK. The numerical code was validated against test cases studied theoretically or numerically in the literature. All parameter files corresponding to these test cases are shared with the toolbox, allowing the user to obtain the results presented in this paper. These parameter files can be used as templates to initiate the study of other BEC states. We considered only one and two component systems with a cubic nonlinearity, but the toolbox could be easily modified to study other configurations such as quartic ± quadratic trapping potentials [START_REF] Bretin | Fast rotation of a Bose-Einstein condensate[END_REF], dipolar interactions [START_REF] Tang | A spectrally accurate numerical method for computing the Bogoliubov-de Gennes excitations of dipolar Bose-Einstein condensates[END_REF] or spinor condensates [START_REF] Mithun | Existence, stability, and dynamics of monopole and Alice ring solutions in antiferromagnetic spinor condensates[END_REF]. Future extensions of the toolbox concern the implementation of domain decomposition methods and the use of PETSc and SLEPc libraries for this problem. The new parallel toolbox will make possible the computation (with a reasonable CPU time), of BdG modes for 3D BEC configurations without axial symmetry, such as BEC with U , S vortices [START_REF] Aftalion | Three-dimensional vortex configurations in a rotating Bose-Einstein condensate[END_REF], giant vortices [START_REF] Danaila | Three-dimensional vortex structure of a fast rotating Bose-Einstein condensate with harmonic-plus-quartic confinement[END_REF] or new computationally discovered exotic states reported in very recent contributions [START_REF] Boullé | Two-component threedimensional atomic Bose-Einstein condensates supporting complex stable patterns[END_REF].

Figure 1 :

 1 Figure 1: 1D dark soliton: a) initial state and stationary solution, b) anomalous mode with ω 4 ≈ ωz √ 2 .

Figure 2 :

 2 Figure 2: 2D dark soliton: a) real part ω r and b) imaginary part ω i of eigenvalues as a function of µ. Solution for µ = 3.007: c) adapted mesh and d) atomic density |φ| 2 .

Figure 3 :

 3 Figure 3: 2D BEC with a centered vortex. Real part ω r of eigenvalues as a function of µ computed a) without and b) with mesh adaptation. c) Atomic density |φ| 2 for µ = 0.45 and µ = 3.

  3 (a) and (b).

Figure 4 :

 4 Figure 4: 2D BEC with a centered vortex: first four BdG modes A and B. Surfaces of modulus coloured by the phase.

Figure 5 :

 5 Figure 5: 3D case: ground state. a) Real part ω r of the eigenvalues as a function of µ, b) illustration of four BdG modes (iso-surfaces of the modulus).

Figure 6 :

 6 Figure 6: 1D two-component case: dark-antidark solitary wave. a) Real part ω r and b) imaginary part ω i of BdG eigenvalues, c) atomic density profiles for three values of β 12 .

Figure 7 :

 7 Figure 7: 2D two-component case: ring-antidark-ring state. a) Real part ω r and b) imaginary part of eigenvalues, c) density profiles for three values of β 12 .

4.

  INIT: directory storing the parameter files for the examples presented in Sect. 5. 5. Figures: directory containing Tecplot layouts used to replot the figures shown in Sect. 5. The main code must be run with the associated example before opening the layout to replot the figure.

Table 2 :

 2 1D ground state: eigenvalues and Krein signatures.

	4	0.025 2.76512e-11 1
	ω 5 ω 6	-0.0433018 -4.41549e-11 1 ω TF 2 ≈ 0.043301270 0.0433018 -1.21387e-11 1
	ω 7 ω 8	-0.0612394 -2.87955e-10 1 ω TF 3 ≈ 0.061237243 0.0612394 1.64467e-10 1
	ω 9 ω 10	-0.0790624 -1.09235e-10 1 ω TF 4 ≈ 0.07905694 0.0790624 8.67993e-11 1

Table 3 :

 3 1D dark soliton: eigenvalues and Krein signatures.

Table 4 :

 4 2D ground state: eigenvalues and Krein signatures.

Table 5 :

 5 Test cases for the two-component BEC. Computational time and mesh size (number of elements). All computation were performed on a Macbook pro M1, 16GB of DDR4 2400 MHz RAM.

		Without mesh adaptation		With mesh adaptation
		CPU time GP CPU time BdG mesh size CPU time GP CPU time BdG mesh size
	1D dark-antidark state	00:00:30	00:08:48	2714
	2D vortex-antidark state	00:07:46	01:02:21	10469	00:14:34	01:27:09	7874
	2D ring-antidark state	00:09:22	01:48:45	10469	00:11:05	01:59:30	9533

Acknowledgements

The authors acknowledge financial support from the French ANR grant ANR-18-CE46-0013 QUTE-HPC. Part of this work used computational resources provided by IDRIS (Institut du développement et des ressources en informatique scientifique) and CRIANN (Centre Régional Informatique et d'Applications Numériques de Normandie). The authors are grateful to P. K. Kevrekidis for stimulating discussions and warmly thank E. Charalampidis for his helpful input in validating numerical codes.

(3) In a two component case, some new parameters are defined in the $case.inc file:

• Parameters used to restart a computation:

• mu1restart, mu2restart: initial values of µ 1 and µ 2 in the case of a restart.

• beta12restart, beta21restart initial values of β 12 and β 21 in the case of a restart.

• Parameters of the continuation:

• startmu1, startmu2: initial values of µ 1 and µ 2 .

• endmu1, endmu2: final values of µ 1 and µ 2 .

• dmu1, dmu2: increments of µ 1 and µ 2 during the continuation.

• startbeta12, startbeta21: initial values of β 12 and β 21 .

• endbeta12, endbeta21: final values of β 12 and β 21 .

• dbeta12, dbeta21: increments of β 12 and β 21 during the continuation.

• Coefficients of the GP equation:

• beta11, beta12: nonlinear coefficients β 11 and β 22 .

• Parameters for the initial solution:

• initname1: the name given to the initial solution for the first component.

• initname2: the name given to the initial solution for the second component.

• initcond: a macro defining the initial solution for [phi1,phi2] variables.

Outputs

When a computation starts, the OUTPUT $case directory is created. It contains up to eight folders. The RUNPARAM directory contains a copy of the code and data files, allowing an easy identification of each case and preparing an eventual rerun of the same case. The other folders contains different output format files of the computed solution, to be visualised with Tecplot, Paraview or Gnuplot. The content of these subfolders depends on the case and on the computation parameters (differences in the two component code are given in parenthesis):


à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

The rhizosphere of Sulla spinosissima growing in abandoned mining soils is a reservoir of heavy metals tolerant plant growth-promoting rhizobacteria Heavy metals tolerant soil bacteria are known to play essential roles in biogeochemical cycles, biotransformation of metals, bioremediation processes and, plant adaptation. The objective of this study was to isolate and characterize the bacterial strains associated with Sulla spinosissima, a native legume species growing in three abandoned mining sites situated in the Oujda region (East Morocco). Globally, more than 370 bulk and rhizospheric soil bacteria were isolated. Their identification by 16S rDNA sequencing showed that dominant phyla were Firmicutes, Actinobacteria, and Proteobacteria, while at the genus level Bacillus dominated, followed by Stenotrophomonas, Arthrobacter, and Rhodococcus. All the soils contained strains possessing plant growth-promoting traits. The best-performing strains were LMR283 for auxin production (145 μg mL -1 ), LMR291 for phosphate solubilization (67.6 mg L -1 ), LMR280 for siderophore pro- duction (92.4%) and LMR326 for ACC deaminase activity (105 nmoL α-ketobutyrate mg -1 h -1 ).

Among tolerant PGP bacteria, 17 isolates showed antagonistic activity against the pathogen Fusarium oxysporum and 26 produced lytic enzymes. It was relevant that the rhizospheric soils prospected compared to bulk soils contained more interesting isolates for all the studied properties, in particular for soils sampled from Oued El Heimer and Sidi Boubker sites (respectively 79 and 63% of performing strains). Data presented here indicate that Sulla spinosissima growing in heavy metal soils is associated with multifarious active bacterial populations that probably sustain plant tolerance and growth under the prevailing stressful conditions. Superior strains identified are good candidates to be used with selected plants in rehabilitation programs in the contaminated ecosystems.

Introduction

Environmental pollution is continuously increasing because of the acceleration of human activities and economic development. In particular, industrial production increases the number of foreign and toxic substances released in the environment. Many of these