Efficient time integration methods for Gross–Pitaevskii equations with rotation term

10/26/2019 ∙ by Philipp Bader, et al. ∙ 0

The objective of this work is the introduction and investigation of favourable time integration methods for the Gross–Pitaevskii equation with rotation term. Employing a reformulation in rotating Lagrangian coordinates, the equation takes the form of a nonlinear Schrödinger equation involving a space-time-dependent potential. A natural approach that combines commutator-free quasi-Magnus exponential integrators with operator splitting methods and Fourier spectral space discretisations is proposed. Furthermore, the special structure of the Hamilton operator permits the design of specifically tailored schemes. Numerical experiments confirm the good performance of the resulting exponential integrators.



There are no comments yet.


page 19

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction


Our main concern is the introduction and investigation of efficient numerical methods for nonlinear evolution equations involving an explicit time-dependency


previous work on non-autonomous linear evolution equations


has enhanced our interest in this subject. In the present work, we center our attention on time-dependent Gross–Pitaevskii equations with rotation term modelling the creation of vortices in Bose–Einstein condensates.

Exponential integrators

Theoretical and numerical studies of time integration methods for systems of non-autonomous linear ordinary differential equations that can be cast into the form (

2) with matrix-valued time-dependent function are found in [Alvermann, Fehske (2011), Alvermann, Fehske, Littlewood(2012), Bader, Blanes, Kopylov (2018), Blanes, Casas, González, Thalhammer (2018), Blanes, Casas, Thalhammer (2017), Blanes, Casas, Thalhammer (2018)], see also references given therein. Numerical comparisons presented within the context of large-scale applications arising in quantum physics clearly indicate that exponential integrators based on the Magnus expansion are favourable compared to standard integrators such as explicit Runge–Kutta methods; moreover, it is demonstrated that an approach which avoids the actual computation of matrix commutators

by considering compositions of matrix exponentials leads to approximations that are superior to those obtained by Magnus integrators. In view of an appropriate extension to partial differential equations and finally to the nonlinear case (

1), it is helpful to understand this approach as a replacement of the exact evolution operator associated with (2) by a sequence of autonomous evolution equations, which are essentially of the form


for certain fixed times ; we refer to this type of time integration methods as commutator-free quasi-Magnus (CFQM) exponential integrators.

Computational aspects

For non-autonomous linear ordinary differential equations or spatial semi-discretisations of non-autonomous linear partial differential equations leading to ordinary differential equations of the form (2), respectively, the realisation of CFQM exponential integrators relies on the numerical computation of matrix exponentials. For non-autonomous linear Schrödinger equations of a special structure

written as abstract evolution equations, an alternative methodology is to employ operator splitting [Blanes, Moan (2002), McLachlan, Quispel (2002)] combined with spectral methods; for instance, in situations where the Hamiltonian comprises the Laplace operator with respect to the spatial variables and a real-valued space-time-dependent potential that satisfies appropriate regularity and periodicity conditions

favourable approximations to (3) are based on suitable compositions of the solutions to the associated autonomous subequations

computed by Fourier spectral space discretisation and pointwise multiplication, respectively.

Modified integrators

For Schrödinger equations defined by a Hamilton operator of the special structure (4), the incorporation of matrix commutators corresponds to modifications of the potential; indeed, straightforward calculations show that certain Lie-brackets of the second space derivative and a multiplication operator are given by


Together with Taylor series expansions and polynomial interpolation, this permits to design higher-order schemes with improved efficiency, see 

[Bader, Blanes, Kopylov (2018)]; we refer to this type of methods as modified CFQM exponential integrators, and with regard to the initials of the authors, to a specific sixth-order scheme as BBK scheme. In combination with operator splitting, this approach leads to autonomous subequations of the form

where  represents certain linear combinations of the values of the potential and  involves first space derivatives of the potential; again, their numerical solution relies on Fourier spectral space discretisation and pointwise multiplication, respectively.

Extension to nonlinear equations

In this work, we provide an expedient extension of CFQM exponential integrators and modified CFQM exponential integrators for non-autonomous linear evolution equations (2) to the substantially more complex nonlinear case (1). With regard to relevant quantum physical applications in the field of Bose–Einstein condensation, modelled by time-dependent Gross–Pitaevskii equations with rotation term, we focus on non-autonomous nonlinear evolution equations of Schrödinger type that can be cast into the form


comparably to (4), represents the Laplace operator, a real-valued space-time-dependent trapping potential, and  a real-valued nonlinear multiplication operator. The restriction to this particular setting permits the implementation of CFQM exponential integrators as well as modified CFQM exponential integrators by operator splitting and Fourier spectral methods.

Strategy for full discretisation

