A structure-preserving Fourier pseudo-spectral linearly implicit scheme for the space-fractional nonlinear Schr\"odinger equation

We propose a Fourier pseudo-spectral scheme for the space-fractional nonlinear Schr\"odinger equation. The proposed scheme has the following features: it is linearly implicit, it preserves two invariants of the equation, its unique solvability is guaranteed without any restrictions on space and time stepsizes. The scheme requires solving a complex symmetric linear system per time step. To solve the system efficiently, we also present a certain variable transformation and preconditioner.


Introduction
The Schrödinger equation, a fundamental equation in quantum mechanics, can be derived by using Feynman path integrals over Brownian trajectories [7]. Around 2000, Laskin coined the fractional Schrödinger equation by replacing the Feynman path integrals by the Lévy ones [13,14,15]. For example, the space-fractional nonlinear Schrödinger (FNLS) equation with the cubic nonlinearity is formulated as The fractional Schrödinger equation has found several applications in physics [11,12,18,32]. Besides, from mathematical viewpoints, several properties such as the well-posedness have been established [2,9,10]. Further, there has been a growing interest in constructing reliable and efficient numerical schemes. For the conventional NLS equation, it has been proved and is now wellknown that structure-preserving schemes such as symplectic integrators, multi-symplectic integrators and invariants-preserving integrators, have substantial advantages over general-purpose methods such as explicit Runge-Kutta-type methods with standard finite difference spatial discretization: structure-preserving schemes often produce qualitatively better numerical solutions over a long-time interval with a relatively large time stepsize (see, e.g. [4,6]). Therefore, it is no wonder that recent papers on the FNLS equation have mainly focused on the construction of structure-preserving schemes. For example, based on a Hamiltonian or generalized multisymplectic structure, symplectic or multi-symplectic schemes were presented in [30]. The fractional Schrödinger equation has mass (probability) and energy conservation laws, and several invariant(s)-preserving schemes have been developed: schemes preserving the mass [16,27,29] and schemes preserving both the mass and energy [5,17,28,31].
In this paper, we are concerned with invariants-preserving numerical schemes. The aforementioned numerical schemes preserving both invariants are nonlinear, and thus computationally expensive. Therefore, it is hoped that those schemes are linearized without deteriorating the mass and energy conservations. However, there are two challenges for this aim. First, standard approaches to the linearization of nonlinear schemes often lose conservation properties. In fact, linearly implicit schemes derived in [16,27] only preserve the mass. Second, even if we are succeeded in designing intended linearly implicit schemes, they result in solving dense linear systems due to the fractional Laplacian, which is nonlocal. Although some of the papers mentioned above contribute to the treatment of the fractional Laplacian in the discrete settings, no matter how we discretize the fractional Laplacian a matrix representing a discrete fractional Laplacian is a dense matrix. Thus, the computation of such dense linear systems is an important issue. In fact, the computational complexity of direct solvers applied to dense linear systems is of O(N 3 ), where N is the size of the system. Even if we consider iterative solvers, each iteration usually requires O(N 2 ) operations.
With these backgrounds, we aim to derive a linearly implicit numerical scheme preserving discrete mass and energy simultaneously, and then discuss the implementation issues for the proposed scheme. In our problem setting, imposing the periodic boundary condition, we mainly focus on the FNLS equation of the form [9] whereû k denotes the Fourier coefficients and µ = 2π/L. The mass (probability) and energy for the FNLS equation (1) are defined by respectively, and these quantities are constant along the solution [9].
We employ the pseudo-spectral method for the space discretization. In this approach, the treatment of the fractional Laplacian is rather straightforward, as will be explained in Section 2.2. Our idea of the time discretization is similar to the one proposed for the coupled FNLS equation [26]. However, our derivation is a bit more systematic, and to illustrate this we also show that invariants-preserving linearly implicit schemes can be constructed for the FNLS equation with stronger nonlinearity. We note that our approach is motivated by the discrete variational derivative method [3,8].
When we implement the proposed scheme, we need to choose a linear solver carefully. Krylov subspace methods seem suitable because the multiplication of a vector by the coefficient matrix can be efficiently computed with O(N log N ) operations thanks to the fast Fourier transform. The coefficient matrix of the linear system appearing in the scheme is found to be complex and symmetric (thus, non-Hermitian), and varies as time passes. For non-Hermitian systems, the Bi-CGSTAB method [24] is a standard choice, but this method requires two matrix-vector products per iteration. As alternative choices, we test the conjugate orthogonal conjugate gradient (COCG) method and conjugate orthogonal conjugate residual (COCR) method, which were specially designed for solving complex symmetric linear systems [23,25] and require only a single matrix-vector product per iteration. We will observe that these methods work preferably especially for a relatively small time stepsize. But recall that one of the advantages of structurepreserving schemes is that they can give qualitatively collect numerical solutions with relatively large time stepsizes. Therefore, the convergence of these methods for the case the time stepsize is not small enough is of interest. Unfortunately, however, our preliminary experiments show that if the time stepsize is relatively large, the coefficient matrix tends to be ill-conditioned, and the number of matrix-vector products required to the convergence tends to become large even if the COCG or COCR method is employed (although the results remain better than those by the Bi-CGSTAB method). To address such situations, we consider preconditioning issues. Note that the coefficient matrix involves the discrete Fourier transform, and thus we should not have each element of the coefficient matrix explicitly, which indicates that we do not have much choice of efficient preconditioners. The coefficient matrix can be seen as a dominant part plus perturbation term, and the diagonal of the dominant part seems the simplest preconditioner. However, this does not improve the convergence significantly since all the diagonal elements of the coefficient matrix are almost the same. Then, before discussing preconditioners, we apply a similarity transformation to the linear system by using a variable transformation which makes use of the discrete Fourier transform. Here, another difficulty arises: the proposed variable transformation makes the coefficient matrix of the transformed system non-symmetric. Therefore, the transformed linear system is out of the range of application of the COCG/COCR method, and we employ the Bi-CGSTAB method. The diagonal of the coefficient matrix of the transformed system is used as a preconditioner. We will observe the outstanding performance of the preconditioned Bi-CGSTAB method. But it is still of interest to study what happens if we aggressively apply the COCG/COCR method to the transformed non-symmetric linear system with the diagonal preconditioner because the non-symmetric part of the coefficient matrix of the transformed system can be regarded as a perturbation. We shall investigate this as well.
The paper is organized as follows. In Section 2, our notation and several preliminary results are summarized. In Section 3, our intended scheme is presented, and its properties are discussed. The preconditioning issues are also addressed. Several numerical studies are conducted in Section 4, and finally concluding remarks are given in Section 5.

