Numerical solution to the time-dependent Gross-Pitaevskii equation

In this work we employ the split-step technique combined with a Legendre pseudospectral representation to solve various time-dependent Gross-Pitaevskii equations (GPE). Our findings based on the numerical accuracy of this approach applied for one-dimensional (1D) and two-dimensional (2D) problems show that it can provide accurate and stable solutions. Moreover, this approach has been applied to study the dynamics of the Bose-Einstein condensate which is modeled with the GPE. The breathing of condensate with an repulsive and attractive interactions trapped in 1D and 2D harmonic potentials has been simulated as well.


I. INTRODUCTION
Gross-Pitaevskii theory [1,2] is one of the successful mathematical models describing the Bose-Einstein condensate (BEC) in which the bosons simultaneously occupy the same quantum state of the lowest energy. The BEC was theoretically predicted by Bose and Einstein [3,4] in 1925, however, it had first been experimentally obtained in 1995 [5][6][7].
Various methods have been tested to obtain numerical solution for the time-dependent Gross-Pitaevskii equation (GPE). Weideman [8] employed a split-step Fourier pseudospectral method. Adhikari [9] used a split-step Crank-Nicolson approach for solving the GPE.
Moreover, Bao [10] and Wang [11,12] have made a lot of contributions for developing numerical techniques for solving a single and coupled GPEs with different trapping potentials.
In this work our goal is to solve the time-dependent GPE using the split-step technique combined with a Legendre pseudospectral representation. To our knowledge, the split-step Legendre pseudospectral approach has not been applied to this problem in a systematic and detailed way, yet. In our calculation, we will initially test a numerical accuracy of the proposed approach for a nonlinear one-dimensional (1D) and two-dimensional (2D) Schrödinger equation for which an analytical solution is known. Then we will apply the approach for solving the time-dependent GPE which simulates the dynamics of the BEC with a repulsive and an attractive interactions trapped in 1D and 2D harmonic potentials at zero temperature. The paper is organized as follows. In Section 2, we present the formalism of the Gross-Pitaevskii theory for the harmonic trap. In Section 3 we give a brief discussion of a split-step technique and a pseudospectral method for the 1D and 2D problems. In Section 4 we present the numerical results and discussions for test problems and the dynamics of condensate trapped in 1D and 2D harmonic potentials. Then a conclusion follows.

