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

Antoine Levitt
email: antoine.levitt@inria.fr

Qinglin Tang

Efficient spectral computation of the stationary states of rotating Bose-Einstein condensates by the preconditioned nonlinear conjugate gradient method

Keywords: Bose-Einstein condensation, rotating Gross-Pitaevskii equation, stationary states, Fourier spectral method, nonlinear conjugate gradient, optimization algorithms on Riemannian manifolds, preconditioner

We propose a preconditioned nonlinear conjugate gradient method coupled with a spectral spatial discretization scheme for computing the ground states (GS) of rotating Bose-Einstein condensates (BEC), modeled by the Gross-Pitaevskii Equation (GPE). We first start by reviewing the classical gradient flow (also known as imaginary time (IMT)) method which considers the problem from the PDE standpoint, leading to numerically solve a dissipative equation. Based on this IMT equation, we analyze the forward Euler (FE), Crank-Nicolson (CN) and the classical backward Euler (BE) schemes for linear problems and recognize classical power iterations, allowing us to derive convergence rates. By considering the alternative point of view of minimization problems, we propose the preconditioned gradient (PG) and conjugate gradient (PCG) methods for the GS computation of the GPE. We investigate the choice of the preconditioner, which plays a key role in the acceleration of the convergence process. The performance of the new algorithms is tested in 1D, 2D and 3D. We conclude that the PCG method outperforms all the previous methods, most particularly for 2D and 3D fast rotating BECs, while being simple to implement.

Introduction

Bose-Einstein Condensates (BECs) were first predicted theoretically by S.N. Bose and A. Einstein, before being realized experimentally 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]. This state of matter has the interesting feature that macroscopic quantum physics properties can emerge and be observed in laboratory experiments. The literature on BECs has grown extremely fast over the last 20 years in atomic, molecular, optics and condensed matter physics, and applications from this new physics are starting to appear in quantum computation for instance [START_REF] Byrnes | Macroscopic quantum computation using Bose-Einstein condensates[END_REF]. Among the most important directions, a particular attention has been paid towards the understanding of the nucleation of vortices [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 the properties of dipolar gases [START_REF] Bao | 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] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF]. At temperatures T which are much smaller than the critical temperature T c , the macroscopic behavior of a BEC can be well described by a condensate wave function ψ which is solution to a Gross-Pitaevskii Equation (GPE). Being able to compute efficiently the numerical solution of such a class of equations is therefore extremely useful. Among the most crucial questions are the calculations of stationary states, i.e. ground/excited states, and of the real-time dynamics [START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF][START_REF] Antoine | Modeling and Computation of Bose-Einstein Condensates: Stationary States, Nucleation, Dynamics, Stochasticity[END_REF][START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF].

To fully analyze a representative and nontrivial example that can be extended to other more general cases, we consider in this paper a BEC that can be modeled by the rotating (dimensionless) GPE. In this setting, the computation of a ground state of a d-dimensional BEC takes the form of a constrained minimization problem:

find φ ∈ L 2 (R d ) such that φ ∈ arg min φ =1 E(φ),
where E is the associated energy functional. Several approaches can be developed for computing the stationary state solution to the rotating GPE. For example, some techniques are based on appropriate discretizations of the continuous normalized gradient flow/imaginary-time formulation [START_REF] Adhikari | Numerical solution of the two-dimensional Gross-Pitaevskii equation for trapped interacting atoms[END_REF][START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF][START_REF] Antoine | Modeling and Computation of Bose-Einstein Condensates: Stationary States, Nucleation, Dynamics, Stochasticity[END_REF][START_REF] Bao | 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] Baye | Resolution of the Gross-Pitaevskii equation with the imaginary-time 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] Zeng | Efficiently computing vortex lattices in rapid rotating Bose-Einstein condensates[END_REF], leading to various iterative algorithms. These approaches are general and can be applied to many situations (dipolar interactions, multi-components GPEs...). We refer for instance to the recent freely distributed Matlab solver GPELab that provides the stationary states computation [START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations I: Computation of stationary solutions[END_REF] (and real-time dynamics [START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations II: Dynamics and stochastic simulations[END_REF]) for a wide variety of GPEs based on the so-called BESP (Backward Euler pseu-doSPectral) scheme [START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF][START_REF] Antoine | Modeling and Computation of Bose-Einstein Condensates: Stationary States, Nucleation, Dynamics, Stochasticity[END_REF][START_REF] Bao | 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] (see also Sections 4 and 6). Other methods are related to the numerical solution of the nonlinear eigenvalue problem [START_REF] Dion | Ground state of the time-independent 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 on 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] 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 of the Gross-Pitaevskii energy with rotation[END_REF]. As we will see below in Section 4, some connections exist between these approaches. Finally, a regularized Newton-type method was proposed recently in [START_REF] Wu | A regularized Newton method for computing ground states of Bose-Einstein condensates[END_REF].

Optimization problems with orthogonal or normalization constraints also occur in different branches of computational science. An elementary but fundamental example is the case of a quadratic energy, where solving the minimization problem is equivalent to finding an eigenvector associated with the lowest eigenvalue of the symmetric matrix representing the quadratic form. A natural generalization is a class of orthogonalized minimization problems, which for a quadratic energy reduce to finding the N first eigenvectors of a matrix. Many problems in electronic structure theory are of this form, including the celebrated Kohn-Sham and Hartree-Fock models [START_REF] Cances | Computational quantum chemistry: a primer[END_REF][START_REF] Saad | Numerical methods for electronic structure calculations of materials[END_REF]. Correspondingly, a large amount of effort has been devoted to finding efficient discretization and minimization schemes. A workhorse of these approaches is the nonlinear preconditioned conjugate gradient method, developed in the 80s [START_REF] Payne | Iterative minimization techniques for ab initio total-energy calculations: molecular dynamics and conjugate gradients[END_REF], as well as several variants of this (the Davidson algorithm, or the LOBPCG method [START_REF] Knyazev | Toward the optimal preconditioned eigensolver: locally optimal block preconditioned conjugate gradient method[END_REF]).

Although similar, there are significant differences between the mathematical structure of the problem in electronic structure theory and the Gross-Pitaevskii equation. In some respects, solving the GPE is easier: there is only one wavefunction (or only a few for multi-species gases), and the nonlinearity is often local (at least when dipolar effects are not taken into account), with a simpler mathematical form than many electronic structure models. On the other hand, the Gross-Pitaevskii equation describes the formation of vortices: the discretization schemes must represent these very accurately, and the energy landscape presents shallower minima, leading to a more difficult optimization problem.

In the present paper, we consider the constrained nonlinear conjugate gradient method for solving the rotating GPE (Section 2) with a pseudospectral discretization scheme (see Section 3). This approach provides an efficient and robust way to solve the minimization problem. Before introducing the algorithm, we review in Section 4 the discretization of the gradient flow/imaginary-time equation by standard schemes (explicit/implicit Euler and Crank-Nicolson schemes). This enables us to make some interesting and meaningful connections between these approaches and some techniques related to eigenvalue problems, such as the power method. In sections 5.1 and 5.2, we introduce the projected preconditioned gradient (PG) and preconditioned conjugate gradient (PCG) methods for solving the minimization problem on the Riemannian manifold defined by the spherical constraints. In particular, we provide some formulae to compute the stepsize arising in such iterative methods to get the energy decay assumption fulfilled. The stopping criteria and convergence analysis are discussed in sections 5.3 and 5.4. We then investigate the design of preconditioners (section 5.5). In particular, we propose a new simple symmetrical combined preconditioner, denoted by P C . In Section 6, we consider the numerical study of the minimization algorithms for the 1D, 2D and 3D GPEs (without and with rotation). We first propose in section 6.1 a thorough analysis in the one-dimensional case. This shows that the PCG approach with combined preconditioner P C and pseudospectral discretization, called PCG C method, outperforms all the other approaches, most particularly for very large nonlinearities. In sections 6.2 and 6.3, we confirm these properties in 2D and 3D, respectively, and show how the PCG C algorithm behaves with respect to increasing the rotation speed. Finally, section 7 provides a conclusion.

Definitions and notations

The problem we consider is the following:

find φ ∈ L 2 (R d ) such that φ ∈ arg min φ =1 E(φ). (2.1)
We write φ = R d |φ| 2 for the standard L 2 -norm and the energy functional E is defined by

E(φ) = R d 1 2 |∇φ| 2 + V (x)|φ| 2 + η 2 |φ| 4 -ωφ * L z φ ,
where V is an external potential, η is the nonlinearity strength, ω is the rotation speed, and L z = i(y∂ x -x∂ y ) is the angular momentum operator.

A direct computation of the gradient of the energy leads to

∇E(φ) = 2H φ φ, with H φ = - 1 2 ∆ + V + η|φ| 2 -ωL z
the mean-field Hamiltonian. We can compute the second-order derivative as

1 2 ∇ 2 E(φ)[f, f ] = f, H φ f + η Re φ 2 , f 2 .
We introduce S = {φ ∈ L 2 (R d ), φ = 1} as the spherical manifold associated to the normalization constraint. Its tangent space at a point φ ∈ S is

T φ S = {h ∈ L 2 (R d ), Re φ, h = 0}
, and the orthogonal projection M φ onto this space is given by

M φ h = h -Re φ, h φ.
The Euler-Lagrange equation (first-order necessary condition) associated with the problem (2.1) states that, at a minimum φ ∈ S, the projection of the gradient on the tangent space is zero, which is equivalent to

H φ φ = λφ,
where λ = H φ φ, φ is the Lagrange multiplier associated to the spherical constraint, and is also known as the chemical potential. Therefore, the minimization problem can be seen as a nonlinear eigenvalue problem. The second-order necessary condition states that, for all h ∈ T φ S,

1 2 ∇ 2 E(φ)[h, h] -λ h 2 ≥ 0.
For a linear problem (η = 0) and for problems where the nonlinearity has a special structure (for instance, the Hartree-Fock model), this implies that λ is the lowest eigenvalue of H φ (a property known as the Aufbau principle in electronic structure theory). This property is not satisfied here.

Discretization

