Computation of Ground States of the Gross-Pitaevskii Functional via Riemannian Optimization

Computation of Ground States of the Gross-Pitaevskii Functional via Riemannian Optimization

Ionut Danaila and Bartosz Protas

Laboratoire de mathématiques Raphaël Salem
Université de Rouen
Technopôle du Madrillet
76801 Saint-Étienne-du-Rouvray, FRANCE

Department of Mathematics & Statistics,
McMaster University
Hamilton, Ontario L8S4K1, CANADA
July 5, 2019

In this paper we combine concepts from Riemannian Optimization [P.-A. Absil, R. Mahony, R. Sepulchre, Optimization algorithms on matrix manifolds, Princeton University Press, 2008] and the theory of Sobolev gradients [J. W. Neuberger, Sobolev gradients, Springer 2010] to derive a new conjugate gradient method for direct minimization of the Gross-Pitaevskii energy functional with rotation. The conservation of the number of particles in the system constraints the minimizers to lie on a Riemannian manifold corresponding to the unit norm. The idea developed here is to transform the original constrained optimization problem to an unconstrained problem on this (spherical) Riemannian manifold, so that faster minimization algorithms can be applied. First, we first obtain Sobolev gradients using an equivalent definition of an inner product which takes into account rotation. Then, the Riemannian gradient (RG) steepest descent method is derived based on projected gradients and retraction of an intermediate solution back to the constraint manifold. Finally, we use the concept of the Riemannian vector transport to propose a new Riemannian conjugate gradient (RCG) method for this problem. It is derived at the continuous level based on the “optimize-then-discretize” paradigm instead of the usual “discretize-then-optimize” approach, as this ensures robustness of the method when adaptive mesh refinement is performed in computations. Numerical tests carried out in the finite-element setting based on Lagrangian piecewise quadratic space discretization demonstrate that the proposed RCG method outperforms the simple gradient descent RG method in terms of rate of convergence. The RCG method is extensively tested by computing complicated vortex configurations in rotating Bose-Einstein condensates, a task made challenging by large values of the non-linear interaction constant and the rotation rate.

Keywords: Gross-Pitaevskii functional, Bose-Einstein condensate, Riemannian optimization, Sobolev gradients, conjugate gradients

1 Introduction

The rotating Bose-Einstein condensate (BEC) represents a highly-controllable quantum system offering an ideal framework to study quantized vortices at a macroscopic level. A rich variety of vortex states, from a single vortex line to dense Abrikosov vortex lattices and giant vortices, were experimentally observed and extensively studied in the last two decades (e. g. [50]). A standard mathematical approach to describe equilibrium configurations with quantized vortices in rotating BEC is the minimization of the Gross-Pitaevskii (GP) energy functional with rotation [44, 38]. In addition to the global minima, the so-called “ground states”, local minima are also of interest as they represent excited, or meta-stable, states which are more likely to be observed in experiments [20]. The minimizers have the form of complex-valued wavefunction fields dependent on the space variable, resulting in an infinite-dimensional minimization problem. The complexity of the minimization problem is further compounded by a constraint imposed on the norm of the minimizers which reflects the conservation of the number of atoms in the condensate. For the mathematical properties of the GP energy with rotation and the corresponding ground states we refer the reader to [3, 6, 38, 15].

In this paper we address the problem of direct minimization of the GP energy functional with rotation when large nonlinear interaction constants and high rotation frequencies are considered. A number of approaches to direct minimization of the GP energy have been proposed based on various standard and emerging mathematical methods for optimization problems in finite dimension: Optimal Damping Algorithm [28, 35], Newton-like method based on Sequential Quadratic Programming (SQP) [21], Interior Point Method (Ipopt) [51], Inertial Proximal Algorithm (iPiano) [10] and regularized Newton method with trust region [52]. Alternative approaches which do not involve direct minimization of energy rely instead on the solution of the corresponding Euler-Lagrange system which has the form of a nonlinear eigenvalue problem. In the latter context, a wide variety of classical integration and iterative techniques have been employed such as Newton’s method [18], Runge-Kutta [22] and continuation methods [23], etc.

Another class of approaches, pioneered in [16], relies on a normalized gradient flow for the GP functional and became popular due to their efficiency and ease of implementation (see also the review papers [13, 15, 9, 14]). These methods consist in first solving the gradient flow equation for the minimization of an unconstrained energy followed by a normalization of this “predictor” solution to bring it back to the constraint manifold. Solution of the gradient flow equation can be viewed as a pseudo-time (or imaginary time) integration of the time-independent GP equation. Discretization of the gradient-flow equations using a (natural) steepest descent method would result in a very inefficient explicit Euler scheme for the (imaginary-) time integration. For this reason in [16] the gradient-flow equation was solved using a semi-implicit backward Euler scheme which proved even more efficient than the classical Crank-Nicolson scheme. The convergence of the original scheme suggested in [16] was recently improved in [11, 12] by using different discrete preconditioners. It is interesting to note that the gradient-flow equation for the GP functional has structure similar to the complex-valued heat equation which makes it amenable to solution with different classical time-integration schemes such as Runge-Kutta-Fehlberg [31], backward Euler [16, 6, 17], second-order Strang time-splitting [16, 6], combined Runge-Kutta-Crank-Nicolson scheme [4, 5, 25], etc.

As regards the development of numerical methods, there are two main paradigms, namely, “optimize-then-discretize” and “discretize-then-optimize”, depending on whether the gradient expressions are derived at the continuous or discrete level. In the first case, the Sobolev gradient approach [41] represents the gradient-flow method formulated with respect to a judiciously selected inner product in a Hilbert space, rather than the classical inner product. The required gradients are obtained via the Riesz representation theorem. When discretized, the Sobolev gradient approach can be also interpreted as suitable preconditioning applied to the gradient [29]. However, the key advantage of working with the “optimize-then-discretize” formulation is that the form of this preconditioning is dictated by the functional (Sobolev) setting of the problem and thus avoids the technically complicated search for a good (and mesh-independent) discrete preconditioner. Sobolev gradient methods were successfully applied to minimize the GP energy in the presence of rotation in [31, 27].