Preliminaries
In this section, our notation and several preliminary results are summarized.
2.1. Mass and energy preservations for the FNLS equation. The mass and energy conservation laws follow in a straightforward way, but we here sketch the proof since it will be mimicked in the discrete settings in Section 3.
Theorem 1 (e.g. [9]). For the solution to the FNLS equation (1), it follows that Proof. First, we prove the mass preservation.
The first equality is just the chain rule. The fourth equality is due to the integration-by-parts formula: for any L-periodic complex-valued functions u and v, we have Next, we prove the energy preservation.
2.2.1. Notation. The period L is divided by N equal grids, which means ∆x = L/N . We denote the numerical solution for U (n∆t, k∆x) by U N −1 ) . To treat the periodic boundary condition, we consider {U k } k∈Z , an infinitely long vector, and then its N -dimensional restriction by the discrete periodic boundary condition: U k mod N for all k ∈ Z. The space to which such periodic vectors belong is denoted by

Discrete fractional Laplacian.
We define a discrete fractional Laplacian and show its properties. For simplicity, we assume that N is an odd number keeping in mind that the following discussion can be extended to an even number N straightforwardly. We employ the Fourier pseudo-spectral approach for the space discretization.
We define a function space S N by where g j (x) is a trigonometric polynomial defined by Note that g j (x l ) = δ j,l , where δ j,l denotes the Kronecker delta. We define the interpolation operator I N : L 2 (T) → S N by and thus where We then define a discrete fractional Laplacian (− ) By using the notation of the discrete Fourier transform and its inverse: the discrete fractional Laplacian can also be expressed as We shall abuse notation and write (− ) The next lemma shows that the discrete fractional Laplacian (− ) α 2 d is a Hermitian operator. This property will be used to prove the unique solvability of our proposed scheme.
where the symbol † denotes the Hermitian adjoint.
Proof. The operator D α defined by (6) is Hermitian (D † α = D α ) since it is a real diagonal matrix. Due to the unitarity of the discrete Fourier transform F d : Note that (− ) α 4 d is also Hermitian, and (− ) Thus, as an immediate consequence of Lemma 1, we have the following summation-by-parts formula that corresponds to the integration-by-parts formula (4).