To find a numerical solution of the minimization problem, the function φ ∈ L 2 (R d ) must be discretized. The presence of vortices in the solution imposes strong constraints on the discretization, which must be accurate enough to resolve fine details. Several discretization schemes have been used to compute the solution to the GPE, including high-order finite difference schemes or finite element schemes with adaptive meshing strategies [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 of the Gross-Pitaevskii energy with rotation[END_REF]. Here, we consider a standard pseudo-spectral discretization based on Fast Fourier Transforms (FFTs) [START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF][START_REF] Antoine | Modeling and Computation of Bose-Einstein Condensates: Stationary States, Nucleation, Dynamics, Stochasticity[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].

We truncate the wave function φ to a square domain [-L, L] d , with periodic boundary conditions, and discretize φ with the same even number of grid points M in any dimension. These two conditions can of course be relaxed to a different domain size L and number of points M in each dimension, at the price of more complex notations. We describe our scheme in 2D, its extension to other dimensions being straightforward. We introduce a uniformly sampled grid:

D M := {x k1,k2 = (x k1 , y k2 )} (k1,k2)∈O M , with O M := {0, . . . , M -1} 2 , x k1+1 -x k1 = y k2+1 -y k2 = h,
with mesh size h = 2L/M , M an even number. We define the discrete Fourier frequencies (ξ p , µ q ), with ξ p = pπ/L, -M/2 ≤ p ≤ M/2 -1, and µ q = qπ/L, -M/2 ≤ q ≤ M/2 -1. The pseudo-spectral approximations φ of the function φ in the xand y-directions are such that

φ(t, x, y) = 1 M M/2-1 p=-M/2 φ p (t, y)e iξp(x+L) , φ(t, x, y) = 1 M M/2-1 q=-M/2 φ q (t, x)e iµq(y+L) ,
where φ p and φ q are respectively the Fourier coefficients in the x-and y-directions

φ p (t, y) = M -1 k1=0 φ k1 (t, y)e -iξp(x k 1 +L) , φ q (t, x) = M -1 k2=0 φ k2 (t, x)e -iµq(y k 2 +L) .
The following notations are used: φ k1 (t, y) = φ(t, x k1 , y) and φ k2 (t, x) = φ(t, x, y k2 ). In order to evaluate the operators, we introduce the matrices

I k1,k2 := δ k1,k2 , [[V ]] k1,k2 := V (x k1,k2 ), [[|φ| 2 ]] k1,k2 = |φ k1,k2 | 2 , for (k 1 , k 2 ) ∈ O M ,
) ∈ O M , ∂ 2 x φ(x k1,k2 ) ≈ ([[∂ 2 x ]] φ) k1,k2 := - 1 M M/2-1 p=-M/2 ξ 2 p ( φ k2 ) p e iξp(x k 1 +L) , ∂ 2 y φ(x k1,k2 ) ≈ ([[∂ 2 y ]] φ) k1,k2 := - 1 M M/2-1 q=-M/2
µ 2 q ( φ k1 ) q e iµq(y k 2 +L) ,

(x∂ y φ)(x k1,k2 ) ≈ (x[[∂ y ]] φ) k1,k2 := 1 M M/2-1 q=-M/2
ix k1 µ q ( φ k1 ) q e iµq(y k 2 +L) ,

(y∂ x φ)(x k1,k2 ) ≈ (y[[∂ x ]] φ) k1,k2 := 1 M M/2-1 p=-M/2 iy k2 ξ p ( φ k2 ) p e iξp(x k 1 +L) . (3.2)
By considering the operators from

C N (N = M 2 (in 2D)) to C given by [[∆]] := [[∂ 2 x ]] + [[∂ 2 y ]] and [[L z ]] := -i(x[[∂ y ]] -y[[∂ x ]]
), we obtain the discretization of the gradient of the energy

∇E(φ) = 2H φ φ, with H φ = - 1 2 [[∆]] + [[V ]] + η[[|φ| 2 ]] -ω[[L z ]].
We set φ := ( φ(x k1,k2 )) (k1,k2)∈O M as the discrete unknown vector in C N . For conciseness, we identify an array φ in the vector space of 2D complex-valued arrays M M (C) (storage according to the 2D grid) and the reshaped vector in C N . Finally, the cost for evaluating the application of a 2D FFT is O(N log N ).

In this discretization, the computations of Section 2 are still valid, with the following differences: φ is an element of C N , the operators ∆, V , |φ| 2 and L z are N × N hermitian matrices, and the inner product is the standard C N inner product. In the following sections, we will assume a discretization like the above one, and drop the brackets in the operators for conciseness.

Review and analysis of classical methods for computing the ground states of GPEs

Once an appropriate discretization is chosen, it remains to compute the solution to the discrete minimization problem

φ ∈ arg min φ∈C N , φ =1 E(φ). (4.3)
Classical methods used to find solutions of (4.3) mainly use the so-called imaginary time equation, which is formally obtained by considering imaginary times in the time-dependent Schrödinger equation. Mathematically, this corresponds to the gradient flow associated with the energy E on the manifold S:

∂ t φ = - 1 2 M φ ∇E(φ) = -(H φ φ -λ(φ)φ). (4.4)
As is well-known, the oscillatory behavior of the eigenmodes of the Schrödinger equation become dampening in imaginary time, thus decreasing the energy. The presence of the Lagrange multiplier λ, coming from the projection on the tangent space on S, ensures the conservation of norm: φ = 1 for all times. This equation can be discretized in time and solved. However, since H φ is an unbounded operator, explicit methods encounter CFL-type conditions that limit the size of their time step, and many authors [START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF][START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations I: Computation of stationary solutions[END_REF][START_REF] Bao | 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][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] use a backward-Euler discretization scheme. We are interested in this section in obtaining the asymptotic convergence rates of various discretizations of this equation, to compare different approaches. To that end, consider here the model problem of finding the first eigenpair of a N × N symmetric matrix H. We label its eigenvalues

λ 1 ≤ λ 2 • • • ≤ λ N ,
and assume that the first eigenvalue λ 1 is simple, so that λ 1 < λ 2 . This model problem is a linearized version of the full nonlinear problem. It is instructive for two reasons: first, any good algorithm for the nonlinear problem must also be a good algorithm for this simplified problem. Second, this model problem allows for a detailed analysis that leads to tractable convergence rates. These convergence rates allow a comparison between different schemes, and are relevant in the asymptotic regime of the nonlinear problem.

We consider the following discretizations of equation (4.4): Forward Euler (FE), Backward Euler (BE), and Crank-Nicolson (CN) schemes:

φFE n+1 -φ n ∆t = -(Hφ n -λ(φ n )φ n ), (4.5)
φBE n+1 -φ n ∆t = -(H φBE n+1 -λ(φ n )φ n ), (4.6)
φCN n+1 -φ n ∆t = - 1 2 (H φCN n+1 -λ(φ n ) φCN n+1 ) - 1 2 (Hφ n -λ(φ n )φ n ). (4.7)
These discretizations all decrease the energy when ∆t > 0 is small enough, but do not preserve the norm: the departure from normalization is of order O(∆t 2 ). Therefore, they are followed by a projection step

φ n+1 = φn+1 φn+1 .
Note that some authors do not include the λ term, choosing instead to work with

φFE n+1 -φ n ∆t = -Hφ n , (4.8) φBE n+1 -φ n ∆t = -H φBE n+1 , (4.9)
φCN n+1 -φ n ∆t = - 1 2 (H φCN n+1 + Hφ n ). (4.10)
These also yield schemes that decrease the energy. However, because of the use of the unprojected gradient Hφ instead of M φ (Hφ) = Hφλφ, the departure from normalization is O(∆t), instead of O(∆t 2 ) for the projected case. The difference between the two approaches is illustrated in Figure 1. φn+1 φn+1 For the FE and BE methods, a simple algebraic manipulation shows that one step of the method with the λ term is equivalent to one step of the method without the λ term, but with an effective ∆t modified as

φ n+1 φ n+1 φ n ∆tHφ n
∆ λ,FE t = ∆t 1 + ∆tλ , (4.11) ∆ λ,BE t = ∆t 1 -∆tλ . (4.12)
This is not true for the CN method, nor it is true when nonlinear terms are included. However, even in this case, the difference between including and not including the λ term is O(∆t 2 ), and their behavior is similar. Since the analysis of unprojected gradient methods is simpler, we focus on this here. Then, these schemes can all be written in the form

φ n+1 = Aφ n Aφ n , (4.13)
where the matrix A is given by

A FE = I -∆tH, A BE = (I + ∆tH) -1 , A CN = (I + ∆t 2 H) -1 (I - ∆t 2 H).
The eigenvalues µ i of A are related to the eigenvalues λ i of H by the following spectral transformation

µ FE i = 1 -∆tλ i , µ BE i = 1 1 + ∆tλ i , µ CI i = 1 -∆t 2 λ i 1 + ∆t 2 λ i .
We call this eigenvalue the amplification factor: if φ n has eigencomponents c n,i = v i , φ n on the eigenvector v i of H and φ 0 is normalized to 1, then the iteration (4.13) can be solved as

c n,i = µ n i N i=1 |µ n i c 0,i | 2 c 0,i .
This iteration converges towards the eigenvector associated to the largest eigenvalue µ N (in modulus) of A, if it is simple, with convergence rate µ N -1 /µ N . This is nothing but the classical power method for the computation of eigenvalues, with a spectral transformation from H to A. Therefore we identify the FE method as a shifted power method, the BE method as a shift-and-invert approach, and the CN uses a generalized Cayley transform [START_REF] Bai | Templates for the solution of algebraic eigenvalue problems: a practical guide[END_REF][START_REF] Saad | Numerical Methods for Large Eigenvalue Problems[END_REF].

From this we can readily see the properties of the different schemes. We make the assumption that either λ 1 is positive or that ∆t < 1 -λ1 (for BE) and ∆t < 2 -λ1 (for CN). If this condition is not verified, then the iteration will generally not converge towards an eigenvector associated with λ 1 because another eigenvalue than λ 1 will have a larger amplification factor. Under this assumption, we see that the BE and CN converge unconditionally, while FE only converges if

∆t < 2 λ N .
This is a CFL-like condition: when H is the discretization of an elliptic operator, λ N will tend to infinity as the size of the basis increases, which will force FE to take smaller time steps.

The asymptotic convergence rate of these methods is µ N -1 µ N . While the FE method has a bounded convergence rate, imposed by λ 1 , λ 2 and λ N , the BE and CN methods can be made to have an arbitrarily small convergence rate, by simply choosing ∆t arbitrarily close to -1 λ1 (BE) or -2 λ1 (CN). Since in practice λ 1 is unknown, it has to be approximated, for instance by λ(φ n ). This yields the classical Rayleigh Quotient Iteration:

φ n+1 = (H -φ n , Hφ n ) -1 φ n (H -φ n , Hφ n ) -1 φ n ,
which is known to converge cubically. This iteration can also be seen as a Newton-like method.

From the previous considerations, it would seem that the BE and CN are infinitely superior to the FE method: even with a fixed stepsize, the BE and CN methods are immune to CFL-like conditions, and with an appropriately chosen stepsize, it can be turned into a superlinearly-convergent scheme. The first difficulty with this approach is that it is a linear strategy, only guaranteed to converge when close to the ground state. As is always the case with Newton-like methods, it requires a globalization strategy to be efficient and robust in the nonlinear setting. The second issue, is that the BE and CN method require the solution of the linear system (Hλ(φ n )) φn+1 = φ n .

The difficulty of solving the system depends on the discretization scheme used. For localized basis schemes like the Finite Element Method, H is sparse, and efficient direct methods for large scale sparse matrices can be used [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. For the Fourier pseudo-spectral scheme, which we use in this paper, H is not sparse, and only matrix-vector products are available efficiently through FFTs (a matrix-free problem). This means that the system has to be solved using an iterative method. Since it is a symmetric but indefinite problem, the solver of choice is MINRES [START_REF] Barrett | Templates for the solution of linear systems: building blocks for iterative methods[END_REF], although the solver BICGSTAB has been used [START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF][START_REF] Antoine | Modeling and Computation of Bose-Einstein Condensates: Stationary States, Nucleation, Dynamics, Stochasticity[END_REF]. The number of iterations of this solver will then grow when the grid spacing tends to zero, which shows that BE also has a CFL-like limitation. However, as is well-known, Krylov methods [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF] only depend on the square root of the condition number for their convergence, as opposed to the condition number itself for fixed-point type methods [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]. This explains why BE with a Krylov solver is preferred to FE in practice [START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF][START_REF] Antoine | Modeling and Computation of Bose-Einstein Condensates: Stationary States, Nucleation, Dynamics, Stochasticity[END_REF].

Furthermore, preconditioners can be used to reduce this number of iterations [START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF][START_REF] Antoine | Modeling and Computation of Bose-Einstein Condensates: Stationary States, Nucleation, Dynamics, Stochasticity[END_REF], e.g. with a simple preconditioner (one that is diagonal either in real or in Fourier space). This method is effective for many problems, but requires a globalization strategy, as well as an appropriate selection of parameters such as ∆t and the precision used to solve the linear system [START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF][START_REF] Antoine | Modeling and Computation of Bose-Einstein Condensates: Stationary States, Nucleation, Dynamics, Stochasticity[END_REF]. Here, we propose a method that has all the advantages of BE (robust, Krylov-like dependence on the square root of the condition number, ability to use a preconditioner), but is explicit, converges faster than BE, and has no free parameter (no fine-tuning is necessary).

The Preconditioned nonlinear Gradient (PG) and Conjugate Gradient (PCG) methods

The gradient method

The previous approaches usually employed in the literature to compute the ground states of the Gross-Pitaevskii equation are all based on implicit discretizations of the imaginary-time equation (4.4). As such, these methods come from PDE theory and lack the insight of minimization algorithms. The difficulty of applying classical minimization algorithms comes from the spherical constraints. However, the general theory of optimization algorithms on Riemannian manifolds has been developed extensively in [START_REF] Absil | Optimization Algorithms on Matrix Manifolds[END_REF][START_REF] Edelman | The geometry of algorithms with orthogonality constraints[END_REF], where the authors derive constrained analogues of gradient, conjugate gradient and Newton algorithms. This is the approach we follow here, and employ a preconditioned conjugate gradient method on the manifold S.

In this section, we work with an arbitrary symmetric positive definite preconditioner P . The choice of P will be discussed later in subsection 5.5. The (projected, preconditioned) gradient method for the minimization of E on S is the update

φn+1 = φ n -α n P H φn φ n -λ n φ n , φ n+1 = φn+1 / φn+1 , (5.14)
where λ n = λ(φ n ). We reformulate this equation as

φ n+1 = cos(θ n )φ n + sin(θ n ) p n p n , with p n = d n -Re d n , φ n φ n , (5.15)
where d n = -P r n is the descent direction, equal to the negative of the preconditioned residual r n = H φn φ n -λ(φ n )φ n . The equations (5.14) and (5.15) are equivalent when θ n or α n is small enough, with a one-to-one correspondance between θ n and α n . To first order, we have:

α n = θ n p n .
Without preconditioner, this method, summarized in Algorithm 1, is identical to the FE method (4.8).

while not converged do

λ n = λ(φ n ) r n = H φn φ n -λ n φ n d n = -P r n p n = d n -Re d n , φ n φ n θ n = arg min θ E cos(θ)φ n + sin(θ)p n / p n φ n+1 = cos(θ n )φ n + sin(θ n )p n / p n n = n + 1 end Algorithm 1:
The gradient method

To choose the parameter θ n , a number of strategies are possible. We first show that, when θ n is small enough, the gradient method decreases the energy.

Expanding φ n+1 up to second-order in θ n , we obtain

φ n+1 = 1 - θ 2 n 2 φ n + θ n p n p n + O(θ 3 n ), (5.16)
and therefore

E(φ n+1 ) = E(φ n ) + θ n p n Re ∇E(φ n ), p n + 1 2 θ 2 n p n 2 ∇ 2 E(φ n )[p n , p n ] -λ n p n 2 + O(θ 3 n ).
(5.17)

We now compute the first-order variation

Re ∇E(φ n ), p n = Re ∇E(φ n ), d n -Re d n , φ n φ n d n -Re d n , φ n φ n = Re r n , d n d n -Re d n , φ n φ n = - r n , P r n d n -Re d n , φ n φ n .
Since P was assumed to be positive definite, this term is always negative so that the algorithm decreases the energy when θ n is chosen small enough. Since p n is orthogonal to φ n , the second-order term ∇ 2 E(φ n )[p n , p n ] -λ n p n 2 is guaranteed to be positive when φ n is close to a minimizer by the second-order optimality conditions. Therefore, a basic strategy is to choose θ n fixed and small enough so that the energy decreases. A better one is to choose θ n adaptively. For instance, we could perform the linesearch

θ n = arg min θ E cos(θ)φ n + sin(θ) p n p n . (5.18)
Since E(θ) is not a quadratic function, this is a nonlinear one-dimensional minimization problem, generally requiring many evaluations of E(θ) to converge to a minimum. However, many of the computations for the evaluation of E(θ), including all that require FFTs, can be pre-computed. Since the FFT step is the dominant one in the computation of the energy, the evaluation of E(θ) at many points is not much more costly than the evaluation at a single point. Therefore it is feasible to use a standard one-dimensional minimization routine. Alternatively, we can obtain a simple and cheap approximation by minimizing the second-order expansion of E in θ n . We expect this to be accurate when θ n is small, which is the case close to a minimizer. Minimizing (5.17) with respect to θ n yields

θ opt n = -Re ∇E(φ n ), p n p n Re ∇ 2 E(φ n )[p n , p n ] -λ n . (5.19)
As we have seen, the numerator is always positive, and the denominator is positive when φ n is close enough to a minimizer. In our implementation, we compute the denominator, and, if it is positive, we use θ opt n as a trial stepsize. If not, we use some default positive value. If the energy of φ n+1 using this trial stepsize is decreased, we accept the step. If the energy is not decreased, we reject the step, decrease the trial stepsize, and try again, until the energy is decreased (which is mathematically ensured when θ n is small enough). Alternatively, we can use Armijo or Wolfe conditions as criterion to accept or reject the stepsize, or even use the full line search (5.18). The evaluation of the energy at multiple values of θ do not require more Fourier transforms than at only one point, but only more computations of the nonlinear term, so that a full line search is not much more costly than the heuristic outlined above. In our tests however, the heuristic above was sufficient to ensure fast convergence, and a full line search only marginally improved the number of iterations. Therefore, we simply use the heuristic in the numerical results of Section 6.

Let us note that under reasonable assumptions on the structure of critical points and on the stepsize choice, there are various results on the convergence of this iteration to a critical point (see [START_REF] Absil | Optimization Algorithms on Matrix Manifolds[END_REF] and references therein).

The conjugate gradient method

The conjugate gradient method is very similar, but uses an update rule of the form

d n = -P r n + β n p n-1 (5.20)
instead of simply d n = -P r n . This is justified when minimizing unconstrained quadratic functionals, where the formula

d n = -P r n + β n d n-1 , with β n = r n , P r n r n-1 , P r n-1 , (5.21)
yields the well-known PCG method to solve linear systems. For nonlinear problems, different update formulas can be used, all equivalent in the linear case. Equation (5.21) is known as the Fletcher-Reeves update. Another popular formula is the Polak-Ribière choice β = max(β PR , 0), where

β PR = r n -r n-1 , P r n r n-1 , P r n-1 . (5.22)
We use β = max(β PR , 0), which is equivalent to restarting the CG method (simply using a gradient step) when β PR < 0 and is a standard choice in nonlinear CG methods. For the justification of the CG method for constrained minimization, see [START_REF] Absil | Optimization Algorithms on Matrix Manifolds[END_REF][START_REF] Edelman | The geometry of algorithms with orthogonality constraints[END_REF]. The CG algorithm is presented in Algorithm 2. In contrast with the gradient algorithm, the quantity Re ∇E(φ n ), p n does not have to be negative, and p n might not be a descent direction: even with a small stepsize, the energy does not have to decrease at each step. To obtain a robust minimization

while not converged do λ n = λ(φ n ) r n = H φn φ n -λ n φ n β n = r n -r n-1 , P r n / r n-1 , P r n-1 β n = max(β n , 0) d n = -P r n + βp n-1 p n = d n -Re d n , φ n φ n θ n = arg min θ E cos(θ)φ n + sin(θ)p n / p n φ n+1 = cos(θ n )φ n + sin(θ n )p n / p n n = n + 1 end Algorithm 2:
The conjugate gradient method method, we enforce energy decrease to guarantee convergence. Therefore, our strategy is to first check if p n is a descent direction by computing Re ∇E(φ n ), p n . If p n is not a descent direction, we revert to a gradient step, which we know will decrease the energy, else, we choose θ n as in (5.19), and use the same stepsize control as in the gradient algorithm.

In our numerical tests, we observe that these precautions of checking the descent direction and using a stepsize control mechanism are useful in the first stage of locating the neighborhood of a minimum. Once a minimum is approximately located, p n is always a descent direction and the stepsize choice (5.19) always decreases the energy.

Stopping criteria

A common way to terminate the iteration (in the BE schemes) is to use the stopping criterion

φ n,∞ err := φ n+1 -φ n ∞ ≤ ε. (5.23)
This can be problematic because the minima are generally not isolated but form a continuum due to symmetries (for instance, complex phase or rotational invariance), and this criterion might be too restrictive. A more robust one is based on the norm of the symmetry-covariant residual

r n,∞ err := r n ∞ = H φn φ n -λ n φ n ∞ ≤ ε. (5.24)
or the symmetry-invariant energy difference

E n err := |E(φ n+1 ) -E(φ n )| ≤ ε. (5.25)
This third one converges more rapidly than the two previous ones: as is standard in optimization, when φ * is a minimum of the constrained minimization problem and φ ∈ S, then

E(φ) -E(φ * ) = O( φ -φ * 2 ).
This is consistent with our results in Figure 7.

In the current paper, we always use the energy based stopping criterion (5.25): for the 2D and 3D cases, a criteria based on φ n,∞ err or r n,∞ err can lead to long computational times, most particularly for large rotations ω, even without changing the energy (see the example in subsection 6.2).

Convergence analysis

A full analysis of the convergence properties of our methods is outside the scope of this paper, but we give in this section some elementary properties, and heuristics to understand their asymptotic performance.

Based on the expansion of the energy (5.17) as a function of θ for the gradient method, it is straightforward to prove that, when E is bounded from below and the stepsize θ is chosen optimally, the norm of the projected gradient H φn φ n -λ n φ n converges to 0. Convergence guarantees for the conjugate gradient method are harder, but can still be proven under a suitable restart strategy that ensures that the energy always decreases fast enough (for instance, the Armijo rule).

With additional assumptions on the non-degeneracy of critical points, we can even prove the convergence of φ n to a critical point, that will generically be a local minimum. However, the question of the precise convergence speed of the gradient and conjugate gradient algorithms we use is problematic, because of three factors: the constraint φ = 1, the non-quadraticity of E, and the presence of a preconditioner. To our knowledge, no asymptotically optimal bound for this problem has been derived. Nevertheless, based on known results about the convergence properties of the conjugate gradient method for preconditioned linear systems on the one hand [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF], and of gradient methods for nonlinear constrained minimization [START_REF] Absil | Optimization Algorithms on Matrix Manifolds[END_REF] on the other, it seems reasonable to expect that the convergence will be influenced by the properties of the operator

M = (1 -φφ * )P (∇ 2 E(φ) -λI)(1 -φφ * ), (5.26)
where φ is the minimum and λ = λ(φ). This operator admits 0 as its lowest eigenvalue, associated with the eigenvector φ. It is reasonable to expect that the convergence rate will be determined by a condition number σ equals to the ratio of the largest to the lowest non-zero eigenvalue of this operator. As is standard for linear systems, we also expect that the number of iterations to achieve a given tolerance will behave like √ σ for the conjugate gradient algorithm, and σ for the gradient algorithm. As we will see in Section 6, this is verified in our tests.

The Hessian operator ∇ 2 E(φ), which includes a Laplace operator, is not bounded. Correspondingly, on a given discretization domain, when the grid is refined, the largest eigenvalues of this operator will tend to +∞. For a linear meshsize h, the eigenvalues of ∇ 2 E(φ) will behave as O(h -2 ). This is another instance of the CFL condition already seen in the discretization of the imaginary time equation. The Hessian ∇ 2 E also includes a potential term V , which is often confining and therefore not bounded, such as the classical harmonic potential V (x) = |x| 2 , or more generally confining potentials whose growth at infinity is like |x| p . for some p > 0 Thus, even with a fixed meshsize h on a domain [-L, L] d , when L is increased, so will the largest eigenvalues of ∇ 2 E(φ), with a O(L p ) growth. When a (conjugate) gradient method is used without preconditioning, the convergence will be dominated by modes associated with largest eigenvalues of ∇ 2 E. This appears in simulations as high-frequency oscillations and localization at the boundary of the domain of the residual r n .

To remedy these problems and achieve a good convergence rate, adequate preconditioning is crucial.

Preconditioners

We consider the question of building preconditioners P for the algorithms presented above. In the schemes based on the discretization of the gradient flow, preconditioning is naturally needed when solving linear systems by iterative methods. In the gradient and conjugate gradient optimization schemes, it appears as a modification of the descent direction to make it point closer to the minimum:

d n := -P (H φn φ n -λ n φ n ).
(5.27)

In both cases, the preconditioning matrix should be an approximation of the inverse of the Hessian matrix of the problem.

Kinetic energy preconditioner.

One of these approximations is to use only the kinetic energy term

P ∆ = (α ∆ -∆/2) -1 , (5.28)
where α ∆ is a positive shifting constant to get an invertible operator, and I is the identity operator. This has been called a "Sobolev gradient" in [START_REF] Danaila | A new Sobolev gradient method for direct minimization of the Gross-Pitaevskii energy with rotation[END_REF] because it is equivalent to taking the gradient of the energy in the Sobolev H 1 -norm (with α ∆ = 1/2). In the framework of the BESP scheme for the GPE with Krylov solver, a similar preconditioner has been proposed in [START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF], α ∆ being the inverse of the time step ∆t of the semi-implicit Euler scheme. A closely-related variant is standard in planewave electronic structure computation [START_REF] Zhou | On the preconditioning function used in planewave dft calculations and its generalization[END_REF], where it is known as the Tetter-Payne-Allan preconditioner [START_REF] Teter | Solution of Schrödinger's equation for large systems[END_REF]. This preconditioner is diagonal in Fourier space and can therefore be applied efficiently in our pseudo-spectral approximation scheme. On a fixed domain [-L, L] d , the effect of this preconditioner is to make the number of iterations independent from the spatial resolution h, because P ∇ 2 E(φ), seen as an operator on the space of functions on [-L, L] d , will be equal to the identity plus a compact operator. This is supported by numerical experiments in Section 6. However, this operator is not bounded in the full domain R d . Therefore, as L increases, so will the largest eigenvalues of P ∇ 2 E(φ). For a potential V that grows at infinity like |x| p , the largest eigenvalues of P ∇ 2 E(φ) are O(L p ), resulting in an inefficient preconditioner. Similarly, when η is large, the nonlinear term becomes dominant, and the kinetic energy preconditioner is inefficient. The choice of α ∆ is a compromise: if α ∆ is too small, then the preconditioner will become close to indefinite, which can produce too small eigenvalues in the matrix (5.26) and hamper convergence. If α ∆ is too big, then the preconditioner does not act until very large frequencies, and large eigenvalues result. We found that a suitable adaptive choice, that has consistently good performance and avoids free parameters, is

α ∆ = λn := 1 2 |∇φ n | 2 + V |φ n | 2 + η|φ n | 4 dx > 0 (5.29)
which is a positive number that represents the characteristic energy of φ n . We use this choice for our numerical simulations.

Potential energy preconditioner. Another natural approach is to use the potential energy term for the preconditioner:

P V = (α V + V + η|φ n | 2 ) -1 . (5. 30
)
This preconditioner is diagonal in real space and can therefore be applied efficiently. Dual to the previous case, this preconditioner has a stable performance when the domain and η are increased, but deteriorates as the spatial resolution is increased. Such a preconditioner has been used in [START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF] when the gradient flow for the GPE is discretized through a BE scheme, leading then to a Thomas-Fermi preconditioner. In this study, the parameter α V was 1/∆t. As in the kinetic energy case, we found it efficient to use α V = λn , and we will only report convergence results for this choice of parameter.

Combined preconditioner. In an attempt to achieve a stable performance independent of the size of the domain or the spatial resolution, we can define the combined preconditioners

P C1 = P V P ∆ , P C2 = P ∆ P V (5.31)
or a symmetrized version

P C = P 1/2 V P ∆ P 1/2 V .
(5.32)

With these preconditioners, P ∇ 2 E(φ) is bounded as an operator on L 2 (R d ) (this can be proven by writing explicitly its kernels in Fourier space and then using Schur's test). However, we found numerically that this operator is not bounded away from zero, and has small eigenvalues of size O(L -p + h 2 ). Therefore, the conditioning deteriorates as both the spatial resolution and the size of the domain increase. In summary, for a spatial resolution h and a domain size L, the asymptotic condition numbers of the preconditioned Hessian with these preconditioners are

κ ∆ = O(L p ), κ V = O(h -2 ), κ C = O 1 L -p +h 2 = O(min(L p , h -2 )).
(5.33) Therefore, the combined preconditioners act asymptotically as the best of both the kinetic and potential preconditioners. However, they might not be more efficient in the pre-asymptotic regime and require additional Fourier transforms.

Computational efficiency. The application the operator P V is almost free (since it only requires a scaling of φ), but the naive application of P ∆ requires a FFT/IFFT pair. However, since we apply the preconditioners after and before an application of the Hamiltonian, we can reuse FFT and IFFT computations, so that the application of P ∆ does not require any additional Fourier transform. Similarly, the use of P C1 and P C2 only require one additional Fourier transform per iteration, and that of the symmetrized version P C two.

In summary, the cost in terms of Fourier transforms per iteration for the rotating GPE model is

• no preconditioner: 3 FFTs/iteration (get the Fourier transform of φ, and two IFFTs to compute ∆φ and L z φ respectively),

• P ∆ or P V : 3 FFTs/iteration,

• non-symmetric combined P C1 or P C2 : 4 FFTs/iteration,

• symmetric combined P C : 5 FFTs/iteration.

Note that this total might be different for another type of GPE model e.g. when a nonlocal dipole-dipole interaction is included [START_REF] Antoine | Modeling and Computation of Bose-Einstein Condensates: Stationary States, Nucleation, Dynamics, Stochasticity[END_REF][START_REF] Bao | Computing the ground state and dynamics of the nonlinear Schrödinger equation with nonlocal interactions via the nonuniform FFT[END_REF]. As we will see in Section 6, all combined preconditioners have very similar performance, but the symmetrized one might be more stable in some circumstances. A theoretical explanation of these observations, and in particular of the effect of a non-symmetric preconditioner is, to the best of our knowledge, still missing.

Numerical results

We first introduce some notations. When we combine one of the preconditioners P ν (ν = I, ∆, V , C, C 1 , C 2 ) (5.28)-(5.32) with the gradient algorithm (Alg: 1), we denote the resulting methods by PG ν . Similarly, we denote by PCG ν if the preconditioned conjugate gradient algorithm (Alg: 2) was applied. In the following, we denote by #iter the number of global iterations for an iterative algorithm to get the converged solution with an a priori tolerance ε with respect to the stopping criterion (5.25).

Concerning the BESP schemes (4.6) and (4.9), at each outer iteration n, one needs to solve an implicit system with the operator (1/∆t + H φn ). We use a Krylov subspace iterative solver (MINRES here) with one of the preconditioners P ν (ν = I, ∆, V, C) (5.28)- (5.31) to accelerate the number of inner iterations [START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF]. The preconditioned BESP schemes is then denoted by BE ν , according to the chosen preconditioner. The number of iterations reported is equal to the sum of the inner iterations over the outer iterations.

In the following numerical experiments, we consider two types of trapping potential V (x): the harmonic plus lattice potential [START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF]

V (x) = V 0 d (x) +        κ x sin 2 (q x x 2 ), ν=x,y κ ν sin 2 (q ν ν 2 ), ν=x,y,z κ ν sin 2 (q ν ν 2 ), with V 0 d (x) =        γ 2 x x 2 , d = 1, ν=x,y γ ν ν 2 , d = 2, ν=x,y,z γ ν ν 2 , d = 3, (6.34)
and the harmonic plus quartic potential for d = 2, 3 [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 of the Gross-Pitaevskii energy with rotation[END_REF][START_REF] Zeng | Efficiently computing vortex lattices in rapid rotating Bose-Einstein condensates[END_REF]]

V (x) = (1 -α)V 0 2 (x) + κ (x 2 + y 2 ) 2 4 + 0, d = 2, γ 2 z z 2 , d = 3. (6.35)
Moreover, unless stated otherwise, we take the initial data as the Thomas Fermi approximation [START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF][START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF]:

φ 0 = φ TF g φ TF g , with φ TF g (x) =    µ TF g -V (x) /η, V (x) < µ TF g , 0, otherwise, (6.36)
where

µ TF g = 1 2      (3ηγ x ) 2/3 , d = 1, (4ηγ x γ y ) 1/2 , d = 2, (15ηγ x γ y γ z ) 2/5 , d = 3. (6.37)
The algorithms were implemented in Matlab (Release 8.5.0).

Numerical results in 1D

Here, V (x) is chosen as the harmonic plus lattice potential (6.34) with γ x = 1, k x = 25 and q x = π 2 . The computational domain and mesh size are respectively denoted as D = [-L, L] and h. In addition, to compare with the common existing method BESP, we choose the stopping criteria (5.25) with ε = 10 -14 all through this section. For BESP, we choose ∆t = 0.01 unless specified otherwise, and fix the error tolerance for the inner loop to 10 -10 . Other values of the error tolerance were also tried, but this choice was found to be representative of the performance of BESP. Example 6.1. We first compare the performance of various solvers without preconditioning in a simple case. We choose L = 16 and η = 250, and varying mesh sizes. We compare the method BE I given by (4.9), BE λ I given (4.6), and the gradient and conjugate gradient algorithms in Figure 2. The difference between methods BE I and BE λ I is the inclusion of the chemical potential in the discretized gradient flow: we showed in Section 4 that both were equivalent for linear problems up to a renormalization in ∆t. We see here that this conclusion approximately holds even in the nonlinear regime (η = 0), with both methods performing very similarly until ∆t becomes large, at which point the BE λ I effectively uses a constant stepsize (see (4.12)), while the large timestep in BE I makes the method inefficient. In this case, λ 1 is positive, so that both methods converge to the ground state even for a very large ∆t. Overall we see that the optimum number of iterations is achieved for a value of ∆t of about 0.01, which we keep in the following tests to ensure a fair comparison. We also use the BE λ variant in the following tests.

For modest values of the discretization parameter h, the Backward Euler methods are less efficient than the gradient method (which can be interpreted as a Forward Euler iteration with adaptive stepsize). As h is decreased, the conditioning of the problem increases as h -2 . The gradient/Forward Euler method is limited by its CFL condition, and its number of iterations grows like h -2 , as can readily be checked in Figure 2. The Backward Euler methods, however, use an efficient Krylov solver that is only sensitive to the square root of the conditioning, and its number of iterations grows only like h -1 . Therefore it become more efficient than the gradient/Forward Euler method.

The conjugate gradient method is always more efficient than the other methods by factors varying between one and two orders of magnitude. Its efficiency can be attributed to the combination of Krylovlike properties (as the Backward Euler method, its iteration count displays only a h -1 growth) and optimal stepsizes. Example 6.2. We compare now the performance of the (conjugate) gradient method with different preconditioners. To this end, we consider the algorithms PG ν and PCG ν with ν = ∆, V, C, C 1 , C 2 . The computational parameters are chosen as L = 128 and h = 1 64 , respectively. Fig. 3 shows the iteration number #iter for these schemes and different values of the nonlinearity strength η. From this figure and other numerical results not shown here, we can see that: (i) For each fixed preconditioner, the PCG schemes works better than the PG schemes; (ii) the combined preconditioners all work equally well, and bring a reduction in the number of iteration, at the price of more Fourier transforms. Example 6.3. In this example, we compare the performance of PG ν , PCG ν and BE ν (ν = I, ∆, V, C) with respect to different domain and mesh sizes. To this end, we fix η = 250. Fig. 4 shows the total 0 2.5 5 7. 5 10 η/10000 iteration number for these schemes with different L and h. From this figure and additional numerical results not shown here for brevity, we see that: (i) Preconditioned solvers outperform unpreconditioned solvers; (ii) The potential preconditioner P V (5.30) makes the solver mainly depend on the spatial resolution h, while the kinetic potential preconditioner P ∆ (5.28) prevents the deterioration as h decreases for a fixed L, consistent with the theoretical analysis in subsection 5.5; (iii) The deterioration is less marked for Krylov-based methods (BE and PCG) than for the PG method, because Krylov methods only depends on the square root of the condition number (iv) The combined preconditioner P C (5.31) makes the solvers almost independent of both the parameters h and L, although we theoretically expect a stronger dependence. We attribute this to the fact that we start with a specific initial guess that does not excite the slowly convergent modes enough to see the dependence on h and L; (v) For each solver, the combined preconditioner P C performs best. Usually, PCG C is the most efficient algorithm, followed by PG C , and finally BE C .

Example 6.4. We investigate further the performance of PCG C and PG C with different nonlinear interaction strenghts η. To this end, we take L = 128 and different discretization parameters h. We vary the nonlinearity from η = 0 to η = 10 5 . Fig. 5 depicts the corresponding iteration numbers to converge. We could clearly see that: (i) The iteration counts for both methods are almost independent on h, but both depend on the nonlinearity η; PCG C depends slightly on η while PG C is more sensitive; (ii) For fixed η and h, PCG C converges much faster than PG C .

From examples 6.2-6.4, we see that PCG C , i.e. PCG with combined symmetric preconditioner P C (5.31) is the best solver. Hereafter, unless stated, we use PCG C as the default solver to compute the ground states.

Numerical results in 2D

Here, we choose V (x) as the harmonic plus quartic potential (6.35) with γ x = γ y = 1, α = 1.2 and κ = 0.3. The computational domain and mesh sizes are chosen respectively as D = [-16, 16] 2 and h = 1 16 . First, we test the evolution of the three errors (5.23)-(5.25) as the algorithm progresses. To this end, we take η = 1000 and ω = 3.5 as example. Fig. 6 plots the φ n,∞ err , r n,∞ err and E n err errors with respect to the iteration number. We can see clearly that E n err converges faster the other two indicators, as expected. Considering φ n,∞ err or r n,∞ err with an improper but relative large tolerance would require a very long computational time to converge even if the energy would not change so much. This is most particularly true for large values of ω. In all the examples below, unless stated, we fix E n err (5.25) with ε = 10 -12 to terminate the code. Example 6.5. In this example, we compare the performance of PCG C and PG C for the 2D rotating case. To this end, V (x) is chosen as the harmonic plus lattice potential (6.34) with γ x = γ y = 1, k x = k y = 25 and q x = q y = π 2 . The computational domain and mesh sizes are chosen respectively as D = [-32, 32] 2 and h = fixed η. Meanwhile, the dependency on η becomes stronger as ω increases. Let us remark here that it would be extremely interesting to build a robust preconditioner including the rotational effects to get a weaker ω-dependence in terms of convergence. Example 6.7. In this example, we apply PCG C to solve some more difficult problems. We compute the ground states φ g of rotating BECs with large values of η and ω. To this end, we take L = 20, h = 1/16 and set the stopping tolerance in (5.25) to ε = 10 -14 . Table 1 lists the CPU times for the PCG C solver to converge while Fig. 9 shows the contour plot of the density function |φ g (x)| 2 for different ω and η.

We can see that the PCG C method converges very fast to the stationary states. Let us remark that, to the best of our knowledge, only a few results were reported for such fast rotating BECs with highly nonlinear (very large η) problems, although they are actually more relevant for real physical problems. Hence, PCG C can tackle efficiently difficult realistic problems on a laptop. Example 6.8. The choice of the initial data also affects the final converged stationary states. Since all the algorithms we discussed are local minimisation algorithms, inappropriate initial guess might lead to local minimum. To illustrate this claim, we take ε = 10 -14 , V (x) = |x| 2 2 , η = 500 and compute the ground states of the rotating GPE for different ω and 10 types of frequently used initial data

(a) φ a (x) = 1 √ π e -(x 2 +y 2 )/2 , (b) φ b (x) = (x + iy)φ a (x), ( b) φb(x) = φb (x), (6.38) (c) φ c = (φ a (x) + φ b (x))/2 (φ a (x) + φ b (x))/2 , (d) φ d = (1 -ω)φ a (x) + ωφ b (x) (1 -ω)φ a (x) + ωφ b (x) , (6.39)
(e) φ e = ωφ a (x) + (1 -ω)φ b (x) ωφ a (x) + (1 -ω)φ b (x) , (c) φ c(x) = φc (x), ( d) φ d(x) = φd (x), (6.40)
(ē) φ ē(x) = φe (x), (f ) Thomas Fermi approximation (6.36). (6.41)

Another approach to prepare some initial data is as follows: we first consider one of the above initial guess (a)-(f), next compute the ground state on a coarse spatial grid, say with a number of grid points M p × M p (with M p = 2 p ), and then denote the corresponding stationary state by φ p g . We next subdivide the grid with M p+1 × M p+1 points, interpolate φ p g on the refined grid M p+1 × M p+1 to get a new initial data and launch the algorithm at level p + 1, and so on until the finest grid with M × M points where the converged solution is still denoted by φ g . Similarly to [START_REF] Wu | A regularized Newton method for computing ground states of Bose-Einstein condensates[END_REF], this multigrid technique is applied here with the coarsest grid based with M 6 = 2 6 and ends with the finest grid M = 2 9 . We use the tolerance parameters ε = 10 -14 for M = 2 9 , and ε = 10 -12 for p = 6, 7, 8.

Tables 2 and3 list the energies obtained by PCG C via the fixed and multigrid approaches, respectively, for different initial data and ω. The stationary states φ g (x) with lowest energies are marked by underlines and the corresponding CPU times are listed in the same Table . Moreover, Fig. 10 shows the contour plots |φ g (x)| 2 of the converged solution with lowest energy obtained by the multigrid approach. Now, let us denote by E p n := E(φ p n ) the evaluated energy at step n for a discretization level p = 6, 7, 8, and let E g = E(φ g ) the energy for the converged stationary state for the finest grid. Then, we represent on Fig. 11 the evolution of log 10 (|E p n -E g |) vs. the CPU time for a rotating velocity ω = 0.95. For comparison, we also show the corresponding evolution obtained by the fixed grid approach. The contour plots of |φ p g (x)| 2 obtained for each intermediate coarse grid for p = 6, 7, 8, and the initial guess are also reported. From these Tables and Figures, we can see that: (i) Usually, the PCG C algorithm with an initial data of type (d) or ( d) converges to the stationary state of lowest energy; (ii) The multigrid approach is more robust than the fixed grid approach in terms of CPU time and possibility to obtain a stationary state with lower energy. Here, we apply the PCG C algorithm to compute some realistic 3D challenging problems. To this end, V (x) is chosen as the harmonic plus quartic potential (6.35), with γ x = γ y = 1, γ z = 3, α = 1.4 and κ = 0.3. The computational domain is D = [-8, 8] 3 and the mesh size is: h = 1 8 . We test four cases: (i) η = 100, ω = 1.4; (ii) η = 100, ω = 1.8; (iii) η = 5000 and ω = 3; (iv) η = 10000 and ω = 3. The initial guess is always taken as the Thomas-Fermi initial data and the multigrid algorithm is used. Fig. 12 shows the isosurface |φ g (x)| 2 = 10 -3 and the surface plot of |φ g (x, y, z = 0)| 2 for the four cases. The CPU times for these four cases are respectively 2256 s, 1403 s, 11694 s and 21971 s.

Conclusion

We have introduced a new preconditioned nonlinear conjugate gradient algorithm to compute the stationary states of the GPE with fast rotation and large nonlinearities that arise in the modeling of Bose-Einstein Condensates. The method, which is simple to implement, appears robust and accurate. In addition, it is far more efficient than standard approaches as shown through numerical examples in 1D, 2D and 3D. Furthermore, a simple multigrid approach can still accelerates the performance of the method and leads to a gain of robustness thanks to the initial data. The extension to much more general systems of GPEs is direct and offers an interesting tool for solving highly nonlinear 3D GPEs, even for very large rotations.

Figure 1 :

 1 Figure 1: Update rule with the projected gradient, with the λ term (red) and with the unprojected gradient, without the λ term (blue).

Figure 2 :

 2 Figure 2: Example 6.1. Number of iterations to converge for different methods and different stepsizes, with different values of the discretization parameter h.

1 8 .Figure 4 :

 84 Figure 4: Example 6.3. Number of iterations to converge for BEν , PGν and PCGν for ν = I, V, ∆, C, vs. the mesh refinement h.

Figure 5 :Figure 6 :

 56 Figure 5: Example 6.4. Number of iterations to converge for PCG C (left) and PG C (right) with L = 128, and various values of h and η.

2 Figure 7 :Example 6 . 6 .Figure 8 :

 27668 Figure 7: Exemple 6.5. Number of iterations for PCG C and PG C for ω = 0 (left) and PCG C for different ω (right) vs. η.

Figure 10 :

 10 Figure 10: Contour plots of |φg(x)| 2 corresponding to the lowest energy levels in Table3.

Figure 11 :

 11 Figure 11: Exemple 6.8. Energy error log 10 (|E p n -Eg|) vs. the accumulated CPU time for ω = 0.95 with initial data (d) in Table 2 (p = 9, upper left) and respectively Table 3 (upper right) as well as the stationary state obtained at each intermediate level (lower, p = 6, 7, 8)

Figure 12 :

 12 Figure 12: Exemple 6.9. Isosurface |φg(x)| 2 = 10 -3 (upper) and surface plot of |φg(x, y, z = 0)| 2 (lower) in example 6.9. The CPU cost for these four cases are respectively 2256 (s), 1403 (s), 11694 (s) and 21971 (s).

  Example 6.2. Number of iterations of PGν and PCGν (ν = ∆, V, C, C 1 , C 2 ) to converge, for different nonlinear strengths η.

		10 4			500					PCG ∆ PCG V
										PCG C
	#iter	10 2 10 3	PG C PG V PG ∆	#iter	200 350					PCG C1 PCG C2
		10 1	PG C1 PG C2		50	0	2	4	6	8	10
								η/10000
	Figure 3:

Table 1 :

 1 CPUs time (seconds) for PCG C to compute the ground states of the GPE with various ω and η.

	η	ω=1	1.5	2	2.5	3	3.5	4	4.5
	1000	493	551	560	2892	2337	720	966	3249
	5000	1006	1706	867	6023	1144	1526	12514	19248
	10000	4347	21525	5511	15913	15909	6340	16804	32583

Figure 9: Exemple 6.7. Corresponding contour plots of the density function |φg(x)|

2 of Table 1.

Table 2 :

 2 Exemple 6.8. Fixed grid approach (with M = 29 ): converged energies and the CPU times (seconds) for the solution with lowest energy (which is underlined). .5118 8.2606 9.2606 8.0246 8.0197 8.0246 8.0197 8.0246 8.0197 8.0246 176.0 0.6 8.5118 8.1606 9.3606 7.5845 7.5910 7.5845 7.5845 7.5845 7.5910 7.5845 310.7 0.7 8.5118 8.0606 9.4606 6.9754 6.9731 6.9792 6.9754 6.9754 6.9792 6.9767 542.4 0.8 8.5118 7.9606 9.5606 6.1016 6.0997 6.1031 6.1031 6.1040 6.1019 6.1016 417.0 0.9 8.5118 7.8606 9.6606 4.7777 4.7777 4.7777 4.7777 4.7777 4.7777 4.7777 1051.1 0.95 8.5118 7.8106 9.7106 3.7414 3.7414 3.7414 3.7414 3.7414 3.7414 3.7414 3280.5

	ω	(a)	(b)	(b2)	(c)	(c2)	(d)	(d2)	(e)	(e2)	(f)	CPU
	0.5 8

Table 3 :

 3 Exemple 6.8. Multigrid approach (starting from the coarsest level p = 6 to the finest level p = 9): converged energies and the CPU times (seconds) for the solution with lowest energy (underline).

  CPU 0.5 8.0246 8.0197 8.0197 8.0197 8.0197 8.0197 8.0197 8.0197 8.0197 8.0257 29.5 0.6 7.5845 7.5845 7.5910 7.5845 7.5910 7.5890 7.5845 7.5845 7.5910 7.5845 32.3 0.7 6.9767 6.9726 6.9792 6.9754 6.9731 6.9731 6.9731 6.9757 6.9731 6.9731 53.3 0.8 6.1019 6.1031 6.1019 6.0997 6.1016 6.1016 6.1016 6.1019 6.1016 6.0997 75.2 0.9 4.7777 4.7777 4.7777 4.7777 4.7777 4.7777 4.7777 4.7777 4.7777 4.7777 238.1 0.95 3.7414 3.7414 3.7414 3.7414 3.7414 3.7414 3.7414 3.7414 3.7414 3.7414 621.9 6.3. Numerical results in 3D Example 6.9.

Acknowledgements

X. Antoine and Q. Tang thank the support of the French ANR grant ANR-12-MONU-0007-02 BE-CASIM ("Modèles Numériques" call).


Somayeh Sobati Moghadam
email: somayeh.sobati-moghadam@univ-lyon2.fr

Gérald Gavin
email: gerald.gavin@univ-lyon1.fr

Jérôme Darmont
email: jerome.darmont@univ-lyon2.fr

A Secure Order-Preserving Indexing Scheme for Outsourced Data

Keywords: Cloud computing, Data privacy, Orderpreserving indexing scheme

Cloud computing offers the opportunity of data outsourcing as well as data management. However, because of various privacy issues, confidential data must be encrypted before being outsourced to the cloud. But query processing over encrypted data without decrypting data is a very challenging task. Property-preserving encryption schemes allow encrypting data while still enabling efficient querying over encrypted data. The inherent merits of property-preserving encryption schemes make them very suitable and efficient for cloud data outsourcing. However, the security of such schemes is still a challenge because they are vulnerable to statistical attacks. We present a new order-preserving scheme for indexing encrypted data, as an alternative to propertypreserving schemes, which hides data frequency to achieve a strictly stronger notion of security. The proposed indexing method is secure against statistical attacks. Hence, data cannot be recovered from indexes. Moreover, our scheme is still efficient for query processing.

I. INTRODUCTION

Nowadays, data outsourcing become casual with the advent of cloud computing. Cloud computing appeals to many organizations because of a wide variety of benefits such as cost savings, higher availability, scalability and effective disaster recovery rather than in-house operations [START_REF] Damiani | Key management for multi-user encrypted databases[END_REF]. One of the most notable cloud outsourcing services is database outsourcing (Database-as-a-Service or DBaaS), where individuals and organizations outsource data storage and management to a Cloud Service Provider (CSP) [START_REF] Xiong | Preserving data privacy in outsourcing data aggregation services[END_REF]. Such services allow storing data on a remote CSP and querying data on demand [START_REF] Amanatidis | Provablysecure schemes for basic query support in outsourced databases[END_REF].

Although cloud data outsourcing induces many benefits, it also brings out security and privacy concerns. Privacy issues arise when sensitive data are stored, maintained and processed by an external third party (honest but curious CSP).

A straightforward solution to preserve privacy is encrypting data before outsourcing to the cloud. Encryption protects the exposure of sensitive data even if the server is compromised and ensures that an adversary will be unable to interpret data.

However, when data are encrypted, query processing is not trivial. To overcome this problem, several solutions have been proposed [START_REF] Bellare | Deterministic and efficiently searchable encryption[END_REF], [START_REF] Popa | Cryptdb: protecting confidentiality with encrypted query processing[END_REF], [START_REF] Tu | Processing analytical queries over encrypted data[END_REF]. These approaches implement cryptographic techniques, e.g., order-preserving encryption or homomorphic encryption, which allow computations to be carried out on encrypted data. In the context of relational databases, state-of-the-art solutions use property-preserving encryption (PPE) schemes. PPE schemes enable processing query over encrypted data without decryption. For instance, order-preserving encryption (OPE) preserves the order of plaintext in ciphertexts. Deterministic (DET) schemes encrypt the same plaintext into identical ciphertexts, thus the equality property is preserved [START_REF] Bellare | Deterministic and efficiently searchable encryption[END_REF]. PPE schemes are undoubtedly efficient, enabling queries to be directly processed over encrypted data [START_REF] Liu | New order preserving encryption model for outsourced databases in cloud environments[END_REF].

The term efficient means here that all computation can be processed in time logarithmic in the size of the database, in contrast to performing linear work on each query, which is prohibitively slow for large databases [START_REF] Boldyreva | Orderpreserving symmetric encryption[END_REF]. Any practical PPE scheme inherently leaks some information about underlying data and are vulnerable to statistical attacks [START_REF] Naveed | Inference attacks on property-preserving encrypted databases[END_REF]. In statistical attacks, an adversary possesses some prior knowledge about the plaintext domain and its frequency distribution to gain access to encrypted data. Statistical attacks do not impose any threat when the underlying data has a uniform distribution [START_REF] Sanamrad | Randomly partitioned encryption for cloud databases[END_REF]. CryptDB [START_REF] Popa | Cryptdb: protecting confidentiality with encrypted query processing[END_REF] is the first practical system that uses PPE schemes to support a wide range of queries to be processed over encrypted data. To the best of our knowledge, other systems such as BigQuery demo [START_REF]Google encrytped big query[END_REF], Always Encrypted [START_REF]Always encrypted[END_REF], CipherBase [START_REF] Arasu | Transaction processing on confidential data using cipherbase[END_REF] and Relational Cloud [START_REF] Curino | Relational cloud: a database service for the cloud[END_REF] use PPE schemes, too. As a result, they are vulnerable to statistical attacks. Naveed et al. demonstrate that a large fraction of the records from PPE encrypted columns can be decrypted by statistical attacks [START_REF] Naveed | Inference attacks on property-preserving encrypted databases[END_REF].

Thus, we aim at proposing a secure and feasible scheme that is practical on cloud outsourced databases. In this paper, we focus on the problem of performing range and exact match queries over encrypted databases. We propose an Order-Preserving Indexing (OPI) scheme to address the vulnerability of PPE schemes. Our scheme bears good performance and leads to minimal change for existing software. The proposed indexing scheme is robust against statistical attacks through uniforming the distribution of underlying data. We extend the proposed scheme to introduce a new wrap-around OPI scheme. waOPI provides higher level of security while hiding the 978-1-5090-1072-1/16/$31.00 c 2016 European Union original order of plaintexts.

In our scheme, the user does not need to model data distribution. Our scheme can be used along any encryption scheme to deal with range and exact match queries over encrypted data. Our scheme prevents the CSP from performing statistical analysis, even with prior knowledge about data.

Note that our scheme does not hide information such as the number of data items. One solution for applications that demand the privacy of data item numbers is injecting dummy data items into small data sets to build equally sized data sets [START_REF] Liu | New order preserving encryption model for outsourced databases in cloud environments[END_REF].

The remainder of this paper is organized as follows. Section II gives some preliminaries. The overview of our system is described in section III. Section IV details our indexing scheme. We enhance the security of our scheme in section V, which is followed in section VI by a discussion. Section VII reviews related works. Finally, section VIII provides final conclusions and directions for future work.

II. PRELIMINARIES

A. Property Preserving Encryption

In the context of DBaaS, state-of-the-art solutions rely on PPE schemes. A PPE scheme leaks certain properties of plaintexts by default, which enables an untrusted CSP to compute over encrypted data. Order-preserving encryption (OPE) and deterministic encryption (DET) are two PPE schemes that are used to handle range and exact match queries over encrypted data.

DET allows equality checking by encrypting a plaintext into the same ciphertexts when using the same key k. Thus:

∀x, y : x = y ⇔ EncDET k (x) = EncDET k (y).
DET allows performing SELECT with equality predicates, equality JOIN, GROUP BY, COUNT, DISTINCT queries [START_REF] Popa | Building practical systems that compute on encrypted data[END_REF]. Blowfish and AES+ECB are examples of widely used DET encryption schemes. DET is secure when there is no data redundancy. In contrast to random encryption schemes (a random scheme encrypts the same plaintext into different ciphertexts using the same key), DET is not robust against statistical attacks.

OPE is a deterministic encryption scheme that preserves the order of plaintexts in ciphertexts, i.e., for any key k,

if x ≤ y, then EncOP E k (x) ≤ EncOP E k (y).
The CSP can perform range queries when given encrypted constants EncOP E k (c 1 ) and EncOP E k (c 2 ) corresponding to range [c 1 , c 2 ]. Aggregation queries MIN, MAX, ORDER BY and SORT can also be computed directly over encrypted data.

OPE is a weaker encryption scheme than DET because it reveals order. This weak form of encryption may provide sufficient security for some applications, e.g., when the adversary does not possess any prior knowledge, while increasing query processing efficiently [START_REF] Özsoyoglu | Anti-tamper databases: Querying encrypted databases[END_REF].

B. Statistical Attacks against PPE Schemes

The most common adversarial model in existing solutions is ciphertext-only. It is assumed that the adversary has access only to ciphertexts without other background information. In this context, PPE schemes are vulnerable to statistical attacks [START_REF] Liu | New order preserving encryption model for outsourced databases in cloud environments[END_REF]. In a statistical attack, the adversary has prior knowledge, e.g., the adversary possesses some statistical information about underlying plaintexts. The adversary can use such knowledge to infer about data or even launch an attack. For example, the adversary can map the distribution of ciphertext and plaintext to find all ciphertexts with high frequency.

The major problem with PPE schemes is that they are deterministic, i.e., plaintext and ciphertext have the same distribution. To address this drawback, a one-to-many mapping can be used to flatten the distribution of plaintexts. In one-to-many mappings, one plaintext is mapped into many ciphertexts belonging to a fixed size interval. As a result, highly frequent plaintexts induce dense intervals (figure 1). Hence, the distribution of plaintexts can easily be estimated by differential attacks [START_REF] Li | Security analysis on one-to-many order preserving encryption-based cloud data search[END_REF].

C. Secure Indexes over Encrypted Attributes

A common technique to speed up the execution of queries in databases is to use a pre-computed index. The purpose of secure indexes is to retrieve requested data without decryption. An index is built for a specific attribute that needs to be accessed to evaluate a query. In fact, the biggest difference between a plain index and a secure index is that a secure index must not expose anything about underlying plaintexts.

A secure index for ordered data should support exact match, range and aggregation queries MAX and MIN. A secure index is built on plaintexts before encryption. Then, the index is stored in an auxiliary attribute along encrypted values at the CSP's.

III. OVERVIEW OF SYSTEM

In this section, we briefly discuss the system model where our scheme is applied, our adversary model and our system goals.

A. Basic Model

In our setting, there are three entities.

• The Data Owner (DO) has a large amount of data to be outsourced in the cloud to eliminate database maintenance. The DO accesses its data remotely without revealing the content of its database. The DO has limited computation and storage resources.

• The CSP provides data storage services and computational resources dynamically. • Data Users (DU) are authorized users. DUs can retrieve encrypted data and decrypt them. Our scheme is implemented in two software layers: a user layer that resides at the DO's and a server layer that resides at the CSP's. When a query is issued, the user layer translates it into a query that can be executed directly over the encrypted database. When a query is received by the server, it is executed by the server layer. The computed result is shipped back to the client layer, which decrypts it (figure 2).

Query translator

Result decryptor

Meta Data +

Query

Transformed query

Encrypted results

Results

Server Layer User layer The user layer stores some metadata such as keys and database schemes to translate queries and decrypt query results.

For ease of presentation, we do not distinguish between the human user (DO or DU) and the user layer. When we say "user", we mean the human user who has access to the user layer and similarly, when we say "CSP", we mean the cloud service provider who has access to the server layer.

B. Security Threats

We assume the following threat model.

• The CSP is honest but curious (semi-honest) and has full access to encrypted data. We limit the attacks by the CSP to passive attacks. It is assumed that the CSP executes submitted queries honestly on encrypted data and sends complete and correct results. However, the CSP is curious and may attempt to increase its knowledge about confidential data.

• The CSP has some prior knowledge about outsourced data. For example, the CSP might possess some statistical information about data, such as the domain of plaintext values, minimum or maximum values, or some information about the frequency of values. • We assume the user's machine is trusted. The communication channel between the user and the CSP uses a standard secure protocol such as SSL or IPsec.

C. System Goals

In order to preserve the privacy of outsourced data, we have the following goals:

• Data confidentiality: ensures that unauthorized parties and the CSP have no access to sensitive data. • Index privacy: the proposed indexing method leaks no information about underlying plaintexts. • Efficiency: the proposed method should be efficient for query processing.

IV. OPI In our new order-preserving scheme, called OPI, the order of plaintext values is preserved in their corresponding ciphertexts, but not the frequency of values.

The intuition of our scheme is simple. First, we split the original plaintext's domain into successive intervals of equal lengths. Secondly, we extend the plaintext domain into a new domain. Finally, we use a mapping function to map the original values into the extended domain.

OPI satisfies two properties. 1) Order-preserving: Let D be the domain of plaintext values. The order of values must be preserved in their corresponding indexes, i.e., ∀v, v ∈ D :

if v < v =⇒ OP I(v) < OP I(v ).
2) Uniform distribution: The distribution of plaintexts must be flattened in their corresponding indexes, i.e., two equal values have different indexes with a high probability.

