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{\"o}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.


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 ∈ (t 0 , T ) ; (1) 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 F : [t 0 , T ] → C M ×M are found in [2,3,8,13,15,16], 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 2 P. BADER, S. BLANES, F. CASAS, M. THALHAMMER 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 t * ∈ [t 0 , T ]; 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 [18,27]  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 [8]; 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 B represents certain linear combinations of the values of the potential and B 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), A represents the Laplace operator, B a real-valued space-timedependent trapping potential, and C 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.
(i) 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 r compositions and a CFQM exponential integrator involving s compositions, this strategy requires the numerical solution of r 2 s linear subequations by Fourier spectral space discretisation as well as r + r 2 s linear subequations by pointwise multiplication, since application of r-stage splitting method: application of s-stage CFQM exponential integrator: application of r-stage splitting method: r 2 s times (Fourier spectral method): i u (t) = A u(t) , r 2 s times (pointwise multiplication): i u (t) = B(t * ) u(t) .
(ii) 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 s * , t * ∈ [t 0 , T ] by Fourier spectral space discretisation and pointwise multiplication, respectively, each rs times, since application of s-stage CFQM exponential integrator: application of r-stage splitting method: rs times (Fourier spectral method): i u (t) = A u(t) , rs times (pointwise multiplication): i u (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 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 [L(τ ), L(σ)] corresponds to a secondorder 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 [9,22] and the transformation to the rotating frame combined with the application of exponential integrators as proposed in [5,6,7,10,11]. 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 [22]; 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 [10] 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 Ω ∈ R denotes the constant angular velocity; by contrast, in Section 3, the function Ω : [t 0 , t 0 + h] → C M ×M denotes the exponent in the Magnus expansion.

Rotational Gross-Pitaevskii equation.
Rotating Bose-Einstein condensates. The phenomenon of Bose-Einstein condensation, for the first time observed in 1995 [4,19,21], 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 [1,25] 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 (d = 3), 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 d = 2. In (7a), denotes the complex-valued macroscopic wave function and ∆ = ∂ 2 x1 + · · · + ∂ 2 x d the Laplace operator with respect to cartesian coordinates x = (x 1 , . . . , x d ) ∈ R d . 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 ω : [t 0 , T ] → R 3 and angular momentum operator defined by the cross product of position and gradient L(x) = − i (x × ∇) ∈ R 3 . Commonly, it is assumed that the axis of rotation coincides with the x 3 -axis; thus, it suffices to consider a scalar time-dependent function ω : [t 0 , T ] → R and set Further decisive quantitites in (7a) are the external potential V : R d → 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 and a cubic nonlinearity 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 V : R d × [t 0 , T ] → R and a nonlinear term of the form f (|ψ(x, t)|). The formulation (7) in particular complies with [10,32] for the special case of constant angular velocity (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 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 V :

Quasi-Magnus exponential integrators.
In this section, we introduce CFQM exponential integrators [2,15] and modified CFQM exponential integrators [8] 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 h > 0.
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 A ∈ C M ×M and a matrix-valued timedependent function B : [t 0 , t 0 + h] → C M ×M . In view of (4), it is justified to assume that the values of B commute whereas the matrix commutator [A, B(t)] 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 nonautonomous 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 [13]. Additional considerations lead to modified CFQM exponential integrators for nonautonomous 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 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 ; 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 [2 3] = 0, whereas α 1 in general does not commute with α 2 or α 3 , respectively. Accordingly, nested matrix commutators are denoted by As shown in [12], 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 hence, it suffices to construct suitable approximations to v(t 0 + h). For a nonautonomous linear differential equation, a formal solution representation by the matrix exponential based on the Magnus expansion [26] 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 H, the arising integrals can be computed analytically and be expressed in terms of α 1 , α 2 , α 3 as well as nested commutators, see also [14,23,28]. As an illustration, we specify the leading contributions; substituting t = τ + t 0 + h 2 and inserting the representation (10b) yields and a similar calculation verifies Including terms up to order six this procedure provides an appropriate approximation to v(t 0 + h) and hence to u(t 0 + h), that is because of the symmetry of the Magnus expansion and the quadrature rule, only odd terms appear in Ω [6] . 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 [2 3] = 0 is only used for modified CFQM exponential integrators. CFQM exponential integrators. A pth-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 J ∈ N, associated coefficients a jk for (j, k) ∈ {1, . . . , J} × {1, 2, 3} are determined such that the corresponding approximate solutions u 1 agree with e Ω [6] u 0 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 [15]. We note that a higher number of exponentials J ∈ N 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 p = 2, 4, 6; for details, we refer to Section 5. Modified CFQM exponential integrators. As indicated in Section 1, for nonautonomous 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 d-dimensional Laplacian ∆ = ∂ 2 ξ1 + · · · + ∂ 2 ξ d and a space-time-dependent multiplication operator 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 [2 3] vanishes; as a consequence, relation (10c) and hence the consistency condition (12) reduce to Moreover, due to the fact that α 2 , α 3 , and the nested commutator [2 1 2] are of similar type, it is advantageous to incorporate terms of the form with appropriately chosen coefficients x 2 , x 3 , y ∈ R into the composition (11).
Sixth-order scheme. A particularly efficient sixth-order modified CFQM exponential integrator has been proposed in [8]; 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 10 , a 11 = 10+  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 (16c) a rigorous error analysis remains an open question.

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 h > 0 and to reformulate it as abstract evolution equation essentially, the values of the solution and the defining operators correspond to (17b) We recall that arguments detailed in Section 1 justify the consideration of a timeindependent linear part A and a time-dependent nonlinear part B; 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 A ∈ C M ×M and a matrixvalued time-dependent function B : [t 0 , t 0 + h] × C M → C M ×M . 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 v 1 (t 0 ) = u 0 and j ∈ {1, . . . , J}, yielding an approximation to the exact solution through see also (11) and (13); denoting the nonlinear evolution equations arising in (19a) corresponds to with ϕ j (·, t) = v j (t) for t ∈ [t 0 , t 0 + h] and j ∈ {1, . . . , J}.
Full discretisation. Favourable space and time discretisation methods for (19a) and (19c), respectively, are provided by high-order operator splitting methods combined with Fourier spectral space discretisations, realised by fast Fourier techniques.
They rely on the repeated numerical solution of subequations of the form ; we include the formulations as abstract evolution equations and the associated partial differential equations. The computational cost for the resolution of subequation (A) amounts to the application of a fast Fourier transform and its inverse; a characteristic invariance property permits to reduce the nonlinear subequation (B) to the linear subequation j (ξ, t 0 ) given , (ξ, t) ∈ R d × (t 0 , t 0 + h) , and to resolve it by pointwise multiplication. Altogether, a pth-order approximation to u(t 0 + h) results from a composition involving the flows of the subequations and suitably chosen coefficients In essence, our implementation of CFQM exponential integrators is based on this formulation, i.e., our algorithm calls subroutines Φ depending on the coefficients (a jk ) K k=1 and b j , which yield approximations to the solutions of the subequations with adjusted time increments α h and β h. Additional information on time-splitting Fourier spectral space discretisation methods for autonomous Gross-Pitaevskii equations is found for instance in [30,31], see also Section 5. Due to the presence of a trapping potential, the values of the macroscopic wave function tend to zero; hence, a suitable restriction of the space domain and the application of the Fourier or, alternatively, Sine spectral method is justified, see [31].
Sixth-order scheme. Based on the above considerations, it is straighforward to extend the sixth-order scheme (16) to the rotational Gross-Pitaevskii equation. We employ the abbreviations 10 , a 11 = 10+ , see also (14); for the special case of a quadratic potential, the arising spatial derivatives are given in Appendix B. As a consequence, we obtain (20b) the solutions to the first and fourth subequations are computed by pointwise multiplication, and the second and third subequation again correspond to autonomous Gross-Pitaevskii equations.
Future work. We believe that it is worth to investigate our approach in the context of ground state computations [20,32]; the use of exponential integrators involving complex coefficients with positive real parts guarantees stability of the approximations.

Numerical results.
In this section, we illustrate the performance of the proposed full discretisation methods for the time-dependent Gross-Pitaevskii equation with rotation term in two and three space dimensions. We provide numerical comparisons for CFQM exponential integrators and a modified CFQM exponential integrator; the numerical realisation relies on operator splitting and Fourier spectral space discretisation methods. We in particular implement the following standard and optimised schemes, see also   Linear and nonlinear test equations. We consider the rotational Gross-Pitaevskii equation (7)-(8) in two and three space dimensions, respectively, with initial and final times in order to retain the nonstiff orders of convergence, we impose a smooth Gaussianlike initial state For the time integration of the rotational Gross-Pitaevskii equation, we apply the above listed exponential integrators of overall orders two, four, and six; a numerical reference solution is determined by a sixth-order exponential integrator for the same spatial resolution and a refined time stepsize. In Figures 2-4, we illustrate the efficiency of the considered exponential time integration methods for constant angular velocity Ω = 0.5 and for increasing values of ϑ; in the special case ϑ = 0, the rotational Gross-Pitaevskii equation reduces to a non-autonomous linear Schrödinger equation. We in particular display the global errors versus the number of fast Fourier transforms and their inverse, which  we consider to be the most costly part of the computations; the slopes of the lines reflect the temporal orders of convergence. The obtained results for two and three space dimensions differ marginally. In all cases, it is observed that the modified CFQM exponential integrator (20) is beneficial, in particular, when high accuracy is desired; this is explained by the fact that optimised CFQM integrators of orders four or six, respectively, require the solution of three or six, respectively, autonomous Gross-Pitaevskii equations in each time step, while the modified CFQM exponential integrator involves only two autonomous Gross-Pitaevskii equations. Solution profile. In order to illustrate the dynamical behaviour of rotating Bose-Einstein condensates, we consider the rotational Gross-Pitaevskii equation (7) Figure 5. Time integration of two-dimensional rotational Gross-Pitaevskii equation (22) by sixth-order modified CFQM exponential integrator (20) and optimised sixth-order splitting method proposed in [18]. Solution profile |ψ(x, t)| 2 displayed for spatial section x ∈ [−5, 5] 2 and time t = 15. Consistent results obtained by time-splitting generalised-Laguerre-Fourier spectral methods [22, Eq. (7)-(8), Fig. 2].
this problem has been studied in [22] in the context of the generalised-Laguerre-Fourier spectral method. As expected, the solution profile displayed in Figure 5 is coincident with the result in [22, Fig. 2]. However, with respect to efficiency, our approach based on exponential integrators and Fourier spectral approximations outperforms the former approach based on the generalised-Laguerre-Fourier spectral method; moreover, its extension to the most relevant case of three space dimensions is straightforward.

Conclusions.
In this work, we proposed efficient exponential time integration methods for the Gross-Pitaevskii equation with rotation term in two and three space dimensions; we used that a transformation to the rotating frame is straightforward by considering the associated classical Hamiltonian and employed a natural approach that combines CFQM exponential integrators with operator splitting and Fourier spectral space discretisation methods. On the basis of a sixth-order modified CFQM exponential integrator, we demonstrated that the special structure of the resulting nonlinear Schrödinger equation allows the design of tailored schemes with reduced computational cost; furthermore, in contrast to time-splitting generalised-Laguerre-Fourier-Hermite spectral methods, our approach permits the implementation of the physically relevant case of three space dimensions as a direct extension of the two-dimensional case. Intrinsic advantages of CFQM exponential integrators and modified CFQM exponential integrators for non-autonomous linear differential equations are their favourable geometric properties, see [24]; in the advanced context of nonautonomous linear Schrödinger equations, this is reflected in an excellent stability behaviour in the underlying Lebesgue space of square-integrable functions. In combination with optimised operator splitting methods and fast Fourier techniques, the proposed exponential integrators lead to efficient high-accuracy discretisations in space and time; moreover, as described and illustrated in the present work, this methodology extends to non-autonomous nonlinear evolution equations.
Proceeding with a rigorous stability and error analysis, both for modified CFQM exponential integrators applied to non-autonomous linear Schrödinger equations and for CFQM exponential integrators applied to non-autonomous nonlinear Schrödinger equations, will be the subject of future work.