We point out that the extension of the afore mentioned concepts to the nonlinear case is not entirely straightforward, since different approaches seem likely and the overall computational cost has to be taken into account.

  1. A strategy for the discretisation of (6), which is in principle practicable, uses operator splitting, that is, a decomposition into linear and nonlinear terms

    The numerical solution of the non-autonomous linear subequation relies on CFQM exponential integrators and once again on operator splitting methods; due to an invariance property that is characteristic for nonlinear Schrödinger equations, the nonlinear subequation reduces to a linear equation

    resolvable by pointwise multiplication. For a splitting method comprising  compositions and a CFQM exponential integrator involving  compositions, this strategy requires the numerical solution of linear subequations by Fourier spectral space discretisation as well as linear subequations by pointwise multiplication, since

  2. In view of CFQM exponential integrators and splitting methods involving a higher number of compositions, potentially yielding higher accuracy and efficiency, it is advantageous to employ a different strategy. At first, a CFQM exponential integrator is applied, which leads to a sequence of autonomous nonlinear evolutions equations

    in accordance with (3); combined with operator splitting, this alternative approach requires the numerical solution of linear subequations of the form

    for certain fixed values by Fourier spectral space discretisation and pointwise multiplication, respectively, each  times, since

Magnus integrators

It is also noteworthy that a formal extension of standard Magnus integrators would lead to involved and impractical numerical approximations. More precisely, rewriting the non-autonomous nonlinear differential equation (1) by means of the formal calculus of Lie-derivatives as non-autonomous linear differential equation

the associated solution representation based on the Magnus expansion

involves nested integrals of Lie-brackets. Specifically, for non-autonomous nonlinear Schrödinger equations (6), the Lie-bracket  corresponds to a second-order differential operator with coefficients depending on the solution values.

High accuracy

A characteristic of rotating Bose–Einstein condensates, modelled by time-dependent Gross–Pitaevskii equations with rotation term, is the creation of vortices; in order to simulate this quantum physical phenomenon correctly, high spatial and temporal resolution is required, and thus it is of relevance to employ accurate and efficient full discretisation methods. Given the relevance of the theme, numerous contributions are concerned with the development of advanced methodologies and illustrate their benefits over standard numerical methods; we in particular mention an approach based on generalised-Laguerre–Fourier–Hermite spectral space discretisation and operator splitting methods that has been studied in [Bao, Li, Shen (2009), Hofstätter, Koch, Thalhammer (2014)] and the transformation to the rotating frame combined with the application of exponential integrators as proposed in [AntoineBesseBao (2013), Bader(2014), Bader (2015), Bao, Marahrens, Tang, Zhang (2013), Besse, Dujardin, Lacroix-Violet (2017)]. We point out that a rigorous stability and local error analysis has shown the reliability of generalised-Laguerre–Fourier–Hermite pseudo-spectral time-splitting methods, see [Hofstätter, Koch, Thalhammer (2014)]; but, due to the fact that the implementation of these specialised spectral methods is involved and computationally expensive, we believe that it is expedient to develop [Bao, Marahrens, Tang, Zhang (2013)] further. In order to design higher-order schemes with improved accuracy and efficiency, we transform the considered time-dependent Gross–Pitaevskii equation to rotating Lagrangian coordinates; this yields a non-autonomous nonlinear Schrödinger equation of the form (6), where the afore sketched strategy can be realised by fast Fourier transforms.


The use of standard notation implicates that the same letter denotes diverse quantities in different contexts. For instance, in Sections 2 and 5, the real number denotes the constant angular velocity; by contrast, in Section (3), the function denotes the exponent in the Magnus expansion.

2 Rotational Gross–Pitaevskii equation

Rotating Bose–Einstein condensates

The phenomenon of Bose–Einstein condensation, for the first time observed in 1995 [Anderson et al. (1995), Bradley et al. (1995), Davis et al. (1995)], has opened new perspectives on quantum physics. Over the last two decades, rotating Bose–Einstein condensates with the creation of quantised vortices as a particular indication of superfluidity [Abo et al. (2001), Madison et al. (2001)] have attracted a lot of interest; various theoretical and experimental studies make a contribution to a deeper understanding of the observed quantum physical effects.

Rotational Gross–Pitaevskii equation

At temperatures significantly below the critical temperature, a prevalent model for the dynamical behaviour of a rotating Bose–Einstein condensate is the time-dependent Gross–Pitaevskii equation with additional angular momentum rotation term. We study the following normalised formulation in three space dimensions (), imposing asymptotic boundary conditions and an initial condition

in view of a numerical illustration for the creation of vortices, it is convenient to consider the reduced two-dimensional model, which is included for . In (7a), denotes the complex-valued macroscopic wave function and the Laplace operator with respect to cartesian coordinates . In general, the angular momentum rotation term has the form

with direction of rotation given by the time derivative of a vector-valued function