The second property implies that the scheme is not deterministic and improves the robustness against statistical attacks.

We design our scheme for numerical values. Other data types should be translated into integers before indexing [START_REF] Liu | Programmable order-preserving secure index for encrypted database query[END_REF].

A. Notations

Consider a database consisting of a relation T . Suppose that T has one attribute, A. T stores N records. We denote by v i the i th value of attribute A. The goal here is storing the encryption of T , T , at the CSP's.

• Let D be the domain of plaintext values of A, consisting of t distinct values {v 1 , . . . , v t } with corresponding frequencies {f 1 , . . . , f t } where v 1 < v 2 < . . . < v t . • Let D =[l , r ) be an extension of D, the domain of indexes, where l , r ∈ Z. Without lose of generality, we consider l = 0. • Let F(i) be a function that outputs a random number in interval d i = [l i , r i ). OPI consists of the following steps.

• Splitting the domain of indexes: The first step is splitting the domain of indexes D into t intervals

d i (1 ≤ i ≤ t) [19] such that d i =[l i , r i )=[( i-1 k=1 f k ). |D | N , ( i k=1 f k ). |D | N ).
Such splitting has the following properties [START_REF] Liu | New order preserving encryption model for outsourced databases in cloud environments[END_REF], [START_REF] Jalili | Secure data outsourcing based on threshold secret sharing: Towards a more practical solution[END_REF]:

1) d i ⊆ D (1 ≤ i ≤ t) 2) t i=1 d i = D 3) ∀d i , d j i = j: d i ∩ d j = Ø. • Mapping: Second, each plaintext value v i (1 ≤ i ≤ t) is mapped to an index I i , where I i =F(i). In other words, each time v i is repeated, F(i) returns a random number in interval d i =[l i , r i ) as I i .

B. Table Structure

To store relation T at the CSP's, column A is presented in two columns Enc(A) and SIdx(A) (Secure orderpreserving Index) in encrypted table T . The values of Enc(A) are encrypted values of A with an additive homomorphic encryption scheme such as Paillier's cryptosystem [START_REF] Paillier | Public-key cryptosystems based on composite degree residuosity classes[END_REF]. Values I i i = 1, ..., N are stored in column SIdx(A). When a range or exact match query is issued, it is rewritten and executed on column SIdx(A). When an aggregation query is issued, it is directly computed over Enc(A).

C. Query Processing

In this section, we discuss how to process queries. Each query is transformed into a form that can be computed by the CSP. We classify queries into two main classes.

1) Range and Exact Match Queries:

• QUERY-I: "A = v i " is an exact match query, e.g., SELECT * FROM T WHERE A=v i . This query is translated into SELECT * FROM T WHERE SIdx(A) ≥ l i AND SIdx(A) < r i , where l i and r i are the boundaries of interval d i .