The purpose of this contribution is to develop and validate an efficient computational approach to minimization of the GP energy by combining the Sobolev gradient method with concepts of Riemannian optimization [1, 48]. This allows us to transform the original constrained optimization problem to an unconstrained problem on a Riemannian manifold with a very simple structure which is amenable to solution using the conjugate gradient approach. We begin by formulating a Riemannian version of the Sobolev gradient approach in which the “retraction” operation ensures the norm constraint is satisfied at all discrete times. Then, convergence is accelerated using a Riemannian version of the conjugate gradient method which relies on the notion of the “vector transport” applied to the gradient and the descent direction. Such approaches are already well established in the context of problems formulated in finite dimensions [19], but have received rather little attention in the context of problems in infinite dimension (an exception being the study [43], see also [7, 8, 40]). To the best of our knowledge, they represent a new direction as regards minimization of the GP energy. We demonstrate that in combination with a flexible finite-element discretization involving adaptive grid refinement [26, 51], this proposed novel approach outperforms standard techniques.

The structure of the paper is as follows: in the next section we state the mathematical model describing minimization of the GP energy; in §3 we recall the Sobolev gradient method with its projected gradient variant, whereas the Riemannian gradient and conjugate gradient methods are introduced in §4.2; numerical discretization based on adaptive finite elements and its software implementation are discussed in §5; in §6 we use the method of the manufactured solutions to estimate the rates of convergence of the different approaches; in §7 we compute a number of challenging BEC configurations with vortices in various arrangements; discussions and conclusions are deferred to §8.

2 Mathematical Model

The energy of a rotating homogeneous BEC at zero temperature is given in terms of the Gross-Pitaevskii (GP) energy functional [44, 38]. After applying standard scaling and dimension reduction [44, 15], its non-dimensional form defined here on a two-dimensional domain becomes


where is a complex-valued wavefunction and its complex conjugate, , is the trapping potential. and are real constants characterizing the strength of the nonlinear interactions and rotation frequency, respectively. The wavefunction vanishes outside the trap and is therefore assumed to satisfy the homogeneous Dirichlet boundary conditions on . The conservation of the number of atoms in the condensate is expressed as


and serves as a constraint on . For the energy functional (1) to be well-defined, the wavefunction must belong to the Sobolev space of functions with square-integrable gradients [2] and vanishing traces on the boundary (precise definitions of the norms in this function space will be provided in §3 below). The constraint (2) may now be interpreted as defining a manifold in the solution space, i. e.


We assume the trapping potential to have the following general form allowing us to represent different trapping potentials used in experiments


for some . Along with the energy (1), another important integral quantity describing the rotating BEC is the total angular momentum


where denotes the real part of a complex number.

Global minimizers of the energy functional (1) defined through


are called ground states. Local minimizers, with energy larger than that of the ground state, are referred to as excited or meta-stable states. The Euler-Lagrange system corresponding to the minimization problem (6) is derived using standard techniques and leads to the stationary Gross-Pitaevskii equation


The ground state and excited states are therefore eigenfunctions of the nonlinear eigenvalue problem (7).

3 Gradient Flows and Steepest Descent Sobolev Gradient Methods

Numerical techniques for the solution of optimization problem (6) can be derived from the general form of the corresponding gradient-flow equation


where is an initial guess and represents the gradient of the GP energy functional (1) at , computed with respect to the topology of the Hilbert space (to be made specific below). The gradient-flow needs to be additionally constrained to ensure that for . As shown below, many different computational approaches can be derived from (8) by making specific choices of (i) the Hilbert space , (ii) discretization of the initial-value problem (8) with respect to pseudo-time and (iii) how the constraint is imposed.

As regards the definition of the gradient , the most common choices of the function space are and the Sobolev space with the inner products


where is the inner product for complex numbers. A new inner product equivalent (in the precise sense of norm equivalence) to (9b) was suggested in [27]


and will be adopted in our considerations below. Its definition was motivated by the following physically more revealing form of the energy functional equivalent to (1)


where the effective non-dimensional trapping potential is obtained from the original potential by subtracting a term representing the centrifugal force [49]


Different gradients can be obtained by making different choices of the Hilbert space and using the Riesz representation theorem [39]. For each , one can find an element of denoted , such that the directional Gâteaux derivative of the energy at in the direction is expressed as , . We call such an element of the gradient of at . The gradient can be obtained from (1) in a straightforward manner as


Then, the gradient, denoted , is obtained using the Riesz representation theorem together with (10) and (13) as the solution of an elliptic boundary-value problem, which we state here in the weak form


3.1 Normalized gradient flow

The gradient flow equation (8) with the gradient (13) was discretized in [16] using a semi-implicit backward Euler (BE) method


where denotes the approximation obtained at the th discrete time level, is an intermediate (predictor) field and is a fixed (imaginary) time step. The approximation at the time level has to satisfy the unit-norm constraint (2) and is therefore obtained by normalizing the predictor solution as


This approach is referred to as the normalized gradient flow method (see also [13, 15, 9, 14]). Different existing variants of this method use different numerical approaches to integrate the gradient-flow equation (8), e. g. Runge-Kutta methods [31, 4, 5, 25], different implicit schemes [16, 6, 17] and Strang-type time-splitting approaches [16, 6]. Even though some of these schemes do possess the energy-decreasing property, they typically do not preserve the gradient-flow structure at the discrete level. Therefore, such imaginary-time methods can be regarded as solving the nonlinear eigenvalue problem (7), rather than directly minimizing the GP energy. Another potential drawback of such approaches is that solutions of (7) are in general critical points of the energy which are not necessarily minima.

3.2 Steepest Descent Sobolev Gradient Methods

Hereafter we focus on techniques which do preserve the gradient-flow structure of (8) on the discrete level while explicitly accounting for the presence of the unit-norm constraint (2). As a starting point, we will thus consider an explicit discretization of (1) in the following generic form


