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.
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.
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.
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.
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
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
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.
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.
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 functionand 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 .
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).
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
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|