Xavier Antoine
email: xavier.antoine@univ-lorraine.fr

Romain Duboscq
email: romain.duboscq@univ-lorraine.fr

Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates

 for computing the stationary states of Bose-Einstein Condensates (BECs) through the Gross-Pitaevskii equation. We show that the fixed point approach introduced in [10] fails to converge for fast rotating BECs. A simple alternative approach based on Krylov subspace solvers with a Laplace or Thomas-Fermi preconditioner is given. Numerical simulations (obtained with the associated freely available Matlab toolbox GPELab) for complex configurations show that the method is accurate, fast and robust for 2D/3D problems and multi-components BECs.

Introduction

Bose-Einstein Condensates (BECs) were first experimentally observed in 1995 [START_REF] Anderson | Observation of Bose-Einstein condensation in a dilute atomic vapor[END_REF][START_REF] Bradley | Evidence of Bose-Einstein condensation in an atomic gas with attractive interactions[END_REF][START_REF] Dalfovo | Theory of Bose-Einstein condensation in trapped gases[END_REF][START_REF] David | Bose-Einstein Condensation in gas of sodium atoms[END_REF] while they were theoretically predicted a long time before by S.N. Bose and A. Einstein. This state of matter allows to study quantum physics at the macroscopic scale. Later on, quantum vortices have been observed [START_REF] Abo-Shaeer | Observation of vortex lattices in Bose-Einstein condensates[END_REF][START_REF] Bretin | Fast rotation of a Bose-Einstein condensate[END_REF][START_REF] Madison | Stationary states of a rotating Bose-Einstein condensate: Routes to vortex nucleation[END_REF][START_REF] Madison | Vortex formation in a stirred Bose-Einstein condensate[END_REF][START_REF] Matthews | Vortices in a Bose-Einstein condensate[END_REF][START_REF] Raman | Vortex nucleation in a stirred Bose-Einstein condensate[END_REF][START_REF] Yuce | Off-axis vortex in a rotating dipolar Bose-Einstein condensate[END_REF] and a growing interest has been directed towards the understanding of rotating BECs. At temperatures T much smaller than the critical temperature T c , the macroscopic behavior of a rotating BEC can be well described by a condensate wave function ψ(t, x) which is solution to the Gross-Pitaevskii Equation (GPE) with rotating term. More precisely, ψ(t, x) is solution to the dimensionless time-dependent GPE [START_REF]Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF]

i∂ t ψ(t, x) = - 1 2 ∆ψ(t, x) + V (x)ψ(t, x) + βf (|ψ(t, x)|)ψ(t, x) -Ω • Lψ(t, x), (1)
for x ∈ R d , d = 2, 3, t > 0. In 3D, the Laplace operator is defined by: ∆ = ∇ 2 , where ∇ := (∂ x , ∂ y , ∂ z ) t is the gradient operator; the spatial variable is x = (x, y, z) t ∈ R 3 (in 2D we have ∇ := (∂ x , ∂ y ) t and x = (x, y) t ∈ R 2 ). Function V is the external (usually confining) potential.

Parameter β is the nonlinearity strength describing the interaction between atoms of the condensate. This parameter is related to the s-scattering length (a s ) and is positive (respectively negative) for a repulsive (respectively attractive) interaction. Function f describes the nonlinearity arising in the problem, which is fixed to the cubic case in the paper: f (|ψ|) = |ψ| 2 (but other cases could be considered like e.g. cubic-quintic problems or integral nonlinearities for dipolar gazes [START_REF]Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF]). For vortices creation, a rotating term is added. The vector Ω is the angular velocity vector and the angular momentum is L = (L x , L y , L z ) = x ∧ P, with the momentum operator P = -i∇. In many situations and all along the paper, the angular velocity is such that Ω = (0, 0, Ω) t leading to

Ω • L = ΩL z = -i(x∂ y -y∂ x ). (2)
In addition, an initial data ψ(t = 0, x) = ψ 0 (x) is prescribed to get a complete system. One important problem in the numerical solution of the GPE is the computation of stationary states which consists in finding a solution ψ(t, x) = e -iµt φ(x),

where µ is the chemical potential of the condensate and φ is a time independent function. This function is given as the solution to the nonlinear elliptic equation

µφ(x) = - 1 2 ∆φ(x) + V (x)φ(x) + β|φ(x)| 2 φ(x) -ΩL z φ(x), (4)
under the normalization constraint

||φ|| 2 0 = R d |φ(x)| 2 dx = 1, (5)
where || • || 0 is the L 2 -norm in R d . This nonlinear eigenvalue problem can be solved by computing the chemical potential

µ β,Ω (φ) = E β,Ω (φ) + β 2 R d |φ(x)| 4 dx, (6)
with

E β,Ω (φ) = R d 1 2 |∇φ| 2 + V |φ| 2 + β 2 |φ| 4 -Ωℜ (φ * L z φ) dx, (7)
where φ * is the complex conjugate of φ. This also means that the eigenvalues are the critical points of the energy functional E β,Ω over the unit sphere: S := {||φ|| 0 = 1}. Furthermore, (4) can be seen as the Euler-Lagrange equation associated with the constraint minimization problem [START_REF] Bao | Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow[END_REF].