where is a suitable descent step size, whereas is a Sobolev gradient defined for or . Below we discuss two ways in which the information about the constraint can be incorporated in the gradient method.

3.2.1 Projected Sobolev Gradient Method

By considering the following identity derived from (17)


we note that using an unconstrained gradient leads to an error in the satisfaction of the constraint (2) at each iteration. Normalization of the solution is then necessary to bring it back onto the manifold (see Figure 1a).

Figure 1: Schematic illustration of the principle of the steepest descent method on the spherical manifold for (a) the simple (“unprojected”) gradient method and (b) the projected gradient (PG) method. In case (b), since the gradient belongs to the subspace tangent to the manifold at , normalization is equivalent to Riemannian retraction (24).

This error can be reduced to second order by requiring that , which is achieved by projecting the gradient onto the subspace


tangent to the constraint manifold at . As was shown in [27, 37], the associated projection operator can be expressed in the following general form


where is a solution of the variational problem


We note that if , in (21) and we recover the well-known explicit expression of the projected -gradient (e. g. [7]).

Hereafter we will set and denote . Replacing with in (17), we obtain the projected gradient (PG) method suggested in [27]


While in [27] a fixed step size was used, here we use an optimal step size found through the solution of a line-minimization problem


An explicit expression for the optimal decent step was derived in [51] based on a particular form of the GP energy. In this study, we prefer to solve problem (23) with a general line-minimization approach such as Brent’s algorithm [45, 42] as it has the advantage of being easily adapted to the Riemannian gradient methods presented in the next section. To mitigate the drift away from the constraint manifold allowed by the (PG) iterations, normalization analogous to (16) may be applied to the iterates after a certain number of steps. The idea of the (PG) approach is illustrated schematically in Figure 1b.

4 Riemannian Optimization

In this section we discuss some basic concepts relevant to optimization on manifolds, known as Riemannian optimization [1, 48]. In contrast to the perspective developed in the previous section, here we pursue a different, “intrinsic” approach where optimization is performed directly on the manifold . The main advantage of such a formulation is that it allows one to treat (6) as an unconstrained optimization problem opening it to application of algorithms such as the method of conjugate gradients which are more efficient than the simple (projected) gradient approach.

In addition to the definition of the projection on the tangent space already introduced above, cf. (20)–(21), we need to introduce two more concepts, namely, the “retraction” (also referred to as “exponential mapping”) and the associated “vector transport”. While in general these operators can have a rather complicated form, in the present problem where the constraint manifold is given by (3), they are reduced to fairly simple expressions. We refer the reader to the monograph [1] for additional details concerning the differential–geometry foundations of this approach.

4.1 Riemannian Gradient Method

Given a tangent vector , where is a state on the manifold, the retraction is defined as the operator


where the norm used in the denominator is the same as the norm defining the constraint manifold in (3). We note that for the spherical manifold the retraction operator is equivalent to normalization (16) already used in the previous sections.

The Riemannian gradient (RG) method is then obtained applying relation (24) to the projected gradient , cf. (20)–(21), used in the (PG) approach which yields


The step size is found optimally by solving a generalization of the line-minimization problem (23) which uses retraction (24) to constrain the samples to manifold , i. e.


We refer to problem (26) as “arc-minimization”. We solve this problem using a straightforward modification of Brent’s algorithm [45, 42]. In addition to application of retraction (24) at every iteration in the latter case, the key difference between the (PG) and (RG) approaches lies in how the optimal step size is determined, cf. (23) vs. (26). The idea of the (RG) method is illustrated schematically in Figure 1b.

4.2 Riemannian Conjugate Gradient Method

As a point of reference, we begin by recalling the nonlinear conjugate gradients method in the Euclidean case. Given a function , this approach finds its local minimum as , with the iterates defined as follows [42]


where is the initial guess and the descent direction is constructed as


in which and is a “momentum” term chosen to enforce the conjugacy of the search directions , . When the objective function is convex-quadratic, i. e. for some positive-definite matrix , approach (27)–(28) reduces to the “linear” conjugate gradient method in which and are given in terms of simple expressions involving [42]. In the non-quadratic setting, which is the case of problem (6), the step size needs to be found via line minimization as described by (26), whereas the momentum term is typically computed using one of the following expressions


where is the inner product defined with respect to the metric (in the simplest case when , for ). The coefficient is periodically reset to zero (say, every 10 or 20 iterations) which is required to ensure convergence for convex, non-quadratic problems [42]. It is well known that for optimization problems which are locally quadratic the conjugate gradient approach exhibits faster (though still linear) convergence than the convergence characterizing the simple gradient method, especially for poorly scaled problem [42]. Similar observations have been also reported for Riemannian optimization problems [1, 48].

We explain below how the conjugate gradients approach can be adapted to the Riemannian case involving the energy functional (1) defined for infinite-dimensional state variables . There are two key issues which must be addressed:

  1. the two term on the right-hand side (RHS) of formula (28) belong to two different linear spaces which are the tangent spaces constructed at two consecutive iterations, i. e. and ; as a result, they cannot be simply added; the same problem also concerns the inner-product expressions in the numerator of the Polak-Ribière momentum term (29b),

  2. while in the finite dimension all norms are equivalent, this is no longer the case in the infinite-dimensional setting where the choice of the metric does play a significant role; in our approach, although the gradient-descent equations are discretized for the purpose of the numerical solution, their specific form is derived in the infinite-dimensional setting (in other words, we follow the “optimize-then-discretize” paradigm [32]); in addition to the momentum term (29), the choice of the metric implied by the inner product also plays a role in the construction of the projection (20)–(21) and the vector transport which will be defined below.

The key concept required in order to address the first issue is the vector transport


where is the tangent bundle, describing how the vector field is transported along the manifold by the field [1]. It therefore generalizes the concept of the parallel translation to the motion on the manifold and is also closely related to the “affine connection” which is one of the key differential-geometric quantities characterizing a manifold. The vector transport thus provides a map between the tangent spaces and obtained at two consecutive iterations, so that algebraic operations can be performed on vectors belonging to these subspaces.

