1 Introduction
We are interested in the steady state solution of the Gross–Pitaevskii equation, which describes the ground state of a quantum system of identical bosons under simplifying hypotheses, and their stability properties. Often referred to as a nonlinear Schrödinger equation, the Gross–Pitaevskii equation models certain classes of BoseEinstein condensates (BECs), a special state of matter formed by bosons at ultralow temperatures. It is well known that the solutions of the Gross–Pitaevskii equation with a parabolic trap in two dimensions exhibit a rich bifurcating behavior [16, 18, 6], which includes symmetrybreaking bifurcations and vortexbearing states [10].
The bifurcation diagram plots the number of bosons in a BEC as a function of the chemical potential. In order to trace a bifurcation diagram, one has to compute multiple solution of a parametrized, nonlinear problem. We will start with a one parameter study (the chemical potential) and then consider a two parameter case (chemical potential and the normalized trap strength). Moreover, we are interested in studying the stability property of each steady solution to the Gross–Pitaevskii equation. We assume a perturbation around the steady state characterized by a “small” perturbation amplitude and a complex eigenfrequency. At order
, we obtain an eigenvalue problem that allows us to characterize the stability of the steady solution and, in the case of an unstable state, to classify the instability. Tracing the bifurcation diagram and the numerical study of the stability of the solutions can be every expensive in terms of computational time. In this work, we propose to apply a Reduced Order Modeling (ROM) technique to reduce the demanding computational costs associated with this stability analysis.
2 The Gross–Pitaevskii equation
The Gross–Pitaevskii equation models certain classes of BoseEinstein condensates. A BoseEinstein condensate (BEC) is a special state of matter formed by an unlimited number of bosons that “condense” into the same energy state at low temperatures. A BEC is formed by cooling a gas of extremely low density, about onehundredthousandth the density of normal air, to ultralow temperatures (close to absolute zero).
A quantum system is the environment to be studied in terms of waveparticle duality (i.e., all particles exhibit a wave nature and viceversa). A quantum system involves the wavefunction and its constituents, such as the momentum and wavelength. The Gross–Pitaevskii equation describes the ground state of a quantum system of identical bosons using two simplifications: the Hartree–Fock approximation and the pseudopotential interaction model. In the Hartree–Fock approximation, the total wavefunction of a system of bosons is taken as a product of singleparticle functions :
where is the coordinate of the th boson. If the singleparticle wavefunction satisfies the Gross–Pitaevskii equation, the total wavefunction minimizes the expectation value (i.e., the probabilistic expected value of the result of an experiment) of the pseudopotential model Hamiltonian under normalization condition:
(1) 
where is the domain under consideration and is interpreted as the particle density. The Gross–Pitaevskii equation reads: Find the singleparticle wavefunction such that
(2) 
where is the imaginary unit, is the radial coordinate, and is the external potential, with being the normalized trap strength, i.e. the ratio of trappings along and transverse to the plane. Here, we choose unless specified otherwise. Notice that we consider a single well potential. Eq. (2) is similar in form to the Ginzburg–Landau equation and is sometimes referred to as a nonlinear Schrödinger equation. Obviously, eq. (2) needs to be supplemented with suitable boundary conditions.
The construction of the steady solution is based on the ansatz:
(3) 
where is the chemical potential, which has to satisfy . By plugging (3) into (2), we obtain nonlinear problem
(4) 
It is well known that the solutions of onedimensional version of problem (4) exhibit a bifurcating behavior [14, 13, 2, 8], which is not particularly rich though. The bifurcations occurring in the twodimensional problem (4) are far more interesting [16, 18, 6]. Indeed, several secondary bifurcations appear, which include symmetrybreaking bifurcations and vortexbearing states [10]. The bifurcation diagram plots the number of bosons in the BEC (1) as a function of the chemical potential . When , the nonlinearity of the problem becomes irrelevant and the states bifurcate from the respective linear limit. Since an arbitrary potential can be approximated as a harmonic potential at the vicinity of a stable equilibrium point, when
we can decompose the linear eigenfunction
in Cartesian form as being proportional towhere is the Hermite polynomial with being the associated quantum number of the harmonic oscillator. The eigenvalue corresponding to linear eigenfunction is . Notice that the characteristic eigenvalue parameter, i.e. the eigenvalue responsible for the bifurcation, is the chemical potential.
Let be a steady state solution to (4). In order to study its stability property we assume a perturbation ansatz around of the form
(5) 
where is the complex eigenfrequency, is a “small” amplitude of the perturbation, and the star stands for complex conjugation. Plugging (5) into eq. (2), at order we obtain an eigenvalue problem that can be written in matrix form as:
with eigenfrequencies
, and matrix blocks given byThe steady state is classified as stable in this Hamiltonian system if all the eigenfrequencies have vanishing imaginary part . On the other hand, when the solution becomes unstable, two types of instabilities can occur: i) exponential instabilities characterized by a pair of imaginary eigenfrequencies with zero real part , and ii) oscillatory instabilities characterized by a complex eigenfrequency quartet.
3 Numerical approximation of the problem
The Gross–Pitaevskii equation (2) describes a complex physical phenomenon with a rich bifurcating behavior, which is due to its nonlinearity. The bifurcation diagrams and the stability of the associated steady states can be studied numerically. The reconstruction of a bifurcation diagram entails computing multiple solution of a parametrized, nonlinear problem, which can be computationally very demanding. For this reason we decided to take advantage of Reduced Order Methods (ROMs).
ROMs have been introduced for parametrized problems requiring realtime capabilities due to a manyquery setting. The goal is to compute reliable results at a fraction of the cost of a conventional (full order) method. A practical way to realize this is to organize the computation in two steps:

An offline phase: full order approximation solutions corresponding to selected representative parameters values/system configurations are computed and stored, together with other information concerning the parametrized problem. This is a computationally expensive step usually performed on high performance computing facilities.

An online phase: the information obtained during the offline phase is used to compute the solution for a newly specified value of the parameters in a short amount of time (ideally in real time), even on a relatively low power device such as a laptop or a smartphone.
These split computational procedures are built in such a way that new parameter dependent quantities are easily and quickly computed online, while representative basis functions for selected parameter values and more demanding quantities are precomputed offline. Among all the possible ROMs, we choose the Reduced Basis technique
The algorithm we use to rapidly plot the bifurcation diagram, presented in [20], is based on three building blocks :

A continuation technique: it permits to reconstruct properly the bifurcation diagram by following each branch and providing a proper initial guess for the nonlinear iterations (next building block). We use a slight modification of the continuation method in [3], which consists in a for loop in the parameter domain where at each new cycle we check if the problem bifurcates. If it does bifurcate, we set that bifurcated solution as the initial guess for the next cycle.

A Galerkin Finite Element discretization: at each step of the Newton iteration, we have to approximate the solution of a new linear weak formulation. For the space discretization, we use the GalerkinFinite Element method. This choice is made also in view of the numerical extension towards the model order reduction.
Algorithm 1 summarizes the steps for a generic nonlinear parametric problem that reads: given a scalar parameter , find , being a given functional space, such that
(6) 
The next subsection is devoted to the description of building blocks 2 and 3.
3.1 Galerkin Finite Element method
As mentioned in the previous subsection, for the space discretization of the Gross–Pitaevskii equation we apply the Galerkin finite element method. Here, we introduce some standard notion for the discretization of generic problem (6) and then write the discrete counterpart of eq. (4).
In order to write the variational formulation the Galerkin Finite Element method builds upon, we introduce a family of finite dimensional spaces , such that . Let . Given , we seek that satisfies
(7) 
We can not directly apply the finite element method to the equation (7) because of the nonlinearity in . So, we first apply the NewtonKantorovich method [7, 21], which reads as follows. Let be an initial guess. For every , we seek the variation such that
(8) 
and then we update the solution for the successive step as , until we reach convergence for a prescribed tolerance .
From the algebraic point of view we denote with a basis for , such that the Newton’s method combined with the Galerkin finite element method and applied to our weak formulation reads: find such that
(9) 
where the Jacobian matrix in is defined as
(10) 
We can now present the weak formulation projected into the Finite Element space for problem (4). We recall that the solution to eq. (4) is a complex function. Let and be its real and imaginary part, respectively. The generic th iteration of the Newton’s method (9) reads: seek , such that
(11)  
where
3.2 Reduced Basis method
In the previous section we derived a parametrized weak formulation that, once projected in a suitable Finite Element space, gives us a huge nonlinear system which has to be solved for every parameter . Thus, the computational cost needed in order to study the bifurcating phenomena of the GrossPitaevskii equation has to be controlled, in such a way that we can explore rapidly and reliably all the branches. A possible strategy is the application of the reduced order methods (ROM), a collection of methodologies used to replace the original high dimension problem, typically called high fidelity approximation, with a reduced order one that is easy to manage.
One of these methodologies is the Reduced Basis method (RB), that roughly speaking consists in a projection of the high fidelity problem on a subspace of smaller dimension, constructed with some properly chosen basis functions.
As we have seen, the preliminary step is the projection of the weak formulation in a discrete space, which results in the Galerkin problem (7). Finding a numerical solution to this problem is very challenging because of the potential high number of degrees of freedom . Thus, we aim at constructing a discrete reduced order manifold, in which the parametrized solutions lie, that efficiently approximate the high fidelity one. This is the description of the first step, namely the offline phase, in which we explore the parameter space in order to construct a basis for the low dimensional manifold, solving times the Galerkin highfidelity problem varying .
On the other side, the second step, called online phase, is the efficient and reliable part where the solution is computed through the projection on the reduced manifold. In this way we do not have to project the equation in the huge Finite Element space for every , rather we can project on the low dimensional space spanned by the selected base. This complexity reduction is based on two main key points: the assumption that it holds the affine decomposition and the fact that . The latter assumption means that we can approximate the discrete manifold with a low number of basis functions.
To be precise we want to construct the reduced problem through the projection on a subspace spanned by a collection of the so called snapshots, i.e. the solution of the full order problem for different values of , obtained by Proper orthogonal decomposition (POD) or Greedy techniques [11, 19, 21].
As well as in the offline phase (7), given , we seek that satisfies
(12) 
At this point we have again to face up with the nonlinearity, and so come back to the NewtonKantorovich method obtaining: given an initial guess , for every we seek the variation such that
(13) 
and then we update the solution as until convergence.
During the offline phase, through the POD sampling, we obtain an orthonormal base given by , such that if we denote we can write every element as
(14) 
denoting with
the reduced solution vector.
Now we can choose the test element as for every , obtaining the algebraic system in given by
(15) 
Moreover, we denote with
the residual reduced vector and with the transformation matrix whose elements are the nodal evaluation of the mth basis function at the jth node. Furthermore, with this new notation we observe that (15) corresponds to the solution of
(16) 
We can now recovery the algebraic low dimensional system by combining the Newton’s method and the Reduced Basis technique. Thus the problem reads as follows: find such that
(17) 
where is the reduced Jacobian matrix defined as
(18) 
Furthermore, once finished the offline phase, we have a base for the space , given by the span of the orthonormalized (through the GramSchmidt procedure) POD modes. In order to construct properly the bifurcation plots, the training set for the POD technique, with the values for which we solved the offline problem, was setted as an ordered sampling of the interval .
Before ending this section we want to analyze in a more detailed way the system we obtained. Now we show how the Reduced Basis method reflects the projection properties of the Galerkin method also in the online phase. In fact, we can project the weak formulation (7) in the reduced manifold , obtaining the following: given and an initial guess for until convergence we seek such that
(19)  
and then set . The matrix formulation follows directly and we obtain exactly the equation (17), where the reduced Jacobian is given by
(20) 
where we have defined the reduced matrices as:
in terms of the high fidelity ones coming from (17).
We have illustrated the online phase, that permits an efficient evaluation of the solution and possibly related outputs for every possible choice of a different parameter . The key point of this time saving is the so called affine decomposition. Indeed we want the computations to be independent form the, usually very high, number of degrees of freedom of the true discrete problem. In general the reduced matrices we have just presented are
dependent and an affinerecovery technique called Empirical Interpolation Method (EIM) is needed
[5].4 Results
All the results reported in this section have been obtained with FEniCS [1, 15, 4] for the high fidelity phase, and RBniCS [22] for the reduced order one.
4.1 Validation of the full order method
To validate the Full Order Method described in Sec. 3.1, we use a test proposed in [6]. We approximate the solution to eq. (4) in domain . Homogeneous Dirichlet boundary conditions are imposed on the entire boundary of . We recall that we set . For the space discretization, we use finite elements and a mesh with 6889 elements.
Fig. 1 displays the bifurcation diagram in the plane (left) and in the plane (right) obtained with the Full Order Method. These diagrams show the first three bifurcation points and the relative nonuniqueness of the solution with respect to the parameter . As is increased, the sequence of events observed in Fig. 1 is as follows. The ground state is the system simplest state. Its linear eigenfunction has corresponding eigenvalue . The ground state is generically stable, thus no further bifurcations occur from this state [12]. As expected, a unique solution branch departs from in Fig. 1. A representative density function for this branch is shown in Fig. 1(a). We see no further bifurcation for . The first interesting events in terms of bifurcation analysis occur for , with : two branches, associated to and , bifurcate from point in the and planes [17, 9]. Indeed, from point in Fig. 1 we observe the two expected branches. Representative density functions for these two branches are reported in Fig. 1(b) and 1(c). The next, more complicated, case of bifurcations emanates from point , with . In Fig. 1, we see that three branches depart from this point, associated to , , and . The corresponding representative densities are shown in Fig. 1(f), 1(e), and 1(d).
The density functions reported in Fig. 2, which are associated to and all the 6 solution branches in Fig. 1, show the richness of density patterns in order of decreasing In particular, we see the ground state in Fig. 1(a), the single charge vortex in Fig. 1(b), the 1dark soliton stripe in Fig. 1(c), the dark soliton cross in Fig. 1(f), the ring dark soliton in Fig. 1(e), the 2dark soliton stripe in Fig. 1(d). Notice that the 6 branches Fig. 1 are related to the first three eigenvalues. For example, the second bifurcation stems from a double eigenvalue and thus we have two branches, this phenomena is called multiple bifurcations. The stability property of these branches are different for each case, i.e the single charge vortex is always stable while the 1dark soliton stripe is subject to multiple secondary bifurcations. These properties can be easily studied using standard techniques and are not among the purpose of this work.
To obtain the 6 branches of bifurcation diagrams in Fig. 1, we used a simple continuation method on parameter , combined with the aforementioned Finite Element discretization technique and the Newton’s method as nonlinear solver. The overall simulation time required to complete the diagrams is roughly 95 minutes with continuation step . Consider that the mesh we used is quite coarse.
4.2 Reduced order approach
Next, we present the results obtained with our reduced order method as described in Sec. 3.2. Fig. 3 shows the first 6 bifurcation branches. We see that the curves obtained with the reduced order method matched very closely the curves computed with the full order method shown in Fig. 1. As further evidence of the accuracy of our ROM approach, we plot in Fig. 4 the difference between ROM and FOM solution for the solution and the density function in the and norms for each branch in Fig. 3.
As we said the result are exactly the same, but we gain a speedup of almost 1.1, so to construct this plots we took 5147 seconds. In fact as we can see in the next figures for the different branch, to obtain this computational advantage we paid a negligible worsening of the accuracy.
Now we can observe that we saved almost 600 seconds, and still our computations are not independent from the number of degrees of freedom of the Finite Element Method. To improve a lot the results in terms of speedup we need to implement some affine recovery technique, such as EIM or DEIM.
So we can present again the bifurcation diagrams and the representative solutions of each branch.
We investigated also the 2 parameters case, considering the trap strength as a parameter varying in . In this case so we want to detect the evolution of the first bifurcations for different trapping configuration.
4.3 Hyperreduction techniques
We used EIM/DEIM for the construction of one of the branch of the second bifurcation and we obtained results comparable with the previous ones as we can see in the next figures.
Moreover, the error plots the results are exactly the same in terms of accuracy with respect to the affine recovery techniques.
In both this cases we obtained a good performance in terms of reduced error plot, but also we improved a lot the computational savings. In fact here the online phase is independent from the number of degrees of freedom of the finite element method, so the RB method is really efficient.
For the EIM we took just 55 seconds to construct the reduced basis bifurcation diagrams with respect to the 246 seconds for the high fidelity one, so obtaining a speedup of almost 5. Moreover we obtained a surprisingly better result using DEIM. In fact it tooks just 7 seconds for the complete construction of the reduced basis bifurcation diagrams, with respect to the 246 seconds for the high fidelity one, so obtaining a speedup of almost 32.
The hyperreduction techniques are crucial in the two parameters test case. In fact, in order to reconstruct the diagram in Fig. 7 we need 14708 seconds while applying the DEIM technique we just need 366 seconds, with a speedup of more then 40 times.
Acknowledgments
This work was supported by European Union Funding for Research and Innovation through the European Research Council (project H2020 ERC CoG 2015 AROMACFD project 681447, P.I. Prof. G. Rozza). This work was also partially supported by NSF through grant DMS1620384 (A. Quaini). We acknowledge fruitful conversations with Dr. Sean Garner.
References
 [1] FEniCS Project. https://fenicsproject.org.
 [2] G L Alfimov and D A Zezyulin. Nonlinear modes for the gross–pitaevskii equation—a demonstrative computation approach. Nonlinearity, 20(9):2075, 2007.
 [3] E. Allgower and K. Georg. Introduction to Numerical Continuation Methods. Society for Industrial and Applied Mathematics, 2003.
 [4] Martin S. Alnæs, Jan Blechta, Johan Hake, August Johansson, Benjamin Kehlet, Anders Logg, Chris Richardson, Johannes Ring, Marie E. Rognes, and Garth N. Wells. The fenics project version 1.5. Archive of Numerical Software, 3(100), 2015.