Various numerical methods [START_REF]Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF] have been designed to compute stationary states of GPEs, most particularly to get the ground state and the excited state solutions. It is admitted that deriving robust and efficient numerical approaches for the stationary state computation is difficult, most particularly when the nonlinearity is large and the rotation velocity is high. More generally, methods are based either on solving the nonlinear eigenvalue problem [START_REF] Dion | Ground state of the time-independent Gross-Pitaevskii equation[END_REF][START_REF] Wang | A spectral-Galerkin continuation method using Chebyshev polynomials for the numerical solutions of the Gross-Pitaevskii equation[END_REF][START_REF] Wang | A two-parameter continuation method for rotating two-component Bose-Einstein condensates in optical lattices[END_REF] or deriving nonlinear optimization techniques under constraints [START_REF] Bao | Ground-state solution of Bose-Einstein condensate by directly minimizing the energy functional[END_REF][START_REF] Caliari | A minimisation approach for computing the ground state of Gross-Pitaevskii systems[END_REF][START_REF] Caliari | GSGPEs: A MATLAB code for computing the ground state of systems of Gross-Pitaevskii equations[END_REF][START_REF] Danaila | A finite element method with mesh adaptivity for computing vortex states in fast-rotating Bose-Einstein condensates[END_REF][START_REF] Danaila | A new Sobolev gradient method for direct minimization fo the Gross-Pitaevskii energy with rotation[END_REF]. An alternative approach is the imaginary time method. It has been used extensively by the Physics community and has proved to be powerful [START_REF] Adhikari | Numerical solution of the two-dimensional Gross-Pitaevskii equation for trapped interacting atoms[END_REF][START_REF] Baye | Resolution of the Gross-Pitaevskii equation with the imaginarytime method on a Lagrange mesh[END_REF][START_REF] Cerimele | Numerical solution of the Gross-Pitaevskii equation using an explicit finite-difference scheme: An application to trapped Bose-Einstein condensates[END_REF][START_REF] Chiofalo | Ground state of trapped interacting Bose-Einstein condensates by an explicit imaginary-time algorithm[END_REF][START_REF] Edwards | Numerical solution of the nonlinear Schrödinger equation for small samples of trapped neutral atoms[END_REF][START_REF] Gammal | Improved numerical approach for the time-independent Gross-Pitaevskii nonlinear Schrödinger equation[END_REF]. From the mathematical point of view, the imaginary time method has been studied in [START_REF]Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF][START_REF] Bao | Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow[END_REF] and written as a gradient flow formulation for the non rotating case. In particular, the authors show that the time discretization of the Gradient Flow formulation must be carefully considered. It is shown that the (semi-implicit) Backward Euler scheme is particularly well adapted since it leads to an energy diminishing formulation without any CFL constraint on the time step. For the rotating case, such a result is not proved. Concerning the spatial discretization, finite difference (or finite element) methods can be used. However, for correctly capturing the nucleation of vortices, high-order discretization schemes or refined (adaptive) meshes must be used. These difficulties clearly complicate the construction of a simple and versatile numerical method, most particularly for two-and three-dimensional problems [START_REF]Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF], when considering dipolar-dipolar interactions [START_REF]Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF][START_REF] Bao | Efficient numerical methods for computing ground states and dynamics of dipolar Bose-Einstein condensates[END_REF] or multi-components BECs [START_REF] Bao | Ground states and dynamics of multi-component Bose-Einstein condensates[END_REF][START_REF] Bao | Ground states of two-component Bose-Einstein condensates with an internal atomic Josephson junction[END_REF][START_REF]Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF]. The direction that we follow here is based on the use of the Fast Fourier Transform (FFT) to accurately discretize the spatial operators, leading therefore to a pseudo SPectral method. Furthermore, the method requires low memory storage. Introduced by Bao et al. in [START_REF] Bao | Efficient and spectrally accurate numerical methods for computing ground and first excited states in Bose-Einstein condensates[END_REF] for non rotating gazes, the application of this method, called Backward Euler SPectral (BESP) method, to rotating BECs has proved to be very accurate over standard finite difference schemes while being extremely efficient [START_REF] Zeng | Efficiently computing vortex lattices in rapid rotating Bose-Einstein condensates[END_REF] (most particularly if one has in mind to develop HPC codes, possibly using speed up GPU computing). However, as we will see in details in Section 3.1, the main drawback of this method is related to the iterative scheme that is used for solving the linear system associated to the BESP scheme. Indeed, as proposed in the literature, the solution which is based on a fixed point approach can be inefficient and can even diverge when the rotation speed is too large. This can be for example the case when considering quadraticquartic potentials and sufficiently large values of Ω. We propose in this paper a robust approach based on the use of Krylov subspace solvers [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF][START_REF] Saad | GMRES -a Generalized Minimal RESidual algorithm for solving nonsymmetric linear systems[END_REF][START_REF] Van Der | BI-CGSTAB -a fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems[END_REF]. In particular, we show that BiCGStab [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF][START_REF] Van Der | BI-CGSTAB -a fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems[END_REF] is particularly robust and very efficient to solve the linear systems. Moreover, we propose to improve the convergence rate of BiCGStab through two simple physics-based preconditioners related to the Laplace and Thomas-Fermi approximations. The resulting method is simple, accurate, efficient and robust. In addition, it can be easily adapted to many situations e.g. multi-components problems, general nonlinearities and potentials with a great potential for an implementation on high performance computers.

The plan of the paper is the following. In Section 2, we derive the BESP scheme for a onecomponent GPE with rotating term. Section 3 is devoted to the numerical solution of the linear systems related to the BESP scheme. We show that the standard successive approximation scheme with relaxation does not converge for large rotations. We propose a robust alternative based on preconditioned Krylov subspace solvers [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. Two preconditioners are considered: one related to the Laplace operator and another one linked to the Thomas-Fermi equation. They appear to be efficient and robust, most particularly with respect to the rotation speed Ω. In Section 4, we provide three additional examples, one in the 2D case, one in the 3D case, and a third one for a two-components system of 2D GPEs. The short Section 5 reports some informations about GPELab (Gross-Pitaevskii Equation Laboratory) which is a freely available Matlab toolbox that in particular proposes methods based on the scheme presented in this paper to compute stationary solutions and dynamics to general multi-components GPEs. Finally, Section 6 concludes.

2 Gradient flow formulation: the BESP discretization

CNGF and Backward Euler time discretization

One standard solution for computing the solution to the minimization problem ( 5)-( 7) is through the projected gradient method which consists in i) computing one step of a gradient method and then ii) project the solution onto the unit sphere S. The method is the so-called imaginary time method usually used in Physics and based on the remark that the real time is replaced by a complex time following t → -it in Eq. ( 1). Let us denote by t 0 = 0 < ... < t n < ... the uniformly spaced discrete times and by ∆t = t n+1 -t n the uniform time step. The Continuous Normalized Gradient Flow (CNGF) [START_REF]Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF][START_REF] Bao | Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow[END_REF][START_REF] Zeng | Efficiently computing vortex lattices in rapid rotating Bose-Einstein condensates[END_REF] is given by

           ∂ t φ = -∇ φ * E β,Ω (φ) = 1 2 ∆φ -V φ -β|φ| 2 φ + ΩL z φ, t n < t < t n+1 , φ(x, t n+1 ) = φ(x, t + n+1 ) = φ(x, t - n+1 ) ||φ(x, t - n+1 )|| 0 , φ(x, 0) = φ 0 (x), x ∈ R d , with ||φ|| 0 = 1. (8)