In general, vector transport is not defined uniquely and in the present case when the manifold is a sphere, the following two definitions lead to expressions particularly simple from the computational point of view. Let and ; the transport of by can be expressed either by

  • vector transport via differentiated retraction

  • or by vector transport on Riemannian submanifolds


where is the orthogonal projector on . We note that these formulas differ only by a scalar factor and the vector transport is linear in the field , but not in . We refer the reader to monograph [1] for details concerning the derivation of formulas (31)–(32). The numerical results presented in §6 and §7 are obtained using expression (32) with the inner product and norm.

With the use of the vector transport (32), the conjugate gradient method (27)–(29) can be rewritten in the Riemannian infinite-dimensional setting as




with the Polak-Ribière momentum term modified as follows (the corresponding term in the Fletcher-Reeves approach remains unchanged)


The optimal descent step in (33) is computed as in (26) by solving the corresponding arc-minimization problem


using a generalization of Brent’s method. We refer to approach (33)–(36) as the Riemannian Conjugate Gradients (RCG) method and its idea is schematically illustrated in Figure 2.

Figure 2: Schematic illustration of the principle of Riemannian conjugate gradients (RCG) method on a spherical manifold. (a) Riemannian vector transport of the anterior conjugate direction ; the transport of the anterior gradient is performed in a similar way. (b) Projection of the new Sobolev gradient onto the tangent subspace resulting in . The linear combination (34) of and the transported anterior direction is computed in .

5 Space Discretization

We use a finite-element approximation constructed as follows. Let be a family of triangulations of the domain . We assume that is a regular family, with belonging to a generalized sequence converging to zero. We denote by the space of polynomial functions of degree not exceeding defined on triangles . We also introduce the finite-element approximation spaces


The finite dimensional space is a subspace of and therefore will be used to approximate the energy functional (1) and the different expressions representing gradients and descent directions in algorithms (PG), (RG) and (RCG). In the following we use (, piecewise quadratic) finite elements to approximate the solution and a representation of the nonlinear terms. In addition, adaptive mesh refinement suggested in [26] and tested in [26, 51] is used to adapt the grid during iterations leading to a significant reduction of the computational time. The approach is implemented in FreeFEM++ [33, 34]. The implementation of the Riemannian retraction (24) and vector transport (32) is straightforward and was found to work very well with arc-minimization (36) and adaptive mesh refinement.

6 Convergence Rates of Different Gradient Methods

To assess the convergence rate of the minimization algorithms (PG), (RG) and (RCG), we use the method of manufactured solutions [46] which is widely used for verification of calculations. By introducing an extra source term, the original system of equations is modified to admit an exact solution given by a convenient analytic expression. Even though in most cases exact solutions constructed in this way are not physically realistic, this approach allows one to rigorously verify computations. Here we manufacture such an exact solution in the form


where are the cylindrical coordinates of the point and is the radius of the circular domain . We note that this solution satisfies constraint (2) and qualitatively resembles a giant vortex in the condensate (see Figure 3 and §7).

Figure 3: Manufactured solution (39) visualized with a 3D-rendering of the modulus color-coded with (a) the modulus itself and (b) the phase of the solution for .

It also satisfies an inhomogeneous form of the nonlinear problem (7), i. e.


and is a critical point of the modified energy functional


For this energy functional, the gradient is expressed as in (13), but with a supplementary term added. Given the form (39) of the manufactured solution and assuming a harmonic trapping potential , from (40) we obtain , where is a polynomial of degree 9. From this and relations (41) and (5) we can deduce exact expressions for the energy and the angular momentum .

The numerical tests are based on the manufactured solution (39) corresponding to the following parameter values: , , , , , where, to make the problem more challenging, large values of the nonlinear interaction constant and rotation frequency are used (cf. Figure 3). Minimization of energy (41) subject to the constraint of unitary norm is performed using the gradient approaches (PG), (RG) and (RCG) which rely on Sobolev gradients defined with respect to the inner product (10), i. e. in (20)–(21). The step size is determined at each iteration by solving the line/arc-minimization problem, cf. (23), (26) and (36).

In order to assess the mesh-independent effect of the Sobolev gradient preconditioning, cf. §5, we perform computations using two grids: Mesh 1 consisting of 24,454 triangles with and Mesh 2 consisting of 99,329 triangles with , where is the smallest grid size. In the present case no mesh refinement was performed during iterations. The initial guess is taken as , whereas iterations are declared converged once the following condition based on the relative energy decrease is satisfied [26, 51]


To assess the rate of convergence of the (PG), (RG) and (RCG) approaches, the quantities and are shown as functions of iterations for the two spatial discretizations in Figures 4a and 4b. In these figures we observe linear convergence followed by a slower convergence at final iterations. The change of the slope of error curves occurs at the level at which the minimization errors are comparable to the errors related to the spatial discretization. In other words, in the “optimize-then-discretize” setting adopted here, gradient expressions derived based on the continuous formulation, cf. (13)–(14), may no longer accurately represent the sensitivity of the discretized objective function, when the difference between and is of the order of the space discretization errors.

Figure 4: Test case based on the manufactured solution (39). Convergence of the (PG), (RG) and (RCG) methods: (a) and (b) are shown as functions of iterations for different discretizations. Dashed and dash-dotted lines indicate the least-squares fits (43).

In Figures 4a and 4b it is evident that the (RCG) method converges much more rapidly (39 iterations) than the (RG) approach (180 iterations). As expected, the convergence of the (PG) method (202 iterations) is similar to that of the (RG) method. To quantify the convergence rates we use the following ansatz to represent the errors:


The values of the parameters and , which represent the factors by which the corresponding errors are reduced between two iterations can be obtained from least-squares fits of the data in Figures 4a,b in the linear regime. These results are collected in Table 1 (the corresponding fits are also indicated in Figures 4a,b). First, these results demonstrate that the rate of convergence is grid-independent as expected from the general theory of Sobolev-gradient descent methods [41]. The data in Table 1 can also be interpreted in terms of the classical theory of conjugate gradient method [42], which for the minimization of quadratic functions predicts that . We see that the data from Table 1 satisfies this relationship with the accuracy of a few percent. For the simple gradient and conjugate gradient methods we furthermore have the approximate relationships and , respectively, where is the condition number characterizing the problem. Using the data from Table 1 (Mesh 2), we infer that for the gradient (RG) and for the conjugate gradient (RCG) method, indicating that the convergence acceleration produced in the present problem by the Riemannian conjugate gradient approach actually exceeds what can be expected from the standard theory.

Mesh 1 Mesh 2
(RG) 0.9167 0.9574 0.9496 0.9268 0.9627 0.9538
(RCG) 0.2909 0.5394 0.5275 0.2924 0.5408 0.5238
Table 1: Parameters characterizing the least-squares fits (43) of the data shown in Figure 4.

Since the exact solution is usually unavailable in physically relevant problems, we now verify that convergence of iterations can be monitored based on quantities which do not involve . Indeed, the evolution of and with iterations shown in Figures 5a and 5b exhibits the same trends as the data shown in Figures 4a and 4b, except for the slowdown observed in the latter case. This demonstrates that either of these two quantities can be used to monitor convergence and, in particular, check the stopping criterion (see also [16, 13, 15, 9, 14]).

Figure 5: Test case based on the manufactured solution (39). Convergence of the different quantities with iterations : (a) , (b) , (c) , cf. (44), drift away from the constraint manifold and (d) , the optimal descent step.

There are two additional aspects of the convergence of the different methods we wish to comment on. In Figure 5c we show the evolution of the “drift” away from the constraint manifold exhibited by the intermediate approximations before retraction (24) is applied


This quantity measures how far the intermediate steps diverge from the constraint manifold. We see that, as compared to the (PG) and (RG) methods, in the (RCG) approach the intermediate approximations always remain closer to . Finally, the step size determined by the gradient approaches (PG), (RG) and (RCG) via line/arc-minimization, cf. (23), (26) and (36), is shown in Figure 5d. We see that the step sizes generated by the simple gradient methods, (PG) and (RG), tend to oscillate between two values, a behavior indicating that the iterations are trapped in narrow “valleys”. This is a common behavior of the steepest descent method when applied to poorly conditioned problems and is not exhibited by the the (RCG) iterations where on average the steps also tend to be longer. Performance of these different methods applied to several realistic problems will be discussed in the next section.

7 Computation of Rotating Bose-Einstein Condensates

In this section we compare the performance of the minimization algorithms (PG), (RG) and (RCG) on a number of test cases involving configurations of rotating BEC with vortices. We consider increasingly complex arrangements ranging from a single vortex to Abrikosov vortex lattices with more than one hundred vortices and to giant vortices. In some cases, we also compare the gradient methods with the (BE) method (15)–(16) implemented using the same finite-element setting. To make these test cases more challenging, we consider large values of the nonlinear interaction constant and large angular frequencies .

7.1 Test Case #1: BEC with a Single Central Vortex

We consider the case of a BEC trapped in a harmonic potential and rotating at low angular velocities: , , . For this case, the Thomas-Fermi (TF) theory [49] offers a good approximation of the atomic density of the condensate with the effective trapping potential given by (12). By imposing , we can derive analytical expressions for the corresponding approximation of the chemical potential [51]. From the same TF approximation, we can estimate the radius of the condensate as . Consequently, we set up the computational domain as a disk of radius . The initial guess is taken in the form of an off-center vortex placed at , see Figure 6a. We use the ansatz , where with representing the polar coordinates centered at and the non-dimensional healing length, which is a good approximation of the vortex radius in rotating BEC [30]. Similar ansatz will also be used in subsequent sections to set up initial guesses with vortices for the calculation of more complicated BEC configurations. The stopping criterion (42) is used with the value . In these calculations the grid remains fixed (i. e. no grid adaptation is performed). In the (BE) method the imaginary time step is chosen as , which proved to be optimal for convergence after testing values of in the range from to .

Figure 6: Computation of a rotating BEC with a single central vortex (cf. §7.1). 3D rendering of the atomic density for: (a) the initial guess and (b) the converged ground state.

For the considered physical parameters the ground state features a vortex centered at the origin. All considered methods, (PG), (RG), (RCG) and (BE), converged to the same ground state shown in Figure 6b. In order to assess their respective rates of convergence, we compute a reference (“exact”) solution using the same grid and starting the minimization algorithms from the initial guess with . The corresponding energy and angular momentum are and . The quantities , and shown in Figures 7a, 7b and 7c as functions of indicate that while the (RG) and (BE) methods converge with similar rates, the (RCG) approach converges much faster requiring only 100 iterations to satisfy the stopping criterion, as opposed to 1220 iterations for the (BE) method, 1308 for the (RG) method and 1399 for the (PG) method. The drift away from the constraint manifold at intermediate steps , cf. (44), during the first 200 iterations is shown for different methods in Figure 7d. We note that, due to the large step size , this quantity is always for the (BE) method. The normalization step is therefore crucial in this approach. On the other hand, is reduced faster in the (RCG) approach where it also attains lower values than in the (PG) and (RG) methods.

Figure 7: Computation of a rotating BEC with a single central vortex (cf. §7.1). Convergence of the different quantities with iterations : (a) , (b) , (c) and (d) , cf. (44), the drift away from the manifold.

7.2 Test Case #2: BEC with dense Abrikosov vortex lattice

We now move on to consider a more challenging test cases corresponding to a harmonic potential (), high rotation rate () and large values of the nonlinear interaction constant ( to ). We note that for the harmonic trapping potential there is a physical limit occurring at the rotation frequency when the trapping is canceled by the centrifugal force (i. e. , see eq. (12)). The next section will consider cases with a modified trapping potential, allowing for higher rotation frequencies.