Lemma 2. For any two vectors
Remark 1. It can be readily checked that F −1 d D α F d is a real matrix, but this property does not hold if we replace D α with a general real diagonal matrix.

Linearly implicit scheme
In this section, we present a linearly implicit scheme for the FNLS equation (1), show several properties of the scheme, and discuss the implementation issues with an emphasis on preconditioning issues.
3.1. Linearly implicit scheme. We propose the following linearly implicit scheme: given an initial approximation U (0) ∈ X d and starting approximation The initial approximation is usually set to U (0) k = u 0 (k∆x). The starting value U (1) can be prepared in several ways, and this issue will be discussed later. The scheme (8) can be written in the following form: where ) needs to be solved per each time step. We note that the coefficient matrix is complex and symmetric because F −1 d D α F d is a real symmetric matrix according to Lemma 1 and Remark 1.
Below, we show that the scheme (8) preserves the following discrete mass and energy: Note that though the discrete energy depends on two solution vectors, it seems a reasonable one since if V = U it becomes a more straightforward definitioñ Theorem 2. The scheme (8) conserves the discrete mass M d and the discrete energy H d in the sense that Proof. First, we prove the discrete mass preservation.
The first equality is just the factorization, which corresponds to the chain rule. In the second equality, we have substituted the scheme (8). The fourth equality is due to the summation-by-parts formula (7).
Next, we prove the discrete energy preservation.
In the second equality, we have used the identity aA − bB = (a + b)(A − B)/2 + (a − b)(A + B)/2, and in the third equality, we have substituted the scheme (8).
Remark 2. In the above discussion, defining the discrete quantities is the most crucial part, since if we employ other definitions, the corresponding schemes vary, which might be nonlinear.
In particular, the key to defining (10) is that it is at most quadratic with respect to both vectors U and V . For more details, we refer the reader to [3,8,19].
In general, while the computational complexity of the time integration of linearly implicit schemes is much cheaper than that of nonlinear schemes, the former often suffers from the strict restriction on stepsizes to ensure the solvability (see, e.g. [20]). Fortunately, the following theorem shows that the scheme (8) is unconditionally solvable without any restriction on the time and space stepsizes.
Theorem 3 (Unique solvability). For n ≥ 1, let U (n−1) ∈ X d and U (n) ∈ X d be given. Then, the scheme (8) has a unique solution U (n+1) without any restriction on the time and space stepsizes.
Proof. We show that the operator (9) is nonsingular, independently of U (n) . Since F −1 d D α F d is Hermitian due to Lemma 1, its eigenvalues are all real. Note that D(U (n) ) is a diagonal matrix and its elements are all real. Therefore, the real part of all eigenvalues of ) is 1, which indicates that the operator is nonsingular.
3.2. Starting procedure. The starting approximation U (1) can be computed in many different ways. Among several possibilities, we here present a fully nonlinear scheme: (k = 0, . . . , N − 1). Theorem 4. The scheme (12) conserves the discrete mass M d and the discrete energyH d in the sense that The proof is similar to that of Theorem 2 and thus is omitted.
As long as this nonlinear scheme (12) or a mass-preserving one-step scheme is used to obtain the starting approximation U (1) , it follows that M d (U (n) ) = M d (U (0) ) (n = 0, 1, 2, . . . ) for the solution to the linearly implicit scheme (8). However, even if the nonlinear scheme (12) is employed to obtain U (1) ,H d (U (n) ) =H d (U (0) ) in general. In Section 4.2, we shall numerically investigate to what extentH d (U (n) ) remains close toH d (U (0) ).

Remark 3.
Nonlinear schemes preserving either the mass or energy can also be derived, as presented in [1] for the case α = 2. There, the idea is to apply the average vector field method [22] to a semi-discrete scheme which is obtained based on a variational structure.