In the above equations, we set: φ(x, t ± n+1 ) := lim t→t ± n φ(x, t). Hence, time marching corresponds to iterations in the projected gradient. In [START_REF] Bao | Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow[END_REF], the CNGF is proved to be normalization conserving and energy diminishing for β = 0 and a positive potential, in the non rotating case. When t tends towards infinity, φ(x, t) gives an approximation of the steady state solution φ(x) which is a critical point of the energy when the assumption on V ≥ 0 is fulfilled.

Concerning the time discretization of system [START_REF]Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF], the application of the Backward Euler (BE) scheme [START_REF] Bao | Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow[END_REF] leads to the semi-discrete semi-implicit (linear) scheme: for 0 ≤ n ≤ N , we compute the function

φ n+1 such that    A BE,n φ(x) = b BE,n (x), x ∈ R d , φ n+1 (x) = φ(x) || φ|| 0 , (9)
where the operator A BE and the right hand side function b BE are given by

A BE,n := ( I ∆t - 1 2 ∆ + V + β|φ n | 2 -ΩL z ), b BE,n = φ n ∆t . (10)
The maximal time of computation T max = N ∆t is fixed by a stopping criterion (see Eq. ( 34)). Therefore, for a considered physical problem, the number of time steps N is not known a priori and depends on the convergence rate of the iterative scheme to get the ground state solution. The initial function is given by: φ(x, 0) = φ 0 (x), x ∈ R d , with ||φ|| 0 = 1, which is generally chosen as an analytical approximate physical solution (like the Thomas-Fermi approximation).

Spatial discretization: pseudo SPectral scheme based on FFTs

Since the ground state evolution is assumed to be localized in a finite region of the space, we consider that the support of the evolving field is inside a box