We start with the case with and for which the ground state features over fifty vortices arranged in a regular triangular lattice called the Abrikosov lattice. The difficulty here is to obtain a very regular lattice, in particular for the vortices located near the border of the condensate where the atomic density is low. This explains why the results previously reported for this case exhibit configurations with somewhat different arrangements of the peripheral vortices which nevertheless have very similar energy levels [53, 27, 11]. These differences can be attributed to the use of different initial guesses . A nearly perfect arrangement of vortices on a triangular/hexagonal lattice is reported in the recent study [52] and will be considered here as a reference result used to validate our methods (the corresponding energy level is ). Details of the computed stationary states depend on the initial guess and we used three distinct forms of : (i) “ansatz d” proposed in [52] to model a central vortex using Gaussian functions, and the Thomas-Fermi approximation described in §7.1 with (ii) one central vortex and (iii) six vortices. The corresponding stationary solutions are shown in Figures 8a, 8b and 8c, where we can see that the central parts of the vortex lattices are in all cases essentially identical (modulo rotation) and some differences are detected among the peripheral vortices. The values of energy corresponding to these configurations differ by less than .

Figure 8: Computation of a rotating BEC with a dense Abrikosov vortex lattice (cf. §7.2). Stationary states obtained using different initial condition : (a) “ansatz d” suggested in [52], (b) Thomas-Fermi atomic density with one central vortex, (c) Thomas-Fermi atomic density with a ring of six vortices. The figures in the first row show contours of the atomic density (normalized by its maximum value ) and in the second row they show the 3D-rendering of the same contours. The corresponding energies are: (a) , (b) , (c) , to be compared to the reference value from [52].

In the following we carry out computations starting from the initial guess (i) and grid adaptation is now performed during iterations. Convergence of the iterations carried out with the (RG), (RCG) and (BE) approaches is presented in Figures 9a and 9b, where one can see that the conjugate gradient (RCG) method converges after 1046 iterations, based on the stopping criterion (42), whereas the gradient (RG) and the backward-Euler (BE) methods show a markedly slower convergence (they were stopped after 3000 iterations). We note that the peaks in the curves shown in Figures 9a and 9b result from reinterpolation of intermediate solutions after grid adaptation [26, 51]. In these figures we see that in all cases convergence slows down at later iterations which is related to the slow rearrangement of vortices near the boundary of the condensate. Another possible reason is that since in [52] a different discretization was used (Fourier spectral approach with periodic boundary conditions), the value of taken from that reference might not exactly correspond to ours.

Figure 9: Computation of a rotating BEC with a dense Abrikosov vortex lattice (cf. §7.2). Convergence of the different quantities with iterations : (a) and (b) . All computations start from the initial guess suggested in [52] ("ansatz d").

Finally, we use the (RCG) method to compute fast rotating BEC () corresponding to large values of the nonlinear interaction constant with varying from 1,000 to 15,000. For these difficult cases, a more physically relevant assessment of the convergence of iterations is provided by the alignment of vortices on parallel lines inside the vortex lattice. Since isocontours of atomic density do not always coincide with these lines, we developed a post-processing approach to identify the centers of vortices by detecting local minima of the function . This post-processing is similar to that used for experimental data [24] or 3D numerical simulations [25] and also allows to build the Delaunay triangulation of the lattice and compute the radius of each vortex. The resulting stationary states are presented in this way in Figure 10. We notice an arrangement of vortices on a nearly perfect lattice for and , and a less regular arrangement for and with the presence of some defects in the lattice. This effect could be related to physical theories addressing the non-uniformity of the inter-vortex spacing in dense Abrikosov lattices [24, 47].

Figure 10: Computation of a fast rotating BEC () in a harmonic trapping potential (cf. §7.2). The Abrikosov vortex lattice is represented using the Delaunay triangulation built from the detected vortex centers. Configurations obtained for large values of the nonlinear interaction constant: (55 vortices), (134 vortices), (193 vortices) and (237 vortices).

7.3 Test Case #3: BEC with giant vortex

To overcome the limit imposed by the harmonic trapping potential, a modified “harmonic-plus-Gaussian” potential was tested in experiments [20]. In [5, 25] this new experimental set-up was modeled as


with the possibility to switch from a “quartic-plus-quadratic” potential (), which corresponds to experiments, to a “quartic-minus-quadratic” potential (), which is experimentally feasible but was never tested. Adapting the analysis from [5] to our 2D case, we obtain three possible regimes depending on the type of potential:
Regime 1: “quartic-plus-quadratic” (or weak attractive) potential obtained for and (see the Thomas-Fermi approximation in §7.1).
Regime 2: weak “quartic-minus-quadratic” (or weak repulsive) potential obtained when and ; this regime appears when .
Regime 3: strong “quartic-minus-quadratic” (or strong repulsive) potential obtained when and ; this regime appears when .

Figure 11: Computation of a rotating BEC with giant vortices (cf. §7.3). 3D-rendering of the atomic density (normalized by its maximum value ) obtained in Regime 1 (a, b, c), Regime 2 (d, e, f) and Regime 3 (g, h, i) for different rotations (first column), (second column) and (third column).

All computations are performed with the (RGC) method and the obtained stationary states are presented in Figure 11. The parameters for these simulations are: , and (Regime 1), (Regime 2), (Regime 3). In the first column of Figure 11 we notice that the atomic density distribution in the condensate without rotation () changes from the classical parabolic profile in Regime 1 to a Mexican-hat type profile in Regime 2 and, finally, to a profile with a central hole in Regime 3. It is then expected that, when rotation is applied, in Regimes 2 and 3 the condensate will develop a central hole (or giant vortex) at lower rotations frequencies than in Regime 1. This prediction is indeed supported by the results in Figure 11 (second and third column). When rotation is increased, the condensate configuration evolves from a classical (Abrikosov) vortex lattice to a vortex lattice with a central depletion and, finally, to a giant vortex surrounded by a ring of individual vortices. The giant vortex is indeed obtained for lowest rotation frequencies in Regimes 2 and 3. The existence of a giant vortex was predicted theoretically (e. g. [30]) and was already observed in 2D [36] and 3D computations [5, 25].

