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

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.

## Authors

• 3 publications
• 5 publications
• 7 publications
• 2 publications
• ### A Deuflhard-type exponential integrator Fourier pseudospectral method for the "Good" Boussinesq equation

We propose an exponential integrator Fourier pseudospectral method DEI-F...
10/08/2019 ∙ by Chunmei Su, et al. ∙ 0

• ### Spectral splitting method for nonlinear Schrödinger equations with quadratic potential

In this paper we propose a modified Lie-type spectral splitting approxim...
10/27/2021 ∙ by Andrea Sacchetti, et al. ∙ 0

• ### Efficient Numerical Algorithms for the Generalized Langevin Equation

We study the design and implementation of numerical methods to solve the...
12/08/2020 ∙ by Benedict Leimkuhler, et al. ∙ 0

• ### Applying splitting methods with complex coefficients to the numerical integration of unitary problems

We explore the applicability of splitting methods involving complex coef...
04/06/2021 ∙ by S. Blanes, et al. ∙ 0

• ### Pseudospectral computational methods for the time-dependent Dirac equation in static curved spaces

Pseudospectral numerical schemes for solving the Dirac equation in gener...
09/06/2019 ∙ by Xavier Antoine, et al. ∙ 0

• ### Maximum bound principles for a class of semilinear parabolic equations and exponential time differencing schemes

The ubiquity of semilinear parabolic equations has been illustrated in t...
05/23/2020 ∙ by Qiang Du, et al. ∙ 0

• ### Adaptive pseudo-time methods for the Poisson-Boltzmann equation with Eulerian solvent excluded surface

This work further improves the pseudo-transient approach for the Poisson...
11/29/2020 ∙ by Benjamin Jones, et al. ∙ 0

##### 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

### Scope

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

 u′(t)=F(t,u(t)),t∈(t0,T); (1)