3.3.
Linearly implicit schemes for the FNLS with strong nonlinearity. In addition to the FNLS equation with cubic nonlinearity, which is the main focus of this paper, we briefly explain that similar structure-preserving linearly implicit schemes can be constructed for the FNLS equation with stronger nonlinearity. We consider where ρ ∈ {1, 2, . . . }. For this equation, while the mass (2) remains an invariant, the definition of the energy is modified to Note that this reduces to (3) when ρ = 1. We have already discussed the case ρ = 1 in the above subsections. Observe that the scheme (8) is a two-step method. If ρ ≥ 2, an intended scheme requires more steps so that the resulting scheme is linear in terms of the unknown solution vector. Our intended scheme is defined by which is a (ρ+1)-step method. This scheme is mass-preserving M d (U (n+1) ) = M d (U (n−ρ) ) and further energy-preserving in the sense that Note that this quantity is quadratic with respect to each vector U (i) (i = n − ρ, . . . , n). The proof of the preservations is similar to that of Theorem 2 and thus is omitted.

3.4.
Preconditioning. The coefficient matrix of the linear system (9) is complex and symmetric because F −1 d D α F d is a real symmetric matrix according to Lemma 1 and Remark 1. It is hoped that the linear system is solved efficiently per each time step.
As explained in Section 1, the multiplication of a vector by the coefficient matrix can be efficiently computed with O(N log N ) operations thanks to the fast Fourier transform, and hence Krylov subspace methods seem suitable choices for solving (9). For non-Hermitian systems, the Bi-CGSTAB method [25] is a standard choice, but this requires two matrix-vector products per iteration. As alternative choices, we test the conjugate orthogonal conjugate gradient (COCG) method [24] and the conjugate orthogonal conjugate residual (COCR) method [23], which were specially designed for solving complex and symmetric linear systems. These methods require only a single matrix-vector product per iteration. As will be seen later, they work better than the Bi-CGSTAB method. However, when relatively large stepsize or large N is used, both the COCG/COCR and Bi-CGSTAB methods tend to require a large number of iterations because the coefficient matrix becomes ill-conditioned, and this behaviour is problematic especially when we consider multi-dimensional problems. Therefore, it is worth considering the preconditioning issues for one-dimensional problems.
We need to understand why the convergence behaviour is deteriorated (some numerical results will be illustrated in Section 4.3). To simplify the notation, we rewrite (9) as where D(U ) = diag(|U 0 | 2 , |U 1 | 2 , . . . , |U N −1 | 2 ) and D α is given in (6). The smallest and largest eigenvalues of the sum of the first two terms I + i∆tF −1 d D α F d are 1 and 1 + i∆t|µ(N − 1)/2| α , respectively, which indicates that the sum of the first two terms tends to be ill-conditioned for large α, N and ∆t. On the other hand, since the third term, i.e. the diagonal matrix i∆tD(U ), represents the shape of the solution, its condition number remains nearly unaffected by the choice of N . Therefore, the third term can be seen as a perturbation. As will be illustrated in Section 4.3, the convergence is deteriorated for large α, N and ∆t, and the sum of the first two terms is certainly the cause of the deterioration. Note that for large ∆t the influence of the last term could be substantial, and this influence further worsens the convergence.
Below, we consider the preconditioning issues. Our idea is a combination of a certain variable transformation and preconditioner.
For Ax = b, the preconditioned COCG and preconditioned Bi-CGSTAB methods with a matrix (preconditioner) M are summarized in Algorithms 1 and 2, where (x, y) = x y. The preconditioner M is ideally chosen such that M −1 A has a smaller condition number than A. The simplest approach is to extract the diagonal part of the coefficient matrix as a preconditioner M . However, in our situation, each diagonal element in the coefficient matrix of (13) has almost the same value, because all diagonal elements of the dominant part I + i∆tF −1 d D α F d are the same. Therefore, this approach does not improve the convergence behaviour (in fact our preliminary numerical experiments support this discussion). Let us consider a variable transformation y = F d x. By a similarity transformation, the linear system (13) can be rewritten as For the solution to the transformed linear system (14), the following property holds. Proposition 1. The relative 2-norm residual for (14) coincides with that for (13).
Proof. Letx andỹ be approximations to (13) and (14), respectively, with the relationỹ = F dx . The exact solutions are denoted by x * and y * . Then, it follows that due to the unitarity of F d . p n+1 = r n + β n (p n − ω n v n ) This property indicates that we only have to monitor the relative error for the transformed linear system (14) during the iterations and calculatex = F −1 dỹ only after the error meets the convergence criteria. As a preconditioner, we use It should be noted that the coefficient matrix of the transformed linear system (14) is still complex but no longer symmetric because F d D(U )F −1 d is not real (it is skew-Hermite). Therefore, the transformed linear system is out of the range of application of the COCG/COCR method. Therefore, it is natural to apply the Bi-CGSTAB method to the transformed linear system. We will observe the performance of the preconditioned Bi-CGSTAB method. We note that although the coefficient matrix of the transformed linear system is no longer symmetric, the sum of the first two terms I + i∆tD α remains symmetric and the third term −i∆tF d D(U )F −1 d could be regarded as a perturbation. It is thus still of interest to investigate what happens if we aggressively apply the COCG/COCR method to the transformed linear system (14) with the preconditioner (15). In Section 4.3, we will also see how this approach works.