O :=] -a x ; a x [×] -a y ; a y [×] - a z ; a z [ (O :=] -a x ; a x [×]
-a y ; a y [ in the 2D case). We now choose the uniform spatial grid:

O J,K,L = {x j,k,ℓ := (x j , y k , z ℓ ); 0 ≤ j ≤ J -1, 0 ≤ k ≤ K -1, 0 ≤ ℓ ≤ L -1}, J
, K, L being three even positive integers. We set:

x j+1 -x j = h x = 2a x /J, y k+1 -y k = h y = 2a y /K, y ℓ+1 -y ℓ = h z = 2a z /L. (11)
Since we assume that φ is compactly supported in O, it satisfies a periodic boundary condition on ∂O (in fact zero) and discrete Fourier transforms can then be used (sine or cosine transforms could also be used according to the boundary condition). The partial Fourier pseudospectral discretizations in the x-, y-and z-directions are respectively given by

φ(x, y, z, t) = 1 J J/2-1 p=-J/2 φp (y, z, t)e iµp(x+ax) , φ(x, y, z, t) = 1 K K/2-1 q=-K/2 φq (x, z, t)e iλq(y+ay) , φ(x, y, z, t) = 1 L L/2-1 r=-L/2 φr (x, y, t)e iξr(z+az) , (12)
with φp , φq , φr respectively the Fourier coefficients in the x-, y-and z-directions

φp (y, z, t) = J-1 j=0 φj (y, z, t)e -iµp(x j +ax) , φq (x, z, t) = K-1 k=0 φk (x, z, t)e -iλq(y k +ay) , φr (x, y, t) = L-1 ℓ=0 φℓ (x, y, t)e -iξr(z ℓ +az) , (13)
and where µ p = πp ax , λ q = πq ay and ξ r = πr az . In the above equations, we set: φj (y, z, t) := φ(x j , y, z, t), φk (x, z, t) := φ(x, y k , z, t) and φℓ (x, y, t) := φ(x, y, z ℓ , t). For the backward Euler scheme, this implies that we have the following spatial approximation

   A BE,n φ = b BE,n , φ n+1 (x) = φ || φ|| 0 , (14)
where φ = ( φ(x j,k,ℓ )) (j,k,ℓ)∈O J,K,L is the discrete unknown vector in C M and the right hand side is b BE,n := φ n /∆t, with φ n = (φ n (x j,k,ℓ )) (j,k,ℓ)∈O J,K,L ∈ C M . For conciseness, let us remark that we do not make the distinction between an array φ in M J×K×L (C) (storage according to the 3D cubic spatial grid) and the corresponding reshaped vector in C M . In the above notation, M J×K×L (C) designates the set of 3D (respectively 2D) arrays with complex coefficients. We also define M = JKL (respectively M = JK) in 3D (respectively in 2D).

The operator A BE,n is given by the map which for any vector

ψ ∈ C M , that is assumed to approximate (ψ(x j,k,ℓ )) ∈ C M for a function ψ, computes a vector Ψ ∈ C M such that Ψ := A BE,n ψ = A BE,n TF ψ + A BE ∆,Ω ψ, A BE,n TF ψ := I ∆t + V + β[[|φ n | 2 ]] ψ, A BE ∆,Ω ψ := - 1 2 [[∆]] -ΩL z ψ. ( 15
)
The evaluation of the different operators is made as follows. For A BE,n TF , the application is direct since it is realized pointwize in the physical space by setting

I j,k,ℓ := δ j,k,ℓ , V j,k,ℓ := V (x j,k,ℓ ), [[|ψ n | 2 ]] j,k,ℓ = |ψ n | 2 (x j,k,ℓ ), (16)
for (j, k, ℓ) ∈ O J,K,L . The symbol δ j,k,ℓ denotes the Dirac delta symbol which is equal to 1 if and only if j = k = ℓ and 0 otherwise. Let us note that the discrete operator A BE,n TF is represented by a diagonal matrix after reshaping. The label "TF" refers to the fact that this operator corresponds to the discretization of the Thomas-Fermi approximation.

For the operator A BE ∆,Ω , we use the three following expressions, for

(j, k, ℓ) ∈ O J,K,L , (- 1 2 ∂ 2 x -iΩy k ∂ x )ψ(x j,k,ℓ ) ≈ (- 1 2 [[∂ 2 x ]] -iΩy k [[∂ x ]])ψ j,k,ℓ := 1 J J/2-1 p=-J/2 ( µ 2 p 2 -Ωy k µ p ) (ψ k,ℓ ) p e iµp(x j +ax) , (- 1 2 ∂ 2 y + iΩx j ∂ y )ψ(x j,k,ℓ ) ≈ (- 1 2 [[∂ 2 y ]] + iΩx j [[∂ y ]])ψ j,k,ℓ := 1 K K/2-1 q=-K/2 ( λ 2 q 2 + Ωx j λ q ) (ψ j,ℓ ) q e iλq(y k +ay) , (- 1 2 ∂ 2 z )ψ(x j,k,ℓ ) ≈ (- 1 2 [[∂ 2 z ]])ψ j,k,ℓ := 1 L L/2-1 r=-L/2 ξ 2 r 2 (ψ j,k ) r e iξr(z ℓ +az) , (17)
and we define the discrete operators

[[∆]] j,k,ℓ := [[∂ 2 x ]] j,k,ℓ + [[∂ 2 y ]] j,k,ℓ + [[∂ 2 z ]] j,k,ℓ , (L z ) j,k,ℓ := -i(x j [[∂ y ]] j,k,ℓ -y k [[∂ x ]] j,k,ℓ ). ( 18
)
The discrete operator

[[∆]] is diagonal in the Fourier space but not L z . Finally, the discrete L 2 -norm || • || 0 is given by ∀φ ∈ C M , ||φ|| 0 := h 1/2 x h 1/2 y h 1/2 z ( (j,k,ℓ)∈O J,K,L |φ j,k,ℓ | 2 ) 1/2 . ( 19
)
3 Robust and efficient preconditioned iterative solvers for BESP

Fixed point approach: limitation for high rotations

Following [START_REF] Bao | Efficient and spectrally accurate numerical methods for computing ground and first excited states in Bose-Einstein condensates[END_REF] and for Ω = 0, the most direct way to solve the linear system in ( 14) is to use a fixed point approach with stabilization parameter ω. However, we will see that this method is not robust since it fails to converge for a high rotation and a stiff nonlinearity. Let us introduce the following operators

A BE ∆,ω = I ∆t - 1 2 [[∆]] + ωI, A BE,n Ω,TF,ω = ΩL z - 1 2 V - 1 2 β[[|φ n | 2 ]] -ωI. (20)
Since the Laplacian operator appearing in A BE ∆,ω is diagonalizable in the Fourier basis and can therefore be directly inverted, a natural method, proposed in [START_REF] Bao | Efficient and spectrally accurate numerical methods for computing ground and first excited states in Bose-Einstein condensates[END_REF] for [START_REF] Baye | Resolution of the Gross-Pitaevskii equation with the imaginarytime method on a Lagrange mesh[END_REF] with Ω = 0, is to compute

the sequence of iterates ( φ(m+1) ) m∈N through    φ(0) = φ n (x), φ(m+1) = A BE ∆,ω -1 A BE,n Ω,TF,ω φ(m) + b BE,n , (21)
to get the solution φ (in fact, an approximation) of the first equation of system [START_REF] Baye | Resolution of the Gross-Pitaevskii equation with the imaginarytime method on a Lagrange mesh[END_REF] for a sufficiently large index m. Since our method is supposed to be spectrally accurate, we need to fix a strong stopping criterion

φ(m+1) - φ(m) ∞ := max (j,k,ℓ)∈O J,K,L φ(m+1) j,k,ℓ - φ(m) j,k,ℓ ≤ ε, (22)
with ε very small (e.g. 10 -12 ). In [START_REF] Bao | Efficient and spectrally accurate numerical methods for computing ground and first excited states in Bose-Einstein condensates[END_REF], the authors prove that the optimal stabilization parameter ω * that minimizes the spectral radius of the iteration matrix (A BE ∆,ω ) -1 A BE,n Ω,TF,ω , is given by

ω * = b max + b min 2 , (23)
where

b max = max (j,k,ℓ)∈O J,K,L 1 2 V j,k,ℓ + 1 2 β[[|φ n | 2 ]] j,k,ℓ (24)
and b min = min

(j,k,ℓ)∈O J,K,L 1 2 V j,k,ℓ + 1 2 β[[|φ n | 2 ]] j,k,ℓ . (25)
The convergence proof of the fixed point method is based on the standard argument that if ρ ω (β, 0) < 1, then the method converges for any β ≥ 0. Denoting by ρ(A) the spectral radius of a matrix A, we define ρ ω (β, Ω) := ρ((A BE ∆,ω ) -1 A BE,n Ω,TF,ω ). This method has been applied in [START_REF] Zeng | Efficiently computing vortex lattices in rapid rotating Bose-Einstein condensates[END_REF] in the non rotating case Ω > 0. However, no proof of convergence is given and only numerical simulations are available.

Let us now numerically illustrate the limitations of the fixed point approach when a rotating term is included. We consider a 2D case for the quadratic-quartic potential

V (x) = (1 -α) x 2 + κ x 4 .
We take the set of parameters α = 1.2 and κ = 0.3 [START_REF] Zeng | Efficiently computing vortex lattices in rapid rotating Bose-Einstein condensates[END_REF]; x is the usual euclidian norm of a vector x ∈ R d . The time step is ∆t = 10 -2 and the square domain is

O =] -15; 15[ 2 , with J = K = 2 9 . The stopping criterion is: either || φ(m+1) - φ(m) || ∞ ≤ 10 -12
or m ≥ 5000. The initial data φ 0 of BESP is given by the Thomas-Fermi approximation

φ 0 (x) := φ TF β (x) = (µ TF β -V (x))/β, if µ TF > V (x), 0, otherwise, (26)
with µ TF β = (4βγ x γ y /π) 1/2 /2 (with γ x = γ y = 1 here), if β = 0, or by the normalized gaussian

φ 0 (x) = (γ x γ y ) 1/4 √ π e -(γxx 2 +γyy 2 )/2 (27)
otherwise. Our example consists in testing numerically the convergence when we iteratively solve the linear system in ( 14) from n = 0 to n = 1 by using the scheme [START_REF] Cerimele | Numerical solution of the Gross-Pitaevskii equation using an explicit finite-difference scheme: An application to trapped Bose-Einstein condensates[END_REF]. We first report on Figure 1(a) the number of iterations #iter of the method (for ω * ) to reach the convergence with respect to β and Ω. As we can see, many iterations are generally required to converge with high precision. Most particularly, for large enough values of Ω (Ω ≥ 1.4 in this example), we observe the divergence of the algorithm, even for small nonlinearities, since #iter = 5000. To clarify the convergence problem, we report on Figure 1(b) the variations of ρ ω * (10 3 , Ω) versus Ω. We can observe that, for a moderate rotating speed Ω (≤ 1.4), then the spectral radius is slightly smaller than 1, implying a slow but converging method. However, the spectral radius ρ ω * (10 3 , Ω) becomes larger than 1 after Ω = 1.4 which makes the method diverges. As a conclusion, even if the iterative method works for small Ω, then the convergence is slow. Furthermore, the method is not robust if sufficiently fast rotations are considered since convergence can break down.

Robust preconditioned Krylov solvers for BESP

Accelerators-preconditioners

For complex problems, it is known that Jacobi, Gauss-Seidel and SOR iterative methods [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF] can suffer from a lack of robustness as this is the case above. More robust iterative solvers are however available like for example for Krylov subspace solvers [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. These methods are well-adapted to efficiently solve large linear systems of the form

Aφ = b, ( 28
)
where A is a matrix from M M (C) (a M × M matrix with complex-valued coefficients) and b ∈ C M . For example, successful solvers are the Conjugate Gradient Squared method (CGS) [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF], the BiConjugate Gradient Stabilized method (BiCGStab) [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF][START_REF] Van Der | BI-CGSTAB -a fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems[END_REF] and the Generalized Minimal RESidual method (GMRES) [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF][START_REF] Saad | GMRES -a Generalized Minimal RESidual algorithm for solving nonsymmetric linear systems[END_REF]. The two first methods rely on the minimization of the error

(φ * -φ (m) ) t A(φ * -φ (m) )
, where φ * is the solution of system (28) and φ (m) is an approximation of φ * in the Krylov subspace K m . These methods can be viewed as the minimization of a quadratic functional in terms of φ (m) , using a "search direction" and "step length" method. The third method is based on the minimization of the residual Aφ (m)b , for φ (m) ∈ K m , using the Arnoldi iteration. All these iterative methods are often called accelerators. If #iter is the number of iterations to reach the solution with a required tolerance ε, then the global computational cost is given by C global = #iter × C iter A , where C iter A is the computational effort for one iteration. This is directly related to the evaluation of matrix-vector products Ψ = Aψ, where ψ is given.

The use of an accelerator alone is generally not sufficient [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. Indeed, the convergence rate of the solver e.g. #iter is guided by the way the eigenvalues of the operator spread out in the complex plane. Clustering is often an expected property to get fast convergence and a small number of iterations to converge. This property cannot hold for the matrices arising in system [START_REF] Baye | Resolution of the Gross-Pitaevskii equation with the imaginarytime method on a Lagrange mesh[END_REF] since they are related to the discretization of second-order elliptic operators. To improve the convergence of the solver and therefore to decrease #iter, one has to precondition the linear system. Preconditioners are generally built in an algebraic way as an approximation of the inverse of A in Eq. ( 28). Among the most robust preconditioners, let us for example mention ILUT methods, SPAI, multigrid techniques [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. The restriction with these approaches is that they require the access to the entries of A. In our case, this cannot be expected since the partial differential operators are efficiently evaluated through FFTs (see Eq. [START_REF] Bretin | Fast rotation of a Bose-Einstein condensate[END_REF]). An alternative is to rather approximate the original partial differential operator involved in the equation by another operator which can be easily inverted. We also have to take into account the computational cost of the application of such preconditioners into the global scheme. For the BESP scheme, the operator that must be solved is given by

A BE,n = I ∆t - 1 2 [[∆]] -ΩL z + V + β[[|φ n | 2 ]]. (29)
When one wants to evaluate Ψ = A BE,n ψ, for ψ ∈ C M , FFT/iFFTs must be used. Following [START_REF] Boyd | Chebyshev and Fourier Spectral Methods[END_REF], the computational cost of a FFT/iFFT for a complex-valued vector ψ is C FFT/iFFT = 5M log 2 (M ). Therefore, the evaluation of A BE,n ψ requires one FFT and two iFFTs resulting in a global cost C iter A BE,n = 3C FFT for the application of the unpreconditioned operator.

A first possibility to precondition our system is to use the operator related to the fixed point approach. This results in the following preconditioned linear system to solve

I + P BE ∆ A BE,n Ω,TF φ = P BE ∆ b BE,n , (30)
where

P BE ∆ = I ∆t - 1 2 [[∆]] -1 and A BE,n Ω,TF = -ΩL z + V + β[[|φ n | 2 ]]. (31)
In the sequel, the preconditioner "P BE ∆ " is called "Laplace" (∆) preconditioner and corresponds to the discretization of the linear heat equation (diffusion part of the operator). Let us remark that a matrix-vector evaluation for the preconditioned operator equation has a global cost C iter ∆ = 4C FFT . A second simple preconditioner can be obtained by only keeping the explicit diagonal terms of the original operator without partial differential operator to avoid any FFT computation and to get an explicit inversion. This gives us the second preconditioned system

I + P BE,n TF A BE ∆,Ω φ = P BE,n TF b BE,n , (32)
where

P BE,n TF = I ∆t + V + β[[|φ n | 2 ]] -1 and A BE,n ∆,Ω = - 1 2 [[∆]] -ΩL z . (33)
All along the paper, the preconditioner P BE,n TF is called "Thomas-Fermi" (TF) preconditioner and is associated with the BE discretization of the Thomas-Fermi equation without rotation. An analysis of the computational cost of a matrix-vector evaluation for this new equation shows that we obtain: C iter TF = 3C FFT , which is the same as for the equation without preconditioner.

Numerical comparison of the preconditioned solvers for BESP

We first compare the different following accelerators: CGS, BiCGStab and GMRES without restart. The tolerance of the Krylov solvers is fixed to ε := 10 -13 all along this section. We consider the same situation and parameters as in subsection 3.1 for the fixed point approach. Let us choose BiCGStab without preconditioner. We solve the BESP system to get the solution from n = 0 to n = 1. We report on Figure 2(a) the number of iterations #iter vs. (β, Ω) for BiCGStab to converge. Compared with the fixed point approach, we can first observe that the method always converges even for large Ω and is therefore robust. Furthermore, the number of iterations is always much smaller than for the fixed point approach to get a very small residual and slightly grows with respect to both increasing values of β and Ω. Hence, the method is efficient. We now compare on Figures 3 the efficiency of the Krylov solvers: GMRES, CGS and BiCGStab. We fix β = 2000 and consider the same situation as above. We report the number of iterations #iter and the CPU time (in seconds) required for the solver to converge. We can observe that the most efficient Krylov solver is BiCGStab. For Ω = 0, CGS fails to converge, thus making it a non robust solver.

We now select the BiCGStab accelerator and compare the efficiency of the Laplace (∆-BiCGStab) and Thomas-Fermi (TF-BiCGStab) preconditioners for the same problem. As we can observe on Figure 4(a), the number of iterations #iter is strongly reduced when one considers the Thomas-Fermi preconditioner and is stable according to Ω. The effect of the Laplace preconditioner is less impressive most particularly for increasing values of Ω. The impact on the CPU time can be directly seen on Figure 4(b). We now fix β = 2000 and Ω = 2. For computing the stationary solution, the global time stopping criterion is fixed (all along the paper) by

φ n+1 -φ n ∞ ≤ ǫ∆t, (34)
with ǫ = 10 -6 . According to [START_REF] Raman | Vortex nucleation in a stirred Bose-Einstein condensate[END_REF], one obtains the computed ground state solution given on Figure 2(b). We now solve the linear systems with and without preconditioner for the first 2 × 10 4 time steps to see the behavior of the solvers over a long time interval. As we can see on Figure 5(a), the number of iterations per time step is strongly reduced all along the computations for the TF-BiCGStab solver. The preconditioner ∆-BiCGStab is also helpful but less effective. Furthermore, the reduction of iterations is quite stable over the time interval even if we can observe that less iterations are needed for the first time steps (see Figure 5(b)). We only retain TF-BiCGStab which has proved to be the most efficient and robust solver. Let us consider the computational domain O =] -15, 15[ 2 with J = K = 2 9 . The time step is fixed to ∆t = 10 -2 and the stopping criterion for BiCGStab is given by formula [START_REF] Chiofalo | Ground state of trapped interacting Bose-Einstein condensates by an explicit imaginary-time algorithm[END_REF] to ε = 10 -6 . Let us analyze the behavior of the preconditioned solver with respect to the nonlinearity strength β. For Ω = 2, we consider five increasing values of β: 100, 500, 1000, 2000 and 5000. Figures 6 report the number of iterations #iter with respect to the imaginary time for the first 10 4 time steps. We observe that the convergence rate is not too strongly affected by the nonlinearity, even for large values of β. Only the first time steps show a different number of iterations according to the increasing nonlinearity strength. We can conclude that the TF-BiCGStab preconditioner which is designed for strong nonlinear problems is indeed robust. Let us analyze the convergence properties with respect to the rotation speed Ω. We fix β = 2000 and consider Ω = 0, 1, 2, 3 and 3.5. As we can see on Figures 7, the number of iterations increases with respect to Ω. This is due to the fact that the effect of the rotation is not included in the Thomas-Fermi preconditioner.

We now consider the influence of the time step ∆t on the total number of iterations needed to reach a ground state. We fix β = 300 and Ω = 0.5. We launch the BESP scheme for different values of ∆t: 0.5, 0.1, 5 × 10 -2 , 10 -2 , 5 × 10 -3 , 10 -3 and 5 × 10 -4 , for J = K = 2 9 , a stopping criterion ε = 10 -6 for BiCGStab and criterion [START_REF] Raman | Vortex nucleation in a stirred Bose-Einstein condensate[END_REF] for the global iterative time scheme related to the CNGF. As we can see on Figures 8, the total number of iterations to reach the same solution is inversely proportional to the time step. An explanation to this phenomenon is that, when the solution reaches an almost steady state, taking a larger time step is more adapted to attain the global minimum of the energy functional. However, in cases where the physical parameters β and Ω are larger, a smaller time step must be chosen since the minimization algorithm would otherwise only lead to a local minimum.

Ω = 0 Ω = 1 Ω = 2 Ω = 2.5 Ω = 3

Three additional examples

From the previous section, we have seen that the TF-BiCGStab preconditioned solver is efficient and robust, most particularly for strong nonlinearities β, strong potentials and large rotation speeds Ω. We propose to investigate the properties of the preconditioned iterative solvers 1) for a problem with a quadratic potential, 2) for a three-dimensional GPE and 3) for a two components BEC.