• QUERY-II: "A ≤ v i " is a range query that is translated into SELECT * FROM T WHERE SIdx(A) < r i .

• QUERY-III: "A < v" is a range query that is translated into SELECT * FROM T WHERE SIdx(A) < l i .

• QUERY-IV: "v i ≤ A ≤ v j " is translated into SELECT * FROM T WHERE SIdx(A) ≥ l i AND SIdx(A) < r j . l i and r j are the left and the right boundary of d i and d j , respectively.

2) Update Queries: In our scheme, the execution of update queries is straightforward. A deletion query such as

DELETE FROM T WHERE A = v i , is translated into DELETE FROM T WHERE SIdx(A) ≥ l i AND SIdx(A) < r i .
To insert a new value v i , the corresponding encryption and index, Enc(v i ) and I i are computed and inserted into T .

For an update query, tuples satisfying the query predicate are returned. Then, the values in the SET predicate are substituted with the new ones.

V. ENHANCING SECURITY

In this section, we introduce waOPI, a light modification of OPI that improves its security. waOPI reveals less information about underlying plaintexts without sacrificing efficiency. OPI Actually reveals the location of plaintext. waOPI is not order-preserving per se, but still allows range and exact match queries to be processed. This modification is simple and generic: we wrap-around the domain of indexes.