8 Conclusions

The difficulty of direct minimization of the Gross-Pitaevskii energy functional with rotation comes from the unit-norm constraint (2). The novel idea proposed here is to transform this problem to an unconstrained Riemannian optimization problem defined on a spherical manifold and then develop a Riemannian conjugate gradient (RCG) method based on classical approaches. The key enablers of this new method are the following: (i) the gradient direction is derived using the theory of Sobolev gradients and relies on a physically-inspired definition of the inner product which accounts for rotation [27], thereby offering a good preconditioning for the problem; (ii) the gradient is projected on the subspace tangent to the spherical manifold before being used in simple gradient or conjugate gradient methods, which ensures the iterates stay close to manifold ; (iii) the conjugate descent direction is computed using classical approaches (i. e. the Polak-Ribière or Fletcher-Reeves variant of the nonlinear conjugate gradient method) and the Riemannian vector transport is used to bring the gradient and descent directions determined at the previous iteration to the current tangent subspace ; (iv) the optimal descent step is computed by solving an arc-minimization problem (in which samples are constrained to lie on the manifold), instead of the classical line-minimization; (v) finally, the updated solution is “retracted” back to the spherical manifold.

As demonstrated by our tests performed in the finite-element setting, several features make the (RCG) method very appealing for practical computations: (i) since the “optimize-then-discretize” paradigm is used, the preconditioning is mesh-independent; (ii) the Riemannian retraction and transport operators are simple to implement; (iii) for the arc-minimization problem a classical approach such as Brent’s method can be easily adapted; (iv) there are no tuning parameters or trust-region tests involved. In addition, general mesh refinement or mesh adaptivity strategies are compatible with the RCG method without any modifications. Our extensive numerical experiments showed a significant improvement of the convergence rate of the RCG method over the simple gradient and imaginary-time methods. Finally, as a challenging test, the (RCG) method was used to compute vortex configurations in rotating BEC with high values of the nonlinear interaction constants and very high rotation rates.

9 Acknowledgements