A problem with a harmonic potential

Here we consider the computation of the stationary state of the GPE with a quadratic potential: V (x) = 1 2 x 2 and a rotation speed ranging from 0 to 0.95 (1 is the critical velocity here). The discretization/approximation parameters are the same as in subsection 3.1. First, let us consider BiCGStab with a tolerance ε equal to 10 -13 and solve the BESP system from n = 0 to n = 1. The initial data is given by the Thomas-Fermi approximation [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF] when β > 0. For β = 0, we consider the normalized gaussian [START_REF] Saad | GMRES -a Generalized Minimal RESidual algorithm for solving nonsymmetric linear systems[END_REF]. From Figure 9(a), we can see that BiCGStab is robust and efficient. Additional numerical simulations show that the unpreconditioned CGS and BiCGStab lead to about the same convergence rate and outperform GMRES. Let us fix β = 2000. On Figure 9(b), we can see that, unlike the previous case, the Thomas-Fermi preconditioner does not reduce #iter while the Laplace preconditioner is extremely efficient. This reduction is furthermore almost not affected by the value of Ω.

A three-dimensional example

We consider now the 3D GPE Eq. ( 1) with potential V (x) = (1 -α) x 2 + κ x 4 (α = 1.2 and κ = 0.3). The time step is ∆t = 10 -2 and the computational domain is ] -30; 30[ 3 , for a spatial discretization with J = K = L = 2 7 . The tolerance ε of the Krylov subspace solvers is 10 -9 (all along this section). The initial data φ 0 for BESP is the Thomas-Fermi approximation

φ TF β (x) = (µ TF β -V (x))/β, if µ TF > V (x), 0, otherwise, (35)
with µ TF β = (15βγ x γ y γ z /(4π)) 2/5 /2 (with γ x = γ y = γ z = 1 here), if β = 0, or the normalized gaussian

φ 0 (x) = (γ x γ y γ z ) 1/4 π 3/4 e -(γxx 2 +γyy 2 +γzz 2 )/2 (36)
otherwise. We report on Figure 10 the number of iterations #iter vs. (β, Ω) for BiCGStab to converge, from t 0 = 0 to t 1 = ∆t. Like the 2D case with quadratic-quartic potential, the method is robust and efficient. Let us now compare the efficiency of the Krylov solvers: GMRES, CGS and BiCGStab. We fix β = 2000 and consider the same situation as above. We report on Figure 11(a) the number of iterations #iter to converge. We conclude that the most efficient solver is BiCGStab.

We now analyze the performance of the preconditioned and unpreconditioned BiCGStab. As we can see on Figure 11(b), the Thomas-Fermi preconditioner reduces a lot the number of iterations #iter. Moreover, this reduction is not affected by the value of Ω. In addition, the effect of the Laplace preconditioner does not provide such a reduction.

A two-components Bose-Einstein condensate

We consider now a two-components time-dependent GPE with rotating term and coupled nonlinearities

[8]        i∂ t ψ 1 (t, x) = - 1 2 ∆ + V (x) -Ω • L ψ 1 (t, x) + βf 1 (|ψ 1 (t, x)|, |ψ 2 (t, x)|)ψ 1 (t, x) -ψ 2 (t, x), i∂ t ψ 2 (t, x) = - 1 2 ∆ + V (x) -Ω • L ψ 2 (t, x) + βf 2 (|ψ 2 (t, x)|, |ψ 1 (t, x)|)ψ 2 (t, x) -ψ 1 (t, x), (37)
for x ∈ R d , t > 0. Functions ψ 1 and ψ 2 are the two components wave functions of the condensate. In our case, we have the following coupled nonlinearities [START_REF]Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF]

