In this work, we are concerned with the following one-dimensional disordered nonlinear Schrödinger equation (D-NLS) on a torus or whole space which contains a spatial random potential [17, 18, 32, 51]:
where or (a torus), is the unknown, is the given initial data, is a given random potential, and is a given parameter with for defocusing nonlinear interaction and for focusing case. The D-NLS (1.1) is of paramount importance in mathematical and physical studies. For example, mathematically, it can be considered as the complex parabolic Anderson model [17, 18]. Physically, it has been considered for modeling the Anderson localization of Bose-Einstein Condensates [38, 43] and for the study of nonlinear dispersive wave dynamics in a disordered medium [15, 27, 32]. The spatial random potential , as considered in many related numerical works and physical applications in the literature, is composed out of a large number of random independent and identically distributed (i.i.d.) Fourier components [4, 13, 23, 32, 34, 38, 43, 52].
The D-NLS (1.1) in fact has received many mathematical attentions and we refer to [5, 12, 17, 18, 19, 26, 30] for intensive theoretical studies from different aspects. It has a discrete analogue known as the nonlinear Schrödinger equation on the disordered lattice [28, 53]:
where is a given parameter and
are a series of i.i.d. random variables. This discrete model has been widely addressed (more than (1.1) so far) in the physical studies, e.g. [2, 22, 25, 41, 46, 47, 50]. As a matter of fact, it shares many physical properties  with the continuous version (1.1). For the linear regime of (1.1) or (1.2), i.e. , it is well-understood as the famous Anderson localization  which tells that spreadings of any localized waves will be suppressed by the random potential. For the nonlinear regime in (1.1) or (1.2
), the competition between the random potential and the nonlinearity makes the long-time dynamics very complicated. For an initially localized wavepacket in the whole space setup, heuristically, the spreadings in the nonlinear system (1.1) or (1.2) could still be suppressed by the random potential, otherwise the nonlinearity would become negligible as time gets large and then the problem goes back to the classical Anderson localization . The theoretical understandings of this issue are quite challenging and so far the rigours results are still limited [11, 22, 51, 53]. Therefore, many existing works rely on the numerical simulations [22, 25, 32, 41, 46, 47, 50], and the majority of them are made on the above discrete version (1.2). While, the numerical research outcomes so far conflict with the heuristical argument and unlimited spreadings were observed [25, 41].
From the numerical point of view, there is some significant difference between solving the continuous D-NLS model (1.1) and the discrete one (1.2). In the existing numerical studies in the literature, e.g. [25, 41, 46], some classical symplectic integrators including the splitting schemes and finite difference schemes of the second or the fourth order accuracy were used to solve the discrete model. The high order accuracy and the symplecticity of the schemes on the ODEs could provide reliable long-term approximations of (1.2). However, it is far from straightforward to extend these numerical works from the discrete model (1.2) to the continuous model (1.1). The challenging issues for solving (1.1) include the domain truncation of whole space problem for long-time computing, the efficient and accurate spatial discretization, the long-time behaviour of the fully discretized integrators on the PDE, and even the accurate temporal approximation up to a fixed finite time, where the last issue is what we shall address in this paper. For attempts to control the temporal approximation error to arbitrary time, we refer to the work . For the studies of the nonlinear Schrödinger equations with time-dependent random/rough potentials, we refer to [4, 12, 13, 16, 34, 42].
For the popular classical numerical integrators such as Strang splitting method or Crank-Nicolson finite difference method, when applied to the continuous D-NLS model (1.1) as we shall see in our numerical experiments, their accuracy order could be much less than the optimal second order rate. This is because the regularity of the spatial noise could be very low (even in negative Sobolev space in some cases [10, 52]). This directly illustrates the problem of using high order splitting schemes for integrating (1.1), that the high order commutators [40, 49] generate high order differentiations on , and so the unboundedness of such terms causes significant accuracy order reduction. For instance, without the boundedness of , the Strang splitting method will fail to offer the second order accuracy. Moreover, the potential will induce roughness to the solution of (1.1). Recently, when the potential
is assumed to be a Gaussian white noise, the well-posedness of the Cauchy problem (1.1) in torus has been established only in space , and its well-posedness in the whole space case has later been given in . We will show in our study of (1.1) that for and smooth initial data , the regularity of the solution is just in . Such solution with low regularity also causes order reductions for classical integrators. For instance, the Crank-Nicolson method requests the boundedness of the six order spatial derivatives of the solution to reach its second order accuracy. For a rigorous analysis of Crank-Nicolson method, we refer to . To overcome the order reduction problem due to the lack of regularity in the solution, some state-of-the-art called the low-regularity Fourier integrators (LRI) have been proposed very recently for the nonlinear Schrödinger equations without the potential [14, 33, 36, 37], where the roughness was introduced through the initial data. The particular efforts of LRI were made to integrate the nonlinearity in physical space by losing two spatial derivatives for the second order accuracy in time. In the contrast to the setup in [14, 33, 36, 37], many physical problems related to (1.1) consider a smooth initial input [32, 38, 43, 47].
The purpose of this work is to develop the LRI method for solving the D-NLS (1.1), and we shall consider the random potentialor . We shall first review some classical methods and address their order reductions on (1.1) with more details. Then to raise the accuracy order, we consider the LRI method from  and extend it to integrate (1.1). Special efforts are made here to integrate the potential term in (1.1) up to the second order accuracy by losing two spatial derivatives, and to keep the final form of the scheme explicitly defined in the physical space so that it can be programmed efficiently under the Fourier pseudo-spectral method for spatial discretization. Meanwhile, efforts are also made to avoid CFL-type stability conditions in the scheme, which is very important since the roughness in space requires very small spatial mesh size for accurate approximations. We will prove the quadratic convergence in time of the proposed LRI in -norm for . For rougher potential case, e.g. , many numerical experiments are done to show that the expectation of the -error of LRI converges in time at a rate between one and two. Comparisons are made with the classical methods.
The rest of the paper is organised as follows. In Section 2, we review some classical integrators for nonlinear Schrödinger equations. In Section 3, we derive the LRI scheme for solving (1.1) with convergence analysis in the case . Numerical experiments with general random potentials are given in Section 4, and some conclusions are drawn in Section 5.
2. Classical numerical methods
In this section, we will briefly review some popular numerical integrators including the splitting methods and the finite difference methods for solving the nonlinear Schrödinger equations. In fact, these methods have been considered to integrate the discrete model (1.2) in the physical literatures, and here we present them for the continuous D-NLS (1.1).
For the computational reason, we would always consider the torus setup of (1.1), i.e. . By choosing the domain size large enough, this is also an approximation of the whole space case when the waves are away from the boundary during the dynamics. We denote as the time step, as the time grids and for with . For simplicity, to present the schemes and some known convergence results, let us consider a sampled or deterministic that could potentially be very rough.
2.1. Splitting scheme
For both of the subproblems, we have the exact solutions:
Then through composition, the classical Lie-Trotter splitting method for (1.1) reads:
which is of the first order accuracy in time under usual smooth enough setup . By further requesting more smoothness on the initial data and the potential in (1.1), with more times of composition, one could in principle get the high order accurate splitting methods [9, 29, 31, 35, 36, 40, 49]. Among the second order methods, the most popular one is the Strang splitting method:
It is well known that the local truncation error of Strang splitting (2.1) is determined by the double Lie commutator of the kinetic part and the potential part of the equation (1.1), i.e. with notations from ,
That is to say the truncation error of Strang splitting for (1.1) contains the terms and . In order to reach the optimal accuracy of the Strang splitting scheme (2.1), one will then ask for enough smoothness on both the potential and the initial data in (1.1). By using the analysis from  without essential modifications, for a given , the global approximation error of Strang splitting scheme (2.1) for solving (1.1) up to a fixed final time reads
For higher order splitting schemes and their convergence results under smooth potential and initial data case, we refer to .
The Strang splitting scheme (2.1) has some analogies such as the SBAB method , which are also second order accurate under smooth setup. These splitting methods have been applied to integrate the discrete model (1.2) with uniform distribution random potential on a very large time interval in [25, 41, 46], and unlimited subdiffusions of initially localized waves were observed. Numerically, the performance of the SBAB method is similar to the Strang splitting, so here we just consider the Strang splitting as a representative from this class of methods. The spatial discretizations of (2.1) can be made by the Fourier pseudo-spectral method  thanks to the periodic setup.
2.2. Finite difference method
As the most traditional numerical methods, different kinds of finite difference discretizations have been proposed for the nonlinear Schrödinger equations [1, 3, 21, 44, 54]. Here, we consider a second order semi-implicit finite difference integrator (FD) for (1.1):
with the starting value obtained from a first order version:
The implicit and explicit treatments in the above scheme are to avoid CFL-type stability conditions and meanwhile to make it efficient for programming under the Fourier pseudo-spectral method for spatial discretizations.
This FD method has also been applied to solve the discrete model (1.2) in [25, 41, 46], where the same long-term dynamics were observed as from the splitting schemes. By Taylor’s expansion, it is clear that the local truncation error in (2.2) depends on and . The boundedness of these two terms, under the smooth potential case, is equivalent to by the equation (1.1). So as has been established in [1, 7, 21, 44, 54], the global error of the FD method (2.2) for solving (1.1) up to a fixed is
For the rough potential case of , a rigorous first order convergence result has been established in .
We will apply the above presented second order classical methods to integrate the continuous D-NLS model (1.1), and study in particular their practical accuracy order in Section 4 by numerical experiments. As we shall see, all of them will suffer from accuracy order reductions, and so this gives us no motivation to include higher order methods here for (1.1), although some fourth order or eighth-order extensions of the presented second order scheme (2.1) have also been considered in the physical studies of the discrete problem (1.2) in .
As can be seen from above, the convergence order reduction problem is essentially due to the lack of enough regularity in the potential as well as in the solution of D-NLS (1.1). In some numerical recent work [33, 36] on nonlinear Schrödinger equations, the low-regularity of the solution and the order reduction problem of classical methods have also been addressed, where the roughness was introduced to the equation through the initial data. Here in D-NLS (1.1), the difference is that the roughness is induced inherently to the solution by the spatial random potential , even under smooth initial input as we shall discuss in details later. Therefore, numerically more special efforts are needed to take care of both the potential and the solution.
3. Low-regularity Fourier integrator
To raise the accuracy of temporal approximation, in this section, we will consider the low-regularity Fourier integrator (LRI) from  to solve the D-NLS (1.1). The LRI has been proposed originally to integrate the cubic nonlinearity in the equation with the second order accuracy by losing two spatial derivatives of the solution. Here, we need some particular efforts to handle the potential term in (1.1). We shall first derive the scheme and then analyze its convergence.
We shall follow some of the notations used in , and for simplicity, the spatial variable of the functions will be omitted, e.g. and , .
3.1. Derivation of scheme
By introducing the twisted variable,
and then by using the Duhamel’s formula, the mild of solution of (1.1) at for some and reads
By denoting for simplicity
we get from (3.2),
It is clear from the above that , where the term does not involve any spatial derivatives. Therefore, by denoting for short and making approximations and in the integrand of (3.4), where the latter clearly costs two spatial derivatives, we get the following expansion/approximation:
The constant in the above truncation term depends on and . Under the boundedness of these terms, globally, its optimal accuracy is at first order.
where denotes the and terms coming out from the cubic nonlinearity:
and denote the terms coming out from the potential part:
In (3.6), the truncation terms are those terms from the multiplications of the cubic nonlinearity and the potential term, where the highest order spatial derivatives evolved in the error constant are still and .
In the following, we shall compute or further approximate the three terms and in a sequel. The goal is to obtain their explicit (approximated) formulas in the physical space and meanwhile reach the global second order accuracy by losing at most two spatial derivatives of the functions including the solution and the potential
. Firstly, let us denote the Fourier transform of some functionas
and we define the operator as
Integrating potential terms.
For , by directly letting in (3.8a) and then integrating the rests exactly, it is easy to see that the following approximation
fulfills our aforementioned goal. The truncation error here is with two spatial derivatives involved in the error constant. Now the difficulty comes to compute . First of all, we write in the Fourier space
It is seen that the two frequencies and are coupled here, which prevents us from evaluating the integral explicitly in the physical space. So some further approximations are needed, otherwise the scheme will be defined in Fourier space and one will have to deal with the convolution. Noting that
we can approximate the integrand in (3.10) as
so that the truncation error here is with the loss of two spatial derivatives. It seems that this approximation fulfills our goal, but in the above formula (3.12) for , we note that the second term is proportional to , which corresponds to an explicit spatial derivative in the physical space. This kind of explicit spatial derivative term will induce stability issue in the final numerical scheme, and one will end up with a CFL constraint on the time step and spatial mesh size for computing. It is known that spatial discretizations such as the Fourier pseudo-spectral method will only have an algebraic convergence rate and the rate could be very low when the solution is less regular . Therefore, one will need to use small mesh size in space for accurate approximations, and then the CFL constraint would lead to severe inefficiency for practical computations in the rough solution situation. To enhance the stability, we introduce a filter :
and the approximation in (3.12) further becomes
where we used the fact . Similar filtering strategy for stability has been proposed in a very recent work  for the nonlinear Schrödinger equation, and we remark here that there are other possible choices of the filters  that could work here. Finding out the best one is beyond the scope of this work. Then in total, (3.10) is approximated as
Here the and can be regarded as two regularized versions of the rough potential , which can be pre-computed. The truncation error here in (3.14) is a combination of error from (3.11) and (3.13), which in total is with the loss of two spatial derivatives. Now the formula (3.14) is defined explicitly in the physical space, and it can be programmed efficiently under the fast Fourier transform (FFT).
Integrating the nonlinearity.
the nonlinear terms in can be written as
The last two terms of the above are of , and we can adopt the approximation similarly as before to get
For the part , in Fourier space we have
Then by the key fact for the phase term in the above :
one finds that
It can be seen that the above approximation is with the cost of two spatial derivatives. After exact integrations in the above, the explicit formulas of and can be obtained in physical space:
Here we refer the readers to  for the detailed calculations of the last two formulas above. With the findings (3.18) and (3.19), we get the total approximation of in physical space from (3.17) and (3.16) as
Plugging the above approximation for into (3.6), we get
As suggested in , the first term in the above approximation could be simplified as
which gives the desired second order approximation of for with the cost of two spatial derivatives.
The LRI scheme (3.20) is fully explicit in time, and it is explicitly defined in physical space. The two regularized potentials and in (3.15) can be computed once for all before the temporal iteration. The spatial discretization in (3.20) can be done by the Fourier pseudo-spectral method  under our periodic setup. In practical implementation, we take as the spatial mesh size with some integer as the total number of space grid points, and (3.20) can be programmed efficiently via FFT with the computational cost per time step .
3.2. Convergence analysis
Here we give a convergence analysis for the proposed LRI (3.20) under assumption that in (1.1) a sampled potential for some . For simplicity of notations, we shall denote for with some generic constant independent of or .
Assume that and for some , then for the numerical solution defined in the LRI (3.20) up to some fixed , there exists constant such that when ,
To prove the theorem, we first establish some lemmas for stability and local error estimates. First of all, we have some formala prior estimate results for the solution of (1.1) when .
(A prior estimate in ) If and for some in (1.1), then for some , we have .