II. GROSS-PITAEVSKII THEORY FOR TRAPPED BOSONS
The dynamics of boson particles confined in the trap (potential) can be modeled with the following time-dependent Gross-Pitaevskii equation [1,2,[10][11][12]: Here ψ(r, t) is the BEC wave function, m is the mass of boson, V is an external confining potential (trap), a s is the s-wave scattering length and N is the number of bosons in the condensate defined as N ≡ ψ(r, t)|ψ(r, t) , and the energy per particle is defined as Usually the external trap potential V is the harmonic trapping potential mω 2 In the trapping potential ω x , ω y and ω z are the frequencies along the x-, y-and z-directions, respectively.
It is convenient to use the equation (1) in the dimensionless form in the calculation. By choosing thet = ω x t and a x = ( /mω x ) 1/2 as the dimensionless time and length units, respectively, we can have the rescaled quantities:r = r/a x ,ψ(r,t) = a 3/2 x /N 1/2 ψ(r, t), andẼ = E/ ω x . After removing all˜, the time-dependent GPE (1) can be written in the dimensionless form [10][11][12]: where the nonlinear constant g = 4πa s N/a x . In this work we use the 1D and 2D harmonic trapping potentials given in the Cartesian coordinates. In 1D, r = (x), and the trapping For the 2D case, r = (x, y) and V (x, y) = x 2 2 + λ 2 y 2 2 , with λ = ωy ωx . We note that the equation (3) has no analytic solution, however, it can be solved employing the numerical techniques which I will discuss in the next section.

III. NUMERICAL PROCEDURE
In this work we consider consider the initial-boundary value problem for the GPE given in equation (3). We can write it in the form: with an initial condition ψ(r, 0) = ψ 0 (r), r ∈ R d (d = 1, 2, 3), and a boundary condition lim |r|→±∞ ψ(r, t) = 0, t ≥ 0. Here d means a dimension of space. A formal solution for this equation can be written as ψ(r, t + ∆t) = e −iH∆t ψ(r, t).
In our numerical calculation we solve the equation (4) for two cases: (i) without a splitstep technique and (ii) with a split-step technique. In both cases when we solve the equation For the first case in which a split-step technique is not used a solution for equation (4) is indeed given in the form (5), and can be found with the CN scheme: where ψ n = ψ(r, t n ) and ∆t > 0 is a time step size with which a time step is given as t n = n∆t, n = 0, 1, . . . Now let's consider the second case in which the split-step technique is employed. The main idea of this approach is that we divide the Hamiltonian operator H into two parts: a linear part A = − 1 2 ∇ 2 r and a non-linear part B = V + g|ψ(r, t)| 2 . Then we solve the equation (4) in terms of two successive steps. In the Step 1, we solve the following a nonlinear operator for the time step ∆t; then Step 2, followed by solving for the same time step.
From the nonlinear subproblem (7) in the Step 1, for t n ≤ t ≤ t n+1 we can have an exact solution in the form In Step 2, we obtain a numerical solution in the following form: which is an exactly same approach as given in (6). Then the solution (5) for the equation can be found in the form: We note that this splitting technique has only first-order in time (equation (11)), so that it can easily be implemented. Note that one may use higher-order splitting approaches [11][12][13] to solve the GPE as well.
For numerical calculation of ∇ 2 r in H and A operators, we used the Legendre differentiation matrix of the second order (d 2 ij ) (details can be found in [14,15]). For 1D case, where I is the unit matrix and ⊗ denotes the Kronecker product [14,15]. We note that this Legendre pseudospectral approach had been employed to obtain a solution for the time-independent GPE for an anisotropic 3D trapping potential in our previous work [15].
Below we have shown the solutions (6) and (11) in the discrete forms which can directly be implemented in coding. For the 1D GPE, the solution (6) can be written as The solution (11) can be obtained as In both (12) and (13) equations, we use the notations: ψ n i ≡ ψ(x i , t n ), I is unit matrix, "diag" means a diagonal matrix, and x i are non-uniformly distributed N x + 1 number of the Legendre-Gauss-Lobatto grid points in an interval [a, b]. Note that in numerical calculation we need to take into account the boundary condition in such a way that two end points of a column vector will not be used for the Dirichlet boundary condition, and for the 2D case in discretized form a main difference comes with the Kronecker product as mentioned above [15]. Example 1. First we consider the following GPE in 1D [11] with the homogeneous Dirichlet condition over [−20, 20]. The analytical solution for this . This is the usual non-linear Schrödinger equation since it does not include the trapping potential V (GPE (1)). We computed the maximum error between the exact solution and the approximated solution which are obtained from the LS and the SSLS methods. Table 1 shows the maximum of |ψ exact (x, t) − ψ approx (x, t)| until t = 4. A number of grid points along x-axis is N x = 128 and the time step size ∆t = 0.01. We note that this problem was also considered using the split-step finite-difference (SSFD) and split-step Fourier spectral (SSFS) methods in Ref. [11], and at t = 4 with ∆t = 0.01 value of this error obtained by the SSFD approach is 3.6620 · 10 −2 , while the SSFS method's estimation for it is 3.663 · 10 −2 ( Table 1 in [11]). In Ref. points. Then a number of grid point used in [11] is M = 4000. This is quite large number compared with the number we used (N x = 128), indeed. We say that a reason why our calculation with relatively few number of grid points works well is that the Legendre-Gauss-Lobatto grid points are distributed non-uniformly in the interval, that is, more points come near two ends of the interval [16].
where V (x, y) = 1 − sin 2 x sin 2 y. An initial condition ψ 0 (x, y) = sin x sin y are used. The exact solution for this problem with the Dirichlet boundary condition is ψ exact (x, y, t) = sin x sin ye −i2t . Table 2 presents the maximum of |ψ exact (x, y, t) − ψ approx (x, y, t)| until time t = 30 using the LS and SSLS methods. The number of grid points along the x-and y-axis is N x = N y = 16 and the time step size ∆t = 0.01. This example problem was also discussed in Ref. [11] where an author used the SSFD and SSFS methods. The maximum errors from the SSFD and SSFS methods with N x = N y = 128 at t = 20 with ∆t = 0.01 are 4.057 · 10 −3 and 1.138·10 −10 , respectively [11]. From the errors estimated, we see that the SSFS approach has shown an extremely good result since its solution can be obtained with an almost analytical expression in terms of the ordinary Fourier pseudospectral representation.  [11]. We also note that more accurate results from the LS and SSLS approaches, especially for large computational domains of time and space can be easily obtained by handling the computational parameters, such as, increasing number of grid points.

B. Numerical applications
In this subsection based on numerical experiences obtained in a preceding subsection we will solve the GPE (3) of which an analytical solution is unknown. We will apply the split-step technique for the following examples which have a trap potential.

Example 3. We solve the GPE (3) in 1D
with an initial condition ψ 0 (x) = π −1/4 e −x 2 /2 , which is a ground state for the trap with g = 0. When g > 0, the interaction between particles is repulsive. In contrast, when g < 0, the interaction is attractive. We solve this problem on [−10, 10] with a Dirichlet boundary condition.
In Figure 1a we have shown the time-evolution of the condensate |ψ(x, t)| 2 for a repulsive interaction (g = 5). The dynamics of a ground state (an initial state) that is left to evolve with a different trapping potential is called a "breathing" [17,18]. The "breathing" of the condensate shown in panel 1a may happen because the initial trapping frequency is reduced, resulting in expansion of the trapping potential [17]. Panel 1b of Figure 1 presents how the nonlinear constant g affects the breathing frequency of the condensate. It's been clearly shown that the density at the center of the trap is lower for a larger value of g, which is expected because of the corresponding higher inter-particle repulsion. Computational parameters: N x = 128, ∆t = 0.001. Figure 2a presents the dynamics of the condensate for which an inter-particle interaction is attractive (g = −2.5). Because of this attractive behavior, density at the center of the trap |ψ(0, t)| 2 is higher for a smaller value of g (Panel 2b). We also observe that the oscillation frequency of the condensate in the trap increases with a smaller value of g since more attractive potential energy is added. In the above two calculations, we chose the number of grid points along x-axis N x = 128 and the time step size ∆t = 0.001. For the calculations the initial condition ψ 0 (x) is used, therefore all curves shown in panels 1b and 2b have the same starting point. We note that all calculations are independent on the numerical parameters, since a calculated discrete error norm ||ψ Nx=128 (t) − ψ Nx=64 (t)|| l 2 at time t = 12 is 10 −6 .
Example 4. We solve the GPE (3) in anisotropic 2D trap The initial condition for this problem ψ 0 (x, y) is taken as the ground state for an isotropic trap (λ 2 = 1). Then we let this initial state ψ 0 evolve in an anisotropic trap with λ 2 = 4.
Panel 3b demonstrates the breathing frequency of this condensate depending on values of the nonlinear constant g. When g increases, the initial density |ψ(x, y, 0)| 2 at time t = 0 lowers due to a repulsive interparticle interaction. However, it has been seen that the nature of all three curves is different from what has been appeared in Figure 2b. This may happen because the trap is contracted along the y-axis and particles could move higher and closer. Figure 4 shows results for the system composed of the particle with a negative scattering length (a s < 0) which makes the system attractive (g < 0). Panel 4a presents the density |ψ(x, y, t)| 2 for the anisotropic trap g = −2. As expected, similar patterns appear in Figure   4b like what was seen in Figure 2b since the anisotropic trap may support more particles in the center of the trap. We note that when an interaction becomes stronger (as g gets larger/smaller value), uneven (up and down) nature appears (red curve in panels 3b and 4b). We also note that results for our 2D harmonic trap calculations are stable for the grid size parameters and the time interval we employed.

V. CONCLUSION
In this paper we have solved the time-dependent Gross-Pitaevskii equation using splitstep technique combined with the Legendre-pseudospectral method. From our numerical experiment for a non-linear Schrödinger equation it's been shown that the split-step procedure is more accurate than an ordinary pseudospectral method. It has also been shown that the split-step Legendre-pseudospectral approach can be competitive one against the splitstep Fourier pseudospectral approach due to the it's non-uniform grid point distribution along with the requirement of relatively few grid points. Furthermore, the SSLS approach has been tested for the one-dimensional and two-dimensional time-dependent GPE with harmonic trap to simulate the dynamics of the Bose-Einstein condensate with repulsive and attractive interactions. Making the expansion of the 1D harmonic trapping potential, the breathing modes of the condensate have been revealed as a function of the time for various values of the nonlinear constant of the condensate. Moreover, the simulation of those modes has been discussed for 2D harmonic trap which contracts from a looser one to a more confining trap.