Numerical experiments
We now test our linearly implicit scheme (8). First, we check the qualitative behaviour of the numerical solutions obtained from the linearly implicit scheme (8) by comparing the results with those from the fully nonlinear scheme (12). We also check to what extentH d (U (n) ) remains close toH d (U (0) ). Next, we discuss the efficiency of the scheme with particular emphasis on linear solvers, and observe how the preconditioned methods work.
All the computations were performed in a computation environment: 1.6 GHz Intel Core i5, 16GB memory, OS X 10.14. We use Julia version 1.1.1.

Numerical behaviour.
We check the numerical behaviour of the proposed scheme (8) and compare the results with reference solutions which are computed by the nonlinear scheme (12). The linear system in (9) is solved by the conjugate orthogonal conjugate gradient (COCG) method [25] in this subsection (Algorithm 1 summarizes the preconditioned COCG method, but it reduces to the COCG method with M = I). As an initial guess of the COCG method, we use the numerical solution at the current time step. The convergence criteria is set to 10 −10 in terms of the relative 2-norm residual. The nonlinear system (12) is solved by nlsolve with the tolerance 10 −10 .
As an example, we set L = 20, N = 101 and ∆t = 0.02. The initial value is set to u 0 (x) = 2 exp(0.5ix) sech( √ 2(x − 10)), which is a snapshot of a solitary wave solution for the case α = 2. Figs. 1, 2 and 3 show the contour of the absolute value of numerical results for several α obtained by the linearly implicit scheme (8) and the nonlinear scheme (12). It is observed that the linearly implicit scheme exhibits qualitatively comparable results to the expensive nonlinear scheme. Figs. 4, 5 and 6 show errors of the discrete mass and energy. For the discrete mass |M d (U (n) ) − M d (U (0) )| is plotted and for the discrete energy |H d (U (n+1) , U (n) ) − H d (U (1) , U (0) )| is plotted. Both discrete quantities are well-preserved as expected (note that the tolerance of the linear and nonlinear solvers are set to 10 −10 and the scheme is computed 1.25 × 10 4 times until t = 250).    valueH d (U (n) ) could be a good barometer when we consider the long-time stability. It seems quite challenging to obtain an a priori estimate for the scheme (8), and thus we consider this numerically. Fig. 7 shows the results for the case α = 1.6. It is observed that |H d (U (n) ) − H d (U (n) )| is bounded by 10 −2 when t ≤ 300. However, when t exceeds 300, the error becomes large with strong oscillation, and thus we cannot expect an error bound forH d (U ) for all t. With other choices of parameters, qualitatively similar behaviour is observed. These observations indicate that the instability might be caused for a very long-time integration. This could be a drawback of the proposed linearly implicit scheme compared with the nonlinear scheme (12). However, we emphasize that sinceH d (U (n) ) is easily monitored during the time integration, and we could easily detect a sign of the instability.