and angular momentum operator defined by the cross product of position and gradient . Commonly, it is assumed that the axis of rotation coincides with the -axis; thus, it suffices to consider a scalar time-dependent function  and set
Further decisive quantitites in (7a) are the external potential describing the confining trap and the nonlinear function characterising the particle interaction; with regard to numerical simulations, we focus on the consideration of a quadratic potential involving positive weights
and a cubic nonlinearity
describing the strength of short-range two-body interactions in the condensate.


Our approach and the proposed exponential time integration methods are directly applicable to the case of a real-valued space-time-dependent potential and a nonlinear term of the form . The formulation (7) in particular complies with [Bao, Marahrens, Tang, Zhang (2013), Zeng, Zhang (2009)] for the special case of constant angular velocity


for simplicity, we do not include additional integral terms describing long-range dipole-dipole interactions.

Reformulation in rotating Lagrangian coordinates

In view of the introduction of space and time discretisation methods for the rotational Gross–Pitaevskii equation (7), we apply a transformation to the rotating frame

and attain the non-autonomous nonlinear Schrödinger equation

detailed calculations based on the corresponding classical Hamiltonian are given in Appendix A. In case of a space-time-dependent potential , the transformed equation instead involves .

3 Quasi-Magnus exponential integrators

In this section, we introduce CFQM exponential integrators [Alvermann, Fehske (2011), Blanes, Casas, Thalhammer (2017)] and modified CFQM exponential integrators [Bader, Blanes, Kopylov (2018)] for non-autonomous linear Schrödinger equations of the form (4); as indicated in Section 1, we favour operator splitting and Fourier spectral methods for their numerical realisation. It suffices to specify the numerical approximation for the first time step of length .

Stepwise generalisation

We first review fundamentals on interpolating functions and associated solution representations based on the Magnus expansion, since these basic principles are practical in the design of optimised high-order schemes; in this context, we restrict ourselves to the consideration of a system of non-autonomous linear ordinary differential equations

involving a constant square complex matrix and a matrix-valued time-dependent function . In view of (4), it is justified to assume that the values of  commute

whereas the matrix commutator in general does not vanish; for instance, diagonal matrices fulfill condition (9b). The formal extension of CFQM exponential integrators for non-autonomous linear ordinary differential equations to non-autonomous linear evolution equations is straightforward; however, the generalisation of a rigorous convergence analysis for differential equations defined by matrices to differential equations defined by unbounded operators is more involved, see [Blanes, Casas, González, Thalhammer (2018)]. Additional considerations lead to modified CFQM exponential integrators for non-autonomous linear Schrödinger equations.


With regard to the introduction of schemes up to order six, which we expect to be favourable in efficiency in connection with non-autonomous linear and nonlinear Schrödinger equations, we focus on the sixth-order Gauss–Legendre quadrature rule with nodes

The Lagrange polynomials through these nodes are defined by
by construction, the following identities are fulfilled
We consider the associated interpolating function
which indeed satisfies the relation
Substituting , a straightforward calculation yields the representation
consequently, Taylor series expansions
and the restriction to a suitably chosen interval verify the relations
For the following considerations, it is useful to employ a short notation for matrix commutators of the decisive quantities
with regard to the design of modified CFQM exponential integrators, we note that condition (9b) implies , whereas in general does not commute with  or , respectively. Accordingly, nested matrix commutators are denoted by
As shown in [Blanes (2018)], the order of the underlying quadrature rule is inherited by the associated differential equation. More precisely, the solutions to the differential equation (9) and the related initial value problem
coincide up to the desired order
hence, it suffices to construct suitable approximations to . For a non-autonomous linear differential equation, a formal solution representation by the matrix exponential based on the Magnus expansion [Magnus (1954)] is valid. In the present situation, the first terms in this series expansion are given by
due to the simple structure of the quadratic function , the arising integrals can be computed analytically and be expressed in terms of as well as nested commutators, see also [Blanes, Casas, Oteo, Ros (2009), Iserles, Nørsett (1999), Munthe-Kaas, Owren (1999)]. As an illustration, we specify the leading contributions; substituting and inserting the representation (10b) yields
and a similar calculation verifies
Including terms up to order six
this procedure provides an appropriate approximation to  and hence to , that is

because of the symmetry of the Magnus expansion and the quadrature rule, only odd terms appear in

. The relations in (10c) are the starting point for the design of CFQM exponential integrators up to order six; we recall that the additional condition is only used for modified CFQM exponential integrators.

CFQM exponential integrators

A th-order CFQM exponential integrator for a non-autonomous linear ordinary differential equation of the form (9) is given by a composition of several matrix exponentials that comprise certain linear combinations of the values of the defining time-dependent matrix at specified nodes