previous work on non-autonomous linear evolution equations

 u′(t)=F(t)u(t),t∈(t0,T), (2)

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

 [F(s),F(t)]=F(s)F(t)−F(t)F(s),s,t∈[t0,T],

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

 u′(t)=F(t∗)u(t),t∈(t0,T), (3)

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

 iu′(t)=iF(t)u(t)=H(t)u(t)=(A+B(t))u(t),t∈(t0,T), (4a) 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 W:Rd×[t0,T]→R that satisfies appropriate regularity and periodicity conditions H(t)=−12Δ+W(⋅,t),A=−12Δ,B(t)=W(⋅,t), (4b)

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

 iu′(t)=Au(t),iu′(t)=B(t∗)u(t),

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

 (5)

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

 iu′(t)=Au(t),iu′(t)=ˆBu(t),iu′(t)=(ˆB+˜B)u(t),

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

 iu′(t)=iF(t,u(t))=H(t,u(t))=(A+B(t)+C(u(t)))u(t),t∈(t0,T); (6)

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

 iu′(t)=Au(t)+B(t)u(t),iu′(t)=C(u(t))u(t).

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

 iu′(t)=C(u(t∗))u(t),

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

 application of r-stage splitting method:{r times:iu′(t)=Au(t)+B(t)u(t),r times (pointwise multiplication):iu′(t)=C(u(t∗))u(t),application of s-stage CFQM exponential integrator:{rs times:iu′(t)=Au(t)+B(t∗)u(t),application of r-stage splitting method:{r2s times (Fourier spectral method):iu′(t)=Au(t),r2s times (pointwise multiplication):iu′(t)=B(t∗)u(t).
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

 iu′(t)=Au(t)+B(t∗)u(t)+C(u(t))u(t),

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

 iu′(t)=Au(t),iu′(t)=B(t∗)u(t)+C(u(s∗))u(t),

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

 application of s-stage CFQM exponential % integrator:{s times:iu′(t)=Au(t)+B(t∗)u(t)+C(u(t))u(t),application of r-stage splitting method:{rs times (Fourier spectral method):i% u′(t)=Au(t),rs times (pointwise multiplication):iu′(t)=B(t∗)u(t)+C(u(s∗))u(t).

### 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

 u′(t)=L(t)u(t),L(t)≡F(t,u(t))δδu,t∈(t0,T),

the associated solution representation based on the Magnus expansion

 u(t)=eΩ(t−t0)u0,Ω(t−t0)=∫tt0L(τ)dτ+12∫tt0∫τt0[L(τ),L(σ)]dσdτ+…,t∈[t0,T],

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.

### Notation

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

 (7a) 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 d=2. In (7a), ψ:Rd×[t0,T]→C:(x,t)→ψ(x,t) denotes the complex-valued macroscopic wave function and Δ=∂2x1+⋯+∂2xd the Laplace operator with respect to cartesian coordinates x=(x1,…,xd)∈Rd. In general, the angular momentum rotation term has the form ω′(t)⋅L(x)with direction of rotation given by the time derivative of a vector-valued function ω:[t0,T]→R3 and angular momentum operator defined by the cross product of position and gradient L(x)=−i(x×∇)∈R3. Commonly, it is assumed that the axis of rotation coincides with the x3-axis; thus, it suffices to consider a scalar time-dependent function ω:[t0,T]→R and set (7b) Further decisive quantitites in (7a) are the external potential V:Rd→R describing the confining trap and the nonlinear function f:R→R characterising the particle interaction; with regard to numerical simulations, we focus on the consideration of a quadratic potential involving positive weights V(x)=12d∑ℓ=1γ2ℓx2ℓ,γℓ>0,ℓ∈{1,…,d},x=(x1,…,xd)∈Rd, (7c) and a cubic nonlinearity f(|ψ(x,t)|)ψ(x,t)=ϑ|ψ(x,t)|2ψ(x,t),ϑ∈R,(x,t)∈Rd×[t0,T], (7d) describing the strength of short-range two-body interactions in the condensate.

### Generalisation

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

 ω(t)=Ωt,t∈[t0,T],Ω∈R; (7e)

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

 d=2:R(t)=(cos(ω(t))sin(ω(t))−sin(ω(t))cos(ω(t))),d=3:R(t)=⎛⎜ ⎜⎝cos(ω(t))sin(ω(t))0−sin(ω(t))cos(ω(t))0001⎞⎟ ⎟⎠,x=R(t)ξ,t∈[t0,T], (8a) and attain the non-autonomous nonlinear Schrödinger equation {i∂tφ(ξ,t)=−12Δφ(ξ,t)+V(R(t)ξ)φ(ξ,t)+f(∣∣φ(ξ,t)∣∣)φ(ξ,t),φ(ξ,t0)=ψ0(ξ),(ξ,t)∈Rd×(t0,T); (8b)

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

 {iu′(t)=H(t)u(t)=(A+B(t))u(t),t∈(t0,t0+h),u(t0)=u0, (9a) involving a constant square complex matrix A∈CM×M and a matrix-valued time-dependent function B:[t0,t0+h]→CM×M. In view of (4), it is justified to assume that the values of B commute [B(s),B(t)]=B(s)B(t)−B(t)B(s)=0,s,t∈[t0,t0+h], (9b)

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.

### Fundamentals

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

 c1=12−√1510,c2=12,c3=12+√1510. (10a) The Lagrange polynomials through these nodes are defined by L1:R⟶R:σ⟼σ−c2c1−c2σ−c3c1−c3,L2:R⟶R:σ⟼σ−c1c2−c1σ−c3c2−c3,L3:R⟶R:σ⟼σ−c1c3−c1σ−c2c3−c2; by construction, the following identities are fulfilled Lk(ck)=1,Lk(cℓ)=0,k,ℓ∈{1,2,3},ℓ≠k. We consider the associated interpolating function ˆH:R⟶CM×M:t⟼L1(t−t0h)H1+L2(t−t0h)H2+L3(t−t0h)H3,Hk=H(t0+ckh)=A+Bk,Bk=B(t0+ckh),k∈{1,2,3}, which indeed satisfies the relation ˆH(t0+ckh)=L1(ck)H1+L2(ck)H2+L3(ck)H3=Hk,k∈{1,2,3}. Substituting τ=t−t0−h2, a straightforward calculation yields the representation ˆH(τ+t0+h2)=A+B2−√153τh(B1−B3)+103τ2h2(B1−2B2+B3); consequently, Taylor series expansions B3−B1=O(h),B3−2B2+B1=O(h2), and the restriction to a suitably chosen interval verify the relations ˆH(τ+t0+h2)=i1h(α1+α2τh+α3τ2h2),τ∈[−h2,h2],α1=−ih(A+B2)=O(h),α2=−ih√153(B3−B1)=O(h2),α3=−ih103(B3−2B2+B1)=O(h3). (10b) For the following considerations, it is useful to employ a short notation for matrix commutators of the decisive quantities [kℓ]=[αk,αℓ]={0,k=ℓ,O(hk+ℓ),k≠ℓ,k,ℓ∈{1,2,3}; with regard to the design of modified CFQM exponential integrators, we note that condition (9b) implies [23]=0, whereas α1 in general does not commute with α2 or α3, respectively. Accordingly, nested matrix commutators are denoted by =[αi,[αj,[…,[αk,αℓ]…]]]=O(hi+j+⋯+k+ℓ),i,j,…,k,ℓ∈{1,2,3}. 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 {iv′(t)=ˆH(t)v(t),t∈(t0,t0+h),v(t0)=u0, coincide up to the desired order v(t0+h)−u(t0+h)=O(h7); hence, it suffices to construct suitable approximations to v(t0+h). 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 v(t0+h)=eΩ(h)u0,Ω(h)=−i∫t0+ht0ˆH(t)% dt−12∫t0+ht0∫tt0[ˆH(t),ˆH(s)]dsdt+…; due to the simple structure of the quadratic function ˆH, the arising integrals can be computed analytically and be expressed in terms of α1,α2,α3 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 t=τ+t0+h2 and inserting the representation (10b) yields −i∫t0+ht0ˆH(t)dt=−i∫h2−h2ˆH(τ+t0+h2)dτ=α1+112α3, and a similar calculation verifies −12∫t0+ht0∫tt0[ˆH(t),ˆH(s)]dsdt=−112[12]+1240[23]. Including terms up to order six Ω[6]=α1+112α3−112[12]+1240[23]+1360[113]−1240[212]+1720[1112], (10c) this procedure provides an appropriate approximation to v(t0+h) and hence to u(t0+h), that is Ω[6]−Ω(h)=O(h7),eΩ[6]u0−u(t0+h)=O(h7); (10d)

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

 u1=exp(−ihK∑k=1aJkH(t0+ckh))×⋯×exp(−ihK∑k=1a1kH(t0+ckh))u0,u1−u(t0+h)=O(hp+1); (11)

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

 u1−eΩ[p]u0=O(hp+1); (12)

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

 {iu′(t)=H(t)u(t),t∈(t0,t0+h),u(t0)=u0,

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

 ⎧⎪ ⎪⎨⎪ ⎪⎩iv′1(t)=K∑k=1a1kH(t0+ckh)v1(t),t∈(t0,t0+h),v1(t0)=u0,⋮⎧⎪ ⎪⎨⎪ ⎪⎩iv′J(t)=K∑k=1aJkH(t0+ckh)vJ(t),t∈(t0,t0+h),vJ(t0)=vJ−1(t0+h), (13a) and yields a numerical approximation through u1=vJ(t0+h),u1−u(t0+h)=O(hp+1). (13b)

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

 [W(⋅,t),W(⋅,s)]=0,[[Δ,W(⋅,t)−W(⋅,s)],W(⋅,t)−W(⋅,s)]=2d∑ℓ=1(∂ξℓ(W(⋅,t)−W(⋅,s)))2,s,t∈[t0,T], (14)

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

 (15)

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

 exp(x2α2+x3α3+y[212])

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

 u1=exp(−x12α2+x13α3+y[212])×exp(x21α1−x22α2+x23α3)×exp(x21α1+x22α2+x23α3)×exp(x12α2+x13α3+y[212])u0,x12=−x13=−160,x21=12,x22=−215,x23=140,y=143200.

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

 u1=e−ih(¯¯¯¯B4+h2˜B)% e−ih2(A+¯¯¯¯B3)e−ih2(A+¯¯¯¯B2)e−ih(¯¯¯¯B1+h2˜B)u0;

the decisive quantities are defined by

 c1=12−√1510,c2=12,c3=12+√1510,a11=10+√15180,a12=−19,a13=10−√15180,a21=15+8√1590,a22=23,a23=15−8√1590,¯¯¯¯B1=a11B(t0+c1h)+a12B(t0+c2h)+a13B(t0+c3h),¯¯¯¯B2=a21B(t0+c1h)+a22B(t0+c2h)+a23B(t0+c3h),¯¯¯¯B3=a23B(t0+c1h)+a22B(t0+c2h)+a21B(t0+c3h),¯¯¯¯B4=a13B(t0+c1h)+a12B(t0+c2h)+a11B(t0+c3h),˜B=125920[[A,B(t0+c3h)−B(t0+c1h)],B(t0+c3h)−B(t0+c1h)]. (16a) Accordingly to (13), it is convenient to reformulate the extension to non-autonomous Schrödinger equations (4) as {iv′1(t)=(¯¯¯¯B1+h2˜B)v1(t),t∈(t0,t0+h),v1(t0)=u0,{iv′2(t)=12(A+¯¯¯¯B2)v2(t),t∈(t0,t0+h),v2(t0)=v1(t0+h),{iv′3(t)=12(A+¯¯¯¯B3)v3(t),t∈(t0,t0+h),v3(t0)=v2(t0+h),{iv′4(t)=(¯¯¯¯B4+h2˜B)v4(t),t∈(t0,t0+h),v4(t0)=v3(t0+h); (16b) here, the matrix commutator is replaced by the Lie-bracket specified in (14). Numerical experiments indicate that this procedure yields a sixth-order approximation u1=v4(t0+h),u1−u(t0+h)=O(h7); (16c)

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

 {iu′(t)=H(t,u(t))=(A+B(t,u(t)))u(t),t∈(t0,t0+h),u(t0)=u0; (17a) essentially, the values of the solution and the defining operators correspond to u(t)=φ(⋅,t),u0=ψ0,A=−12Δ,W(⋅,t)=V(R(t)(⋅)),B(t,u(t))=W(⋅,t)+f(∣∣φ(⋅,t)∣∣). (17b)

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

### Approach

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

 [B(s,u(t)),B(σ,u(τ))]=0,s,t,σ,τ∈[t0,t0+h]. (18)

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

 ⎧⎪ ⎪⎨⎪ ⎪⎩iv′j(t)=K∑k=1ajkH(t0+ckh,vj(t)),t∈(t0,t0+h),vj(t0)=vj−1(t0+h), (19a) where v1(t0)=u0 and j∈{1,…,J}, yielding an approximation to the exact solution through u1=vJ(t0+h),u1−u(t0+h)=O(hp+1), (19b) see also (11) and (13); denoting bj=K∑k=1ajk,j∈{1,…,J}, the nonlinear evolution equations arising in (19a) corresponds to (19c)

with for