f 1 (|ψ 1 (t, x)|, |ψ 2 (t, x)|) = β 1,1 |ψ 1 (t, x)| 2 + β 1,2 |ψ 2 (t, x)| 2
f 2 (|ψ 2 (t, x)|, |ψ 1 (t, x)|) = β 2,2 |ψ 2 (t, x)| 2 + β 2,1 |ψ 1 (t, x)| 2 ,
where β 1,1 = 1, β 2,2 = 0.97 and β 1,2 = β 2,1 = 0.94. Setting Ψ = (ψ 1 , ψ 2 ), the energy of this system is defined by

E β,Ω (Ψ) = 2 j=1 R d 1 2 |∇ψ j | 2 + V (x)|ψ j (x)| 2 + β β j,j 2 |ψ j | 4 -Ωℜ ψ * j L z ψ j dx + R d β β 1,2 + β 2,1 2 |ψ 1 (x)| 2 |ψ 2 (x)| 2 dx - R d 2ℜ(ψ 1 (x)ψ 2 (x))dx. (38)
The mass conservation property writes down

||Ψ|| 2 0 = R d (|ψ 1 (x)| 2 + |ψ 2 (x)| 2 )dx = 1. (39)
To compute a stationary state of this system, we adapt the CNGF method [START_REF] Bao | Ground states and dynamics of multi-component Bose-Einstein condensates[END_REF][START_REF]Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF] which is energy diminishing. Being given a uniform time step ∆t = t n+1 -t n , ∀n ∈ N, the CNGF formulation leads to computing Φ