A. waOPI: wrap-around OPI First, a random value, wa, is chosen randomly from D such that l < wa < r . Then, in waOPI the values of wrap-around index waI, waI i i = 1, ..., N , are defined based on algorithm1.

Algorithm 1 Calculating waI

Input: wa Output: waI i , i = 1, ..., N for i = 1 to N do if (i > wa) then waI i-wa ← I i else waI i+N -wa ← I i end if return waI i end for
B. Query Processing waOPI processes range and exact match queries in the same way as OPI. The only difference is that when a query is issued, it must be transformed into a form that can be executed over wrap-around indexes.

Considering a query that asks all data in interval [v i , v j ). This interval is transformed into [l i , l j ), as described in section IV-C1. Then, [l i , l j ) is transformed into [c i , c j ) where:

• if l i ≤ l j ≤ wa then [c i , c j ) = [waI l i +N -wa , waI l j +N -wa ) • if wa ≤ l i ≤ l j then [c i , c j ) = [waI l i -wa , waI l j -wa ) • if l i ≤ wa < l j then [c i , c j ) = [waI N -wa+l i , waI N ] ∪ (0, waI l j -wa ).

VI. SECURITY DISCUSSION

In OPI, we uniformly distribute original values across a domain into another broader domain. Having a uniform distribution of values minimizes information leakage to an adversary who would observe the outsourced data. Overhead comes from the new domain size, which is bigger than the original domain's. Theoretically, we need domain D to be of size |D | |D| to uniformly distribute all original values across D . Choosing an appropriate domain size for D is a trade-off between security and storage cost [START_REF] Liu | New order preserving encryption model for outsourced databases in cloud environments[END_REF], [START_REF] Jalili | Secure data outsourcing based on threshold secret sharing: Towards a more practical solution[END_REF].