for instance, for fourth-order or sixth-order approximations based on three Gauss–Legendre quadrature nodes (10a) and positive integers , associated coefficients for are determined such that the corresponding approximate solutions  agree with  up to terms of order four or six, respectively. Generally, these coefficients are obtained from a set of order conditions, which imply


a practicable approach for the derivation of non-redundant conditions is described in [Blanes, Casas, Thalhammer (2017)]. We note that a higher number of exponentials potentially offers the possibility to design optimised schemes. The extension of (11) to non-autonomous linear Schrödinger equations

and in particular to (4) relies on the composition of the solutions to the autonomous linear Schrödinger equations

and yields a numerical approximation through

In our numerical experiments, we apply standard and optimised CFQM exponential integrators of orders ; for details, we refer to Section 5.

Modified CFQM exponential integrators

As indicated in Section 1, for non-autonomous linear Schrödinger equations with Hamilton operator given by the Laplacian and a real-valued space-time-dependent potential, more efficient modified CFQM exponential integrators can be designed from a reduced set of order conditions; more precisely, it is used that certain Lie-brackets of the -dimensional Laplacian and a space-time-dependent multiplication operator simplify as follows


see (4) and (5). In order to sketch the approach, we reconsider the non-autonomous linear ordinary differential equation (9) and in particular assume that condition (9b) holds; a rigorous derivation in the context of partial differential equations will be the subject of future work. In the present situation, the matrix commutator  vanishes; as a consequence, relation (10c) and hence the consistency condition (12) reduce to


Moreover, due to the fact that , , and the nested commutator  are of similar type, it is advantageous to incorporate terms of the form

with appropriately chosen coefficients into the composition (11).

Sixth-order scheme

A particularly efficient sixth-order modified CFQM exponential integrator has been proposed in [Bader, Blanes, Kopylov (2018)]; for a non-autonomous linear ordinary differential equation (9), it is given by

Recalling (10a) and (10b), a straightforward calculation shows that this composition of matrix exponentials is equivalent to

the decisive quantities are defined by

Accordingly to (13), it is convenient to reformulate the extension to non-autonomous Schrödinger equations (4) as
here, the matrix commutator is replaced by the Lie-bracket specified in (14). Numerical experiments indicate that this procedure yields a sixth-order approximation

a rigorous error analysis remains an open question.

4 Application to rotational Gross-Pitaevskii equation

In this section, we describe our methodology for the space and time discretisation of the Gross–Pitaevskii equation with rotation term, see (7) and (8). We indicate how to adapt the basic principles introduced in Section 3 for non-autonomous linear Schrödinger equations (4) to the more complex nonlinear case; in particular, we generalise CFQM exponential integrators (13) as well as the sixth-order modified CFQM exponential integrator (16). We sketch the implementation of CFQM exponential integrators based on operator splitting methods, see also Section 1; the realisation of the sixth-order modified CFQM exponential integrator is in the same line. A rigorous justification and convergence analysis of CFQM exponential integrators for non-autonomous nonlinear Schrödinger equations is out of the scope of the present contribution and will be the subject of future work; however, numerical illustrations presented in Section 5 confirm the appropriateness of our approach.

Abstract evolution equation

In order to identify similarities to non-autonomous linear ordinary differential equations (9) and non-autonomous linear Schrödinger equations (4), it is useful to consider the rotational Gross–Pitaevskii equation (8) on a subinterval of length and to reformulate it as abstract evolution equation

essentially, the values of the solution and the defining operators correspond to

We recall that arguments detailed in Section 1 justify the consideration of a time-independent linear part  and a time-dependent nonlinear part ; a decomposition into three operators as in (6) is no longer needed, see paragraph Strategy for full discretisation (ii).


The calculus of Lie-derivatives provides powerful tools for the formal extension of the linear to the nonlinear case; however, its application in connection with unbounded operators is not straightforward. Thus, as a first step, it is helpful to consider a non-autonomous nonlinear ordinary differential equation of the form (17a) involving a constant square complex matrix and a matrix-valued time-dependent function . Moreover, it is important to observe that the space-time-dependent potential and the nonlinearity act as multiplication operators and commute; that is, an analogous relation to (9b) holds with Lie-brackets replacing matrix commutators


Suitable adaptations of the arguments given in Section 3 yield the analogues of the fundamental expansions (10c) and (15), and in succession, the analogoues of CFQM exponential integrators (13) and the modified CFQM exponential integrator (16), which we describe next.

CFQM exponential integrators

The generalisation of CFQM exponential integrators to the rotational Gross–Pitaevskii equation (17) requires the numerical solution of a sequence of autonomous Gross–Pitaevskii equations

where and , yielding an approximation to the exact solution through
see also (11) and (13); denoting
the nonlinear evolution equations arising in (19a) corresponds to

with for