= (φ 1 , φ 2 ) solution to                  ∂ t φ 1 = -∇ φ * 1 E β,Ω (Φ) = 1 2 ∆φ 1 -V φ 1 -βf 1 (|φ 1 |, |φ 2 |)φ 1 + ΩL z φ 1 + φ 2 , t n < t < t n+1 , ∂ t φ 2 = -∇ φ * 2 E β,Ω (Φ) = 1 2 ∆φ 2 -V φ 2 -βf 2 (|φ 2 |, |φ 1 |)φ 2 + ΩL z φ 2 + φ 1 , t n < t < t n+1 , Φ(x, t n+1 ) = Φ(x, t + n+1 ) = Φ(x, t - n+1 ) ||Φ(x, t - n+1 )|| 0 , Φ(x, 0) = Φ 0 (x), x ∈ R d , with ||Φ|| 0 = 1. (40)
We assume that the two components evolve in O :=] -a x ; a x [×] -a y ; a y [×] -a z ; a z [ (for the 3D case). Using the same notations as in the one-component case, the BESP scheme leads to solving

   A BE,n Φ = b BE,n , Φ n+1 = Φ || Φ|| 0 , (41)
where Φ = ( φ1 (x j,k,ℓ ), φ2 (x j,k,ℓ )) (j,k,ℓ)∈O J,K,L is a discrete vector in C 2M . If Φ n = (φ n 1 (x j,k,ℓ ), φ n 2 (x j,k,ℓ )) (j,k,ℓ)∈O J,K,L in C 2M is an approximation of the solution Φ at time t n on the spatial grid O, then the right hand side is given by b BE := Φ n /∆t. The operator A BE,n is defined by the map which for any vector Φ

= (φ 1 , φ 2 ) ∈ C 2M associates Ψ ∈ C 2M such that Ψ := A BE,n Φ = A BE,n TF Φ + A BE ∆,Ω Φ, A BE,n TF Φ := 1 ∆t I 0 0 I + V -1 -1 V + β F 1 (Φ n ) 0 0 F 2 (Φ n ) Φ, A BE ∆,Ω Φ := - 1 2 [[∆]] 0 0 [[∆]] -Ω L z 0 0 L z Φ, (42)
where

F 1 (Φ n ) = β 1,1 [[|φ n 1 | 2 ]] + β 1,2 [[|φ n 2 | 2 ]], F 2 (Φ n ) = β 2,2 [[|φ n 2 | 2 ]] + β 2,1 [[|φ n 1 | 2 ]]. (43)
The operators above are evaluated by extending the strategy developed in the one-component case. Therefore, Krylov subspace solvers can be directly applied here. Concerning the preconditioners, the extension of the "Laplace" preconditioner gives

P BE ∆ = 1 ∆t I 0 0 I - 1 2 [[∆]] 0 0 [[∆]] -1 . (44)
For the "Thomas-Fermi" preconditioner, we have the diagonal "Thomas-Fermi" preconditioner (TF-Diag in Figure 12(b))

P BE,n TF,Diag = 1 ∆t

I 0 0 I + V 0 0 V + β F 1 (Φ n ) 0 0 F 2 (Φ n ) -1 (45)
and the full "Thomas-Fermi" (TF-Full in Figure 12(b)) preconditioner P BE,n TF,Full = 1 ∆t

I 0 0 I + V -1 -1 V + β F 1 (Φ n ) 0 0 F 2 (Φ n ) -1 . ( 46
)
To analyze the algorithm, we consider the quadratic potential V (x) = 1 2 |x| 2 . In the CNGF, we choose the time step ∆t = 10 -2 . The square domain is ] -15, 15[ 2 , with J = K = 2 9 . The tolerance ε of the Krylov solvers is set to 10 -13 . The initial data Φ 0 is: Φ 0 = (φ 0 , φ 0 ), where φ 0 is either given by the Thomas-Fermi approximation [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF] when β > 0, or by the normalized gaussian [START_REF] Saad | GMRES -a Generalized Minimal RESidual algorithm for solving nonsymmetric linear systems[END_REF] when β = 0. On Figure 12(a), we report the number of iterations #iter vs. (β, Ω) for BiCGStab to converge, (steps n = 0 to n = 1). Like in the one-component case (see Section 4.1), we observe that BiCGStab is robust and efficient. Extensive numerical simulations show that BiCGStab has the best convergence rate compared to GMRES and CGS while being robust. Let us fix β = 2000. We report on Figure 12(b) #iter for the different preconditioned BiCGStab. We see that the "Laplace" preconditioner leads to the best performance. There is no gain when considering the diagonal and full TF preconditioners.

GPELab: an associated Matlab toolbox

The methods presented in the paper are the basic methods implemented in the freely available Gross-Pitaevskii Equation Laboratory * (GPELab) Matlab toolbox. The aim of this collection of algorithms is in particular to compute ground states of multi-components BECs with rotating terms * http://www.iecn.u-nancy.fr/~duboscq/GPELab.html and general nonlinearities based on the preconditioned BESP scheme. Therefore, all the test cases reported here can be directly reproduced by using the associated Matlab scripts [START_REF] Antoine | GPELab: A Matlab toolbox for computing stationary solutions and dynamics of Gross-Pitaevskii Equations (GPE), User-guide[END_REF]. Furthermore, GPELab allows to compute the (stochastic) dynamics of BECs based on GPEs by using a ReSP [START_REF] Besse | A relaxation scheme for the nonlinear Schrödinger equation[END_REF] or second-and fourth-orders time splitting spectral schemes [START_REF]Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF][START_REF] Bao | An efficient and spectrally accurate numerical method for computing dynamics of rotating Bose-Einstein condensates[END_REF] (that are explicit and so do not need any linear system solution).

Conclusion

We presented an accurate, efficient and robust method based on BESP for computing the stationary states of strongly nonlinear and fast rotating GPEs. The method combines a Backward Euler scheme in time and FFT-based discretization in space. The associated linear systems are solved by the BiCGStab solver in conjunction with a Laplace or Thomas-Fermi preconditioner. The method is supported by a full numerical study with various one-and two-components examples in 2D and 3D. Finally, an associated Matlab code, GPELab [START_REF] Antoine | GPELab: A Matlab toolbox for computing stationary solutions and dynamics of Gross-Pitaevskii Equations (GPE), User-guide[END_REF], can be freely downloaded for testing the method for a large class of problems related to the GPE.

Considering a Relaxation scheme, we can also extend [START_REF] Duboscq | Analyse et Simulation d'Equations de Schrödinger Déterministes et Stochastiques[END_REF] the method for computing the dynamics of vortices [START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrodinger/Gross-Pitaevskii equations[END_REF][START_REF] Zhang | Numerical Study of Vortex Interactions in Bose-Einstein Condensation[END_REF]. It appears that the preconditioned CGS is the most robust and efficient iterative method in the case of the ReSP scheme.

Figure 1 :

 1 Figure1: 2D quadratic-quartic potential: convergence/divergence of the iterative scheme[START_REF] Cerimele | Numerical solution of the Gross-Pitaevskii equation using an explicit finite-difference scheme: An application to trapped Bose-Einstein condensates[END_REF] with respect to Ω.

  #iter vs. (β, Ω) for BICGStab without preconditioner to converge. (b) Converged ground state for Ω = 2 and β = 2000.

Figure 2 :Figure 3 :

 23 Figure 2: 2D quadratic-quartic potential: #iter vs. (β, Ω) for BiCGStab without preconditioner (left) and converged solution for β = 2000 (right).

Figure 4 :

 4 Figure 4: 2D quadratic-quartic potential: #iter (left) and CPU time (right) vs. Ω for the first time step of BiCGStab with and without preconditioner.

Figure 5 :

 5 Figure 5: 2D quadratic-quartic potential: #iter for the first 2 × 10 4 time steps (left) and 6 × 10 2 time steps (right) vs. the imaginary time by using BiCGStab with and without preconditioner.

  #iter vs. imaginary time for 10 4 time steps #iter vs. imaginary time for 3 × 10 2 time steps

Figure 6 :

 6 Figure 6: 2D quadratic-quartic potential and TF-BiCGStab: #iter vs. different values of β for the first 10 4 (left) and 300 (right) time steps (Ω = 2).

  (a) #iter vs. imaginary time for 10 4 time steps #iter vs. imaginary time for 3 × 10 2 time steps

Figure 7 :

 7 Figure 7: 2D quadratic-quartic potential and TF-BiCGStab: #iter vs. different values of the rotation speed Ω for the first 10 4 (left) and 300 (right) time steps (β = 2000).

  ∆ t = 0.01 ∆ t = 0.005 ∆ t = 0.001 ∆ t = 0.0005 (b) #iter vs. imaginary time

Figure 8 :

 8 Figure 8: 2D quadratic-quartic potential and TF-BiCGStab: total #iter vs. different time steps ∆t (left) and #iter vs. imaginary time for a complete computation with different ∆t (right).

  #iter vs. (β, Ω) (without preconditioner) #iter vs. Ω (with preconditioner)

Figure 9 :

 9 Figure 9: 2D harmonic potential: #iter vs. (β, Ω) for BiCGStab without preconditioner (left) and #iter vs. Ω for the preconditioned BiCGStab (right).

Figure 10 :

 10 Figure 10: 3D quadratic-quartic potential: #iter vs. (β, Ω) for BICGStab without preconditioner to converge.

  #iter vs. Ω (without preconditioner) #iter vs. Ω (with preconditioner)

Figure 11 :

 11 Figure 11: 3D quadratic-quartic potential: #iter vs. Ω for the first time step of the Krylov solvers without (left) and with (right) preconditioner.

  #iter vs. Ω (without preconditioner) #iter vs. Ω (with preconditioner)

Figure 12 :

 12 Figure 12: 2D harmonic potential for the 2-components case: #iter vs. (β, Ω) for BiCGStab (without preconditioner) and #iter vs. Ω for the preconditioned BiCGStab (right) to converge.

Acknowledgements. This work was partially supported by the French ANR grants MicroWave NT09_460489 ("Programme Blanc" call) and ANR-12-MONU-0007-02 BECASIM ("Modèles Numériques" call).