I. D. acknowledges the support through the ANR (France) grant ANR-12-MONU-0007-01 BECASIM (call “Modéles Numériques”). B. P. acknowledges the support through an NSERC (Canada) Discovery Grant. Computational resources were provided by CRIANN (Centre Régional Informatique et d’Applications Numériques de Normandie, France) under the project 2015001. The authors acknowledge the generous hospitality of the Fields Institute in Toronto during the Thematic Program on Multiscale Scientific Computing (January–April, 2016).


  • [1] P.-A. Absil, R. Mahony, and R. Sepulchre. Optimization Algorithms on Matrix Manifolds. Princeton University Press, 2008.
  • [2] R. A. Adams and J. F. Fournier. Sobolev Spaces. Elsevier, 2005.
  • [3] A. Aftalion. Vortices in Bose-Einstein Condensates. Birkhauser, 2006.
  • [4] A. Aftalion and I. Danaila. Three-dimensional vortex configurations in a rotating Bose-Einstein condensate. Physical Review A, 68:023603(1–6), 2003.
  • [5] A. Aftalion and I. Danaila. Giant vortices in combined harmonic and quartic traps. Physical Review A, 69:033608(1–6), 2004.
  • [6] A. Aftalion and Q. Du. Vortices in a rotating Bose-Einstein condensate: critical angular velocities and energy diagrams in the Thomas-Fermi regime. Phys. Rev. A, 64:063603, 2001.
  • [7] F. Alouges. A new algorithm for computing liquid crystal stable configurations: The harmonic mapping case. SIAM Journal on Numerical Analysis, 34(5):1708–1726, 1997.
  • [8] F. Alouges and C. Audouze. Gradient flows, and application to the Hartree-Fock functional. Numer. Meth. for PDE, 25:380–400, 2009.
  • [9] X. Antoine, C. Besse, and W. Bao. Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations. Comput. Phys. Comm., 184(12):2621–2633, 2013.
  • [10] X. Antoine, C. Besse, R. Duboscq, and V. Rispoli. Acceleration of the imaginary time method for spectrally computing the stationary states of Gross-Pitaevskii equations. hal-01356227, 2016.
  • [11] X. Antoine and R. Duboscq. Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates. J. Comput. Physics, 258(0):509–523, 2014.
  • [12] X. Antoine, A. Levitt, and Q. Tang. Efficient spectral computation of the stationary states of rotating Bose-Einstein condensates by the preconditioned nonlinear conjugate gradient method. arxiv:1611.02045, 2016.
  • [13] W. Bao. Ground states and dynamics of rotating Bose-Einstein condensates, pages 215–255. Modeling and Simulation in Science, Engineering and Technology. Birkhauser, 2006.
  • [14] W. Bao. Mathematical models and numerical methods for Bose-Einstein condensation. Proc. of the International Congress of Mathematicians (Seoul 2014), IV:971–996, 2014.
  • [15] W. Bao and Y. Cai. Mathematical theory and numerical methods for Bose-Einstein condensation. Kinetic and related models, 6:1–135, 2013.
  • [16] W. Bao and Q. Du. Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow. Siam J. Sci. Comput., 25:1674, 2004.
  • [17] W. Bao and J. Shen. A generalized Laguerre-Hermite pseudospectral method for computing symmetric and central vortex states in Bose-Einstein condensates. J. Comput. Physics, 227:9778–9793, 2008.
  • [18] W. Bao and W. Tang. Ground-state solution of Bose-Einstein condensate by directly minimizing the energy functional. J. Comp. Physics, 187:230–254, 2003.
  • [19] N. Boumal, B. Mishra, P.-A. Absil, and R. Sepulchre. Manopt, a Matlab toolbox for optimization on manifolds. J. of Machine Learning Research, 15:1455–1459, 2014.
  • [20] V. Bretin, S. Stock, Y. Seurin, and J. Dalibard. Fast rotation of a Bose-Einstein condensate. Phys. Rev. Lett., 92:050403, 2004.
  • [21] M. Caliari and S. Rainer. GSGPEs: A Matlab code for computing the ground state of systems of Gross-Pitaevskii equations. Comput. Phys. Comm., 184(3):812 – 823, 2013.
  • [22] R.M. Caplan. NLSEmagic: Nonlinear Schrödinger equation multi-dimensional Matlab-based GPU-accelerated integrators using compact high-order schemes. Comput. Phys. Comm., 184(4):1250–1271, 2013.
  • [23] S.-L. Chang, C.-S. Chien, and B.-W. Jeng. Computing wave functions of nonlinear Schrödinger equations: A time-independent approach. J. Comput. Physics, 226(1):104 – 130, 2007.
  • [24] I. Coddington, P. C. Haljan, P. Engels, V. Schweikhard, S. Tung, and E. A. Cornell. Experimental studies of equilibrium vortex properties in a Bose-condensed gas. Phys. Rev. A, 70:063607, 2004.
  • [25] I. Danaila. Three-dimensional vortex structure of a fast rotating Bose-Einstein condensate with harmonic-plus-quartic confinement. Phys. Review A, 72:013605(1–6), 2005.
  • [26] I. Danaila and F. Hecht. A finite element method with mesh adaptivity for computing vortex states in fast-rotating Bose-Einstein condensates. J. Comput. Physics, 229:6946–6960, 2010.
  • [27] I. Danaila and P. Kazemi. A new Sobolev gradient method for direct minimization of the Gross–Pitaevskii energy with rotation. SIAM J. Sci. Computing, 32:2447–2467, 2010.
  • [28] C. M. Dion and E. Cancès. Ground state of the time-independent Gross-Pitaevskii equation. Comput. Phys. Comm., 177:787–798, 2007.
  • [29] I. Faragò and J. Karàtson. Numerical solution of nonlinear elliptic problems via preconditioning operators: Theory and applications. NOVA Science Publishers, New York, 2002.
  • [30] Alexander L. Fetter. Rotating vortex lattice in a Bose-Einstein condensate trapped in combined quadratic and quartic radial potentials. Phys. Rev. A, 64:063608, 2001.
  • [31] J. J. García-Ripoll and V. M. Pérez-García. Optimizing Schrödinger functionals using Sobolev gradients: application to quantum mechanics and nonlinear optics. SIAM J. Sci. Comput., 23:1315–1333, 2001.
  • [32] M. D. Gunzburger. Perspectives in Flow Control and Optimization. SIAM, 2003.
  • [33] F. Hecht. New developments in Freefem++. Journal of Numerical Mathematics, 20:251–266, 2012.
  • [34] F. Hecht, O. Pironneau, A. Le Hyaric, and K. Ohtsuke. FreeFem++ (manual)., 2007.
  • [35] U. Hohenester. OCTBEC a Matlab toolbox for optimal quantum control of Bose-Einstein condensates. Comput. Phys. Comm., 185(1):194–216, 2014.
  • [36] K. Kasamatsu, M. Tsubota, and M. Ueda. Giant hole and circular superflow in a fast rotating Bose-Einstein condensate. Phys. Rev. A, 66:053606, 2002.
  • [37] P. Kazemi and M. Eckart. Minimizing the Gross-Pitaevskii energy functional with the Sobolev gradient– analytical and numerical results. Int. J. of Comput. Methods, 7, 2009.
  • [38] E. H. Lieb and R. Seiringer. Derivation of the Gross-Pitaevskii equation for rotating Bose gases. Comm. in Math. Physics, 264:505–537, 2006.
  • [39] D. Luenberger. Optimization by Vector Space Methods. John Wiley and Sons, 1969.
  • [40] B. Merlet and T. N. Nguyen. Convergence to equilibrium for discretizations of gradient-like flows on Riemannian manifolds. Differential Integral Equations, 26:571–602, 2013.
  • [41] J. W. Neuberger. Sobolev Gradients and Differential Equations. Lecture notes in mathematics, 1670. Springer, Berlin/Heidelberg, 1997 (2nd Edition 2010).
  • [42] J. Nocedal and S. Wright. Numerical Optimization. Springer, 2002.
  • [43] M. Pierre. Newton and conjugate gradient for harmonic maps from the disc into the sphere. ESAIM Control Optim. Calc. Var., 10:142–167, 2004.
  • [44] L. P. Pitaevskii and S. Stringari. Bose-Einstein condensation. Clarendon Press, Oxford, 2003.
  • [45] W. H. Press, B. P. Flanner, S. A. Teukolsky, and W. T. Vetterling. Numerical Recipes: the Art of Scientific Computations. Cambridge University Press, 1986.
  • [46] P. J. Roache. Verification and Validation in Computational Science and Engineering. Hermosa Publishers, 1998.
  • [47] D. E. Sheehy and L. Radzihovsky. Vortices in spatially inhomogeneous superfluids. Phys. Rev. A, 70:063620, 2004.
  • [48] S. T. Smith. Optimization techniques on Riemannian manifolds. Fields Institute Communications, 3, 1994.
  • [49] S. Stringari. Phase diagram of quantized vortices in a trapped Bose-Einstein condensed gas. Phys. Rev. Lett., 82:4371, 1999.
  • [50] Makoto Tsubota. Quantized vortices in superfluid helium and Bose-Einstein condensates. Journal of Physics: Conference Series, 31(1):88, 2006.
  • [51] G. Vergez, I. Danaila, S. Auliac, and F. Hecht. A finite-element toolbox for the stationary Gross-Pitaevskii equation with rotation. Comput. Phys. Comm., 209:144–162, 2016.
  • [52] X. Wu, Z. Wen, and W. Bao. A regularized Newton method for computing ground states of Bose-Einstein condensates. arXiv: 1504.02891, 2015.
  • [53] R. Zeng and Y. Zhang. Efficiently computing vortex lattices in rapid rotating Bose-Einstein condensates. Comput. Phys. Comm., 180:854–860, 2009.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description