[5]
M. Barrault, N. C. Nguyen, Y. Maday, and A. T. Patera.
An “empirical interpolation” method: Application to efficient reducedbasis discretization of partial differential equations.
C. R. Acad. Sci. Paris, Série I., 339:667–672, 2004.  [6] E.G. Charalampidis, P.G. Kevrekidis, and P.E. Farrell. Computing stationary solutions of the twodimensional gross–pitaevskii equation with deflated continuation. Communications in Nonlinear Science and Numerical Simulation, 54:482 – 499, 2018.
 [7] P. G. Ciarlet. Linear and Nonlinear Functional Analysis with Applications:. Other Titles in Applied Mathematics. Society for Industrial and Applied Mathematics, 2013.
 [8] M P Coles, D E Pelinovsky, and P G Kevrekidis. Excited states in the large density limit: a variational approach. Nonlinearity, 23(8):1753, 2010.
 [9] Andres Contreras and Carlos GarciaAzpeitia. Global bifurcation of vortex and dipole solutions in boseeinstein condensates. Comptes Rendus Mathematique, 354(3):265 – 269, 2016.
 [10] C. GarcíaAzpeitia and D. E. Pelinovsky. Bifurcations of multivortex configurations in rotating bose–einstein condensates. Milan Journal of Mathematics, 85(2):331–367, Dec 2017.
 [11] J. Hesthaven, G. Rozza, and B. Stamm. Certified Reduced Basis Methods for Parametrized Partial Differential Equations. Springer Briefs in Mathematics, 2015.
 [12] P. Kevrekidis, D. Frantzeskakis, and R. CarreteroGonzález. The Defocusing Nonlinear Schrödinger Equation. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2015.
 [13] P G Kevrekidis, V V Konotop, A Rodrigues, and D J Frantzeskakis. Dynamic generation of matter solitons from linear states via timedependent scattering lengths. Journal of Physics B: Atomic, Molecular and Optical Physics, 38(8):1173, 2005.
 [14] Yuri S. Kivshar, Tristram J. Alexander, and Sergey K. Turitsyn. Nonlinear modes of a macroscopic quantum oscillator. Physics Letters A, 278(4):225 – 230, 2001.
 [15] Anders Logg, KentAndre Mardal, Garth N. Wells, et al. Automated Solution of Differential Equations by the Finite Element Method. Springer, 2012.
 [16] S. Middelkamp, P. G. Kevrekidis, D. J. Frantzeskakis, R. CarreteroGonzález, and P. Schmelcher. Bifurcations, stability, and dynamics of multiple matterwave vortex states. Phys. Rev. A, 82:013646, Jul 2010.
 [17] S. Middelkamp, P. G. Kevrekidis, D. J. Frantzeskakis, R. CarreteroGonzález, and P. Schmelcher. Bifurcations, stability, and dynamics of multiple matterwave vortex states. Phys. Rev. A, 82:013646, Jul 2010.
 [18] S. Middelkamp, P.G. Kevrekidis, D.J. Frantzeskakis, R. CarreteroGonzález, and P. Schmelcher. Emergence and stability of vortex clusters in boseeinstein condensates: A bifurcation approach near the linear limit. Physica D: Nonlinear Phenomena, 240(18):1449 – 1459, 2011.

[19]
A.T. Patera and G. Rozza.
Reduced basis approximation and A posteriori error estimation for Parametrized Partial Differential Equation
. MIT Pappalardo Monographs in Mechanical Engineering, Copyright MIT (20072010). http://augustine.mit.edu.  [20] F. Pichi and G. Rozza. Reduced basis approaches for parametrized bifurcation problems held by nonlinear Von Kármán equations. Journal of Scientific Computing, 2019. doi: 10.1007/s10915019010033.
 [21] A. Quarteroni, A. Manzoni, and F. Negri. Reduced Basis Methods for Partial Differential Equations: An Introduction. UNITEXT. Springer International Publishing, 2015.
 [22] RBniCS. http://mathlab.sissa.it/rbnics.
Comments
There are no comments yet.