4.3.
Performance of the preconditioning. We here discuss the performance of linear solvers for solving (9).
First, we consider solving the original system (9) by the COCG method. Fig. 8 shows the number of iterations required to the convergence for several settings. It is observed that more iterations are required for large α, N , ∆t. In particular, there is a significant gap when we change N or ∆t. However, we would also like to emphasize that even in the worst case (α = 2, N = 401 and ∆t = 0.02), the result is much better than that by the Bi-CGSTAB method, which is illustrated in Fig. 9. We thus conclude that when we solve the original system (9) directly, the COCG method seems an appropriate choice. We also note that the COCR method gives comparable results to the COCG method.   Although the COCG method is preferred for solving (9), it is hoped that the convergence behaviour is improved. Thus, we next discuss how the variable transformation and preconditioner proposed in Section 3.4 work. In the following numerical experiments, we consider the case α = 2. Table 1 shows the maximum, minimum and average number of iterations of the preconditioned Bi-CGSTAB method for several N , where the time stepsize is set to ∆t = 0.02 and the initial value u 0 (x) = 2 exp(0.5ix) sech( √ 2(x − 10)). The convergence behaviour is outstandingly improved for the case N = 401 compared with Fig. 9, and furthermore, the results are notable in that for all cases the preconditioned Bi-CGSTAB method requires only three iterations. Let us change the time stepsize to ∆t = 0.2. The results are shown in Table 2. By comparing Table 2 with Table 1, we observe that the convergence behaviour depends on ∆t, but the results still remain outstanding. Let us also change the initial value. The results are shown in Table 3, which indicate that the dependency on the shape of solutions is subtle. Table 1. The maximum, minimum and average number of iterations of the preconditioned Bi-CGSTAB method: the time stepsize is set to ∆t = 0.02, and the initial value u 0 (x) = 2 exp(0.5ix) sech( √ 2(x − 10)).
N 401 1001 4001 maximum 6 6 6 minimum 5 5 5 average 5.020 5.020 5.020 Table 3. The maximum, minimum and average number of iterations of the preconditioned Bi-CGSTAB method: the time stepsize is set to ∆t = 0.02, and the initial value u 0 (x) = 2 exp(0.5ix) sech(x − 10). As discussed in Section 3.4, it is also of interest to investigate the behaviour when the COCG method is aggressively applied to the transformed system (14) with the preconditioner (15), since the coefficient matrix in (14) can be seen as a complex symmetric matrix plus a perturbation. Fig. 10 shows the results. From the left figure, it is observed that, when ∆t = 0.01 and 0.02, the preconditioned COCG method actually work and the results are significantly improved compared with the bottom figure in Fig. 8.
Unfortunately, however, if we use a larger time stepsize ∆t = 0.05, the preconditioned COCG method requires 50 iterations at the 7th time step as shown in the right figure of Fig. 10, and the iteration does not converge within 1, 000 iterations at the 8th time step. This observation indicates that with the stepsize ∆t = 0.05 the influence of the perturbation term is not negligible.
Let us change the initial condition to u 0 (x) = 2 exp(0.5ix) sech(4(x − 10)). This function has a steeper slope. The results are displayed in Fig. 11. It is observed that even if a much larger time stepsize ∆t = 0.2 is employed, the preconditioned COCG method works fine. Conversely, a more gradual initial condition u 0 (x) = 2 exp(0.5ix) sech(x − 10) with the time stepsize ∆t = 0.02 was also considered as our preliminary experiments, and it was observed that in this case the preconditioned COCG method did not converge at 7th time step. These observations indicate that the convergence of the preconditioned COCG method strongly depends on ∆t and the shape of the solution (in other words, the influence of D(U (n) )). They make the effect of the perturbation term in (14) significant.
For the problem considered in this paper, it is highly recommended to use the Bi-CGSTAB method with the proposed variable transformation and preconditioner, but it is of interest to understand the behaviour of the preconditioned COCG method, which will be investigated in our future work.  Figure 11. The number of iterations for the preconditioned COCG method applied to (14) with the matrix M = I + i∆tD α . The parameters are set to α = 2 and N = 401.

Concluding remarks
In this paper, we proposed the linearly implicit scheme (8) for the FNLS equation preserving two invariants: mass and energy. The scheme exhibited qualitatively comparable results to the expensive nonlinear scheme (12). The preconditioning issues were also discussed: the preconditioned Bi-CGSTAB method is a preferable choice.
We note several directions for future work.
• It is hoped that the proposed scheme is used to investigate more challenging problems such as multi-dimensional problems. When we consider the multi-dimensional problems, the computational complexity and preconditioning issues become increasingly important, and thus the discussion on the preconditioning considered for the one-dimensional problem would be helpful. • The linearly implicit scheme does not preserveH d (U ), which is defined on a numerical solution of a single time step. The results presented in Fig. 7 should be theoretically investigated in more detail. Note that structure-preserving linearly implicit schemes preserving a certain quantity have also been proposed for other partial differential equations (see, e.g. [3,20,21]), and similar behaviour might also have to be reconsidered as well. • The discussion in Section 4.3 indicates that the COCG/COCR method is applicable to complex but non-symmetric matrices if the non-symmetric term can be regarded as a perturbation in some sense. This behaviour will be investigated theoretically in more detail.