Since the kinetic and the potential energy term of the real time nonlinear Schr\"odinger equation can each be solved exactly, the entire equation can be solved to any order via splitting algorithms. We verified the fourth-order convergence of some well known algorithms by solving the Gross-Pitaevskii equation numerically. All such splitting algorithms suffer from a latent numerical instability even when the total energy is very well conserved. A detail error analysis reveals that the noise, or elementary excitations of the nonlinear Schr\"odinger, obeys the Bogoliubov spectrum and the instability is due to the exponential growth of high wave number noises caused by the splitting process. For a continuum wave function, this instability is unavoidable no matter how small the time step. For a discrete wave function, the instability can be avoided only for $\dt k_{max}^2{<\atop\sim}2 \pi$, where $k_{max}=\pi/\Delta x$.
Taha and Ablowitz 1 have shown for some time that the first order pseudo-spectral, split-operator method is a very fast way of solving the nonlinear Schrödinger equation. Bandrauk and Shen 2 later applied higher-order splitting algorithms with negative coefficients to solve the same equation. They regarded the nonlinear potential as timedependent. Since they can only estimate the intermediate-time nonlinear potential to second order, it is not proven that their higher-order algorithms actually converge at fourth or sixth-order. Recently Javanainen and Ruostekoski 3 have shown by symbolic calculations that fourth-order algorithms are possible by use of the "latest" intermediate wave function in evaluating the nonlinear potential. Strauch 4 , by constructing a special operator that correctly propagates the nonlinear potential term, proved that this use of the "latest" intermediate wave function is valid.
This work shows that: 1) Javanainen and Ruostekoski’s finding is a direct consequence of Taha and Ablowitz’ original work and a much simpler proof than that of Strauch is possible. 2) The time-dependent potential method of Bandrauk-Shen and the time-independent approach suggested by Javanainen and Ruostekoski both yielded identical second-order algorithms but different higher-order algorithms. 3) Verified numerically that algorithms derived by the time-independent method do converge to fourth-order when solving the Gross-Pitaevskii equation. 4) All such splitting algorithms possess a latent numerical instability which causes the wave function to blow up despite excellent total energy conservation. 5) The instability is shown to be due to the exponential growth of high wave number noises intrinsic to the splitting process. For a continuum wave function, this instability is unavoidable no matter how small is the time step. For a discrete wave function, this can only be avoided if ∆t < ∼ 2π/k 2 max , which forces ∆t to be very small if the discretization is very fine with a large k mas = π/∆x. The next three sections summarize how higher order algorithms can be systematically derived and Section V discusses the instability in detail.
Consider the nonlinear Schrödinger equation defined by
(2.1)
The free particle propagation can be solved exactly in operator form
where the operator T = -1 2 ∇ 2 . Since T is diagonal in k-space, (2.2) is usually solved by Fast Fourier Transforms (FFT). Surprisingly, as shown by Taha and Ablowitz, the potential part of the equation
can also be solved exactly
This is because (2.3) exactly conserves |ψ| 2 (multiply (2.3) by ψ * , the complex conjugated equation by ψ and subtract) and the nonlinear potential is just a constant in (2.3). This is also clear from (2.4),
since ψ(0) is only multiplied by a phase. Eq.(2.2) and (2.4) are the basic building blocks for constructing splitting algorithms for solving the nonlinear Schrödinger equation. Eq.(2.4) is the fundamental justification for using the “latest” wave function in computing the nonlinear potential 3 . (See also below). Define a time-independent operator V such that
Note that V only acts on |ψ(t) and does not act on its own eigenvalue g|ψ(t)| 2 . It follows that
The crucial point here is that V has no time-dependence, when it acts on any |ψ(t) , it produce the eigenvalue g|ψ(t)| 2 . The resulting time-dependence of the nonlinear potential is due entirely to the state vector |ψ(t) and not to the operator V . The exact solution can then be written in operator form as
For our purpose here, we only need to know (2.7) and not the explicit form of V . For an elegant, but rather abstract construction of V , see Strauch’s 4 recent work.
To solve (2.8) by splitting algorithms, one factorizes the evolution operator to any order with a suitable set of
where we have denoted ε = -i∆t. For example, we can have the second order algorithm 2A as
where according to (2.4) or (2.7), we must take
Algorithm 2A only requires one-pair of FFT (forward and backward) to achieve second-order accuracy, which is the same number of FFT needed for a first-order algorithm. If the nonlinear potential is treated as a time-dependent potential, as done by Bandrauk and Shen 2 , then we would have the algorithm 5,6 ψ(∆t) = e -i 1 2 ∆tV (∆t) e -i∆t T e -i 1 2 ∆tV (0) ψ(0).
In this case, since the last factor is only a phase,
the result is the same as (3.2). If one ignores the time-dependence 7 and uses V (∆t) = V (0) = g|ψ(0)| 2 , then algorithm (3.4) is degraded to first order. Similarly one has the second-order algorithm 2B,
where here
In the time-dependent potential approach, one would have instead,
One must now evaluate V (∆t/2) = g|ψ(∆t/2)| 2 . Since the algorithm is only second order, one can simply approximate the midpoint wave function to first order,
and therefore
Again, the result is the same as (3.6) For fourth and higher order algorithms, the time-dependent potential approach cannot be easily implemented. It is much more efficient to u
This content is AI-processed based on open access ArXiv data.