waOPI provides a substantial improvement in security. Given the wrap-around indexes, the adversary cannot deduce any information about the real location of underlying plaintexts.

However, our schemes are possibly vulnerable when executing queries. The adversary or the CSP is able to infer the frequency of attribute values by observing frequent intervals accessed at query time. The adversary might also be able to infer random value wa. Yet, it has been proven that no scheme is secure against an adversary that has access to both plaintexts and prior knowledge about query workloads [START_REF] Wong | Secure query processing with data interoperability in a cloud database environment[END_REF], [START_REF] Boldyreva | Orderpreserving encryption revisited: Improved security analysis and alternative solutions[END_REF]. To prevent such inference, Obvious RAM (ORAM) [START_REF] Goldreich | Software protection and simulation on oblivious rams[END_REF], [START_REF] Shi | Oblivious RAM with o((logn)3) worst-case cost[END_REF], [START_REF] Stefanov | Oblivistore: High performance oblivious cloud storage[END_REF], [START_REF] Stefanov | Path ORAM: an extremely simple oblivious RAM protocol[END_REF] can be used. ORAM hides access patterns during query processing, but ORAM techniques are currently computationally prohibitive and cannot be used in practice. Another solution is sending fake queries. Fake queries must be frequently submitted. The goal is misleading the adversary to infer a wrong wa. Submitting fake queries makes waOPI more robust.

VII. RELATED WORKS

OPE was introduced in the database community by Agrawal et al. [START_REF] Agrawal | Orderpreserving encryption for numeric data[END_REF] as a tool to support efficient range queries over encrypted data [START_REF] Boldyreva | Orderpreserving encryption revisited: Improved security analysis and alternative solutions[END_REF]. This scheme maps each value of the plaintext domain to one value in the ciphertext domain, such that the relative ordering of plaintexts is preserved after encryption. This scheme bears weak privacy protection because the CSP can statistically estimate the original data values [START_REF] Liu | New order preserving encryption model for outsourced databases in cloud environments[END_REF]. Boldyreva et al. define the first formal ideal-security definition of OPE, called indistinguishability under ordered chosen-plaintext attack (IND-OCPA) [START_REF] Boldyreva | Orderpreserving symmetric encryption[END_REF]. Informally, an IND-OCPA scheme reveals no additional information about the underlying plaintexts beside their order, which is the minimum requirement of the order-preserving property. IND-OCPA [START_REF] Boldyreva | Orderpreserving symmetric encryption[END_REF] introduces a random mapping that preserves order, but leaks at least half of the plaintext bits (i.e., more information than OPE) and shows poor efficiency [START_REF] Popa | An idealsecurity protocol for order-preserving encoding[END_REF].

Popa et al. introduce the first practical IND-OCPA scheme, mutable order-preserving encoding (mOPE) [START_REF] Popa | An idealsecurity protocol for order-preserving encoding[END_REF]. mOPE requires an interactive protocol for query processing. Additionally, mOPE relies on user defined functions (UDFs) for query processing, which makes it unsuitable for cloud outsourcing.

Liu et al. introduce an OPE scheme that randomly splits the original plaintext domain into successive intervals with different lengths [START_REF] Liu | New order preserving encryption model for outsourced databases in cloud environments[END_REF]. Then, an extended ciphertext domain is selected and split into the same number of intervals. Finally, nonlinear mapping functions map the original plaintexts into ciphertexts in the extended domain.

While Liu et al.'s scheme is efficient for query processing, it partially destroys the distribution of original data because the splitting method is random and does not rely on data frequency.

Table I shows the comparison between existing OPE schemes. It highlights that mOPE has the lowest performance because it requires the user to interact with the server to retrieve requested data. Our schemes do not require such interaction and present the same efficiency as Liu et al.'s.

In IND-OCPA approach, only the order of data is revealed, while our scheme reveals nothing about the original data. Compared with Liu et al.'s scheme, waOPI can achieve security against statistical attacks, because we use a density and frequency-aware domain splitting method.

VIII. CONCLUSION

In this paper, we propose an order-preserving indexing scheme, OPI, in which the ordering relation between plaintexts is preserved in their corresponding indexes. In OPI, the distribution of index values is different from the plaintext distribution. We extend our scheme in order to achieve robustness against statistical attacks. In waOPI, the security of plaintexts is indeed improved by hiding the rule of data distribution and data frequency, without sacrificing efficiency.

Nevertheless, we should prove the security of OPI and waOPI. We also plan to investigate secure indexing mechanisms with additive homomorphic property to efficiently support aggregation queries over encrypted data.

Figure 1 .

 1 Figure 1. Comparison between plaintext and ciphertext distribution. (a) Plaintext distribution. (b) Ciphertext distribution of deterministic OPE. (c) Ciphertext distribution of one-to-many OPE [17].

Figure 2 .

 2 Figure 2. System model.