Conservative modified Crank–Nicolson and time-splitting wavelet methods for modeling Bose–Einstein condensates in delta potentials

Abstract

This paper explores two wavelet-based energy-conserving algorithms for the Gross–Pitaevskii equation with delta potentials in Bose–Einstein condensates, named modified Crank–Nicolson wavelet method and time-splitting wavelet method, respectively. Both proposed methods can preserve the intrinsic properties of original problems as much as possible. Meanwhile, the rigorous error estimates and some conservative properties are investigated. They are proved to preserve the charge conservation exactly. The global energy conservation laws can be preserved under several conditions. In practical computations, to avoid a large drift in energy values caused by discontinuous potential well, an improved discrete delta function is implemented. Numerical experiments for attractive and repulsive cases are conducted during long time computations to show the performances of the proposed methods and verify the theoretical analysis.

Introduction

The nonlinear Schrödinger equation is an important model in many branches of physics and applied mathematics, such as nonlinear quantum field theory, condensed matter, nonlinear optics, hydrodynamics, self-focusing in laser pulse, thermodynamic process in meso scale systems, plasma and so on [1], [2], [3], [4]. We focus on Bose–Einstein condensate (BEC), which was proposed theoretically in 1924–1925 and achieved experimentally in 1995. Undoubtedly, due to the remarkable discovery of BEC, it gives us approach and methodology to explore the extremely complex quantum world.

In this paper, we consider the following nonlinear Schrödinger equation (NLSE) with cubic nonlinearity, known as the Gross–Pitaevskii equation (GPE) [5]:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill i\hslash {\partial }_t\Phi (\mathbf{r},t)=(-\frac{{\hslash }^2}{2m}{\nabla }^2+V_{\mathrm{ext}}\left(\mathbf{r}\right)+g{\left|\Phi (\mathbf{r},t)\right|}^2)\Phi (\mathbf{r},t).\end{array}\end{array}$$

Here, $\Phi (\mathbf{r},t)$ is a complex time dependent function, V_ext is an external potential, and the coupling constant $g=4\pi {\hslash }^{2}a/m,$ which is related to the scattering length a. Eq. (1) requires the s-wave scattering length should be much smaller than the average distance between atoms, and the number of atoms in the condensate should be much larger than 1. Thus, mean field approximation (1) is only valid for dilute boson gases (or usually called weakly interacting boson gases) at extremely low temperature such that almost all particles are in the same states.

The NLSE or GPE (1) admits following conserved quantities:

(1) The total number: the boson gas is only investigated non-relativistically, and no real particles can be created or annihilated, so the total number of atoms

$$\mathcal{N}={\int \left|\Phi (\mathbf{r},t)\right|}^{2}d\mathbf{r}$$

should be conserved in the dynamics (with Φ normalized by $\mathcal{N}$);

(2) The global energy: The gas system is an isolated system, so the total energy

$$\mathcal{E}\left(\Phi \right)=\int [\frac{{\hslash }^2}{2m}{\left|\nabla \Phi (\mathbf{r},t)\right|}^2+V_{\mathrm{ext}}{\left(\mathbf{r}\right)\left|\Phi (\mathbf{r},t)\right|}^2+\frac{g}{2}{\left|\Phi (\mathbf{r},t)\right|}^4]d\mathbf{r}$$

should also be conserved.

Concerning the NLSE or GPE (1) at some finite time, some solitary and stationary solutions were obtained [5]. Since available theoretical solutions for NLSE or GPE, especially for higher-dimensional cases, are still very limited, investigating numerically is an important tool to understand physical behavior of this system. There are many numerical methods for modeling this system [6], [7], [8]. Especially, Delfour et al. [9] proposed a general finite difference method; Daǧ [10] proposed a B-spline finite element method; Chang et al. [11] discussed several different schemes such as Hospscotch scheme, split-step Fourier scheme, pseudo-spectral scheme for generalized NLSE; Bao et al. [12] proposed a time-splitting Laguerre–Hermite pseudo-spectral method and a Fourier spectral method; Compact finite difference schemes for one-dimensional case are constructed in [13]. In addition, numerous symplectic and multi-symplectic methods have been constructed to simulate Schrödinger system and other Hamiltonian equations due to their long time simulation property and good preservation property of conservative quantities of original problem. Theoretically, all most real physical process with negligible dissipation can be cast in suitable Hamiltonian formulation in phase space with symplectic structure, which means the symplectic-preserving algorithms can hold these properties spontaneously. A bunch of works have verified its effectiveness during long-time simulations. Therefore, symplectic integrator is the most popular approach to structure-preservation. Recently, Chen et al. [14] proposed several symplectic and multi-symplectic methods for NLSE. The multi-symplectic Runge–Kutta and Fourier spectral methods were employed to solve the fourth-order Schrödinger equations with trapped term by Hong el at. in [15]. Chen et al. [16], [17] took the splitting technique into the multi-symplectic integrator for solving coupled NLSE. Cai el at. [18], [19] proposed several local structure-preserving schemes for coupled NLSE. But the above methods are almost to solve NLSE or GPE with $V_{\mathrm{ext}}=0$ or continuous potentials. For the system with discontinuous potential well, the common way is to modify the finite difference methods to cope with the general cases. Unfortunately, it is difficult for the finite difference methods to meet requirements of high accuracy and high resolution of complex physical process, since they have not spectral accuracy and structure-preserving properties.

Our main aim is to construct efficient and partially conservative method for the NLSE or GPE with delta potentials without considerations of symplectic and multi-symplectic formulations. Especially in this paper, we consider

$$i\hslash {\partial }_t\Phi (x,t)=-\frac{{\hslash }^2}{2m}{\partial }_{xx}\Phi (x,t)+V_{\mathrm{ext}}\left(x\right)\Phi (x,t)+g{\left|\Phi (x,t)\right|}^2\Phi (x,t),$$

with a single delta potential for modeling a short range interaction. Due to the discontinuity of delta potential, the traditional numerical discrete framework, e.g. finite difference method, symplectic/multi-symplectic method, discontinuous Galerkin method, can not be implemented directly. Therefore, based on the wavelet collocation scheme [20], [21], we propose a simple modified Crank–Nicolson wavelet method, and the analogous time-splitting wavelet method. Wavelet collocation method makes the corresponding spatial differentiation matrix sparse and demands less computations. For the discontinuous potential function, we apply a technique of discrete delta function to the original schemes, which makes them adjust to the discontinuous problems. Meanwhile, the convergence property is discussed. It is proved to preserve the charge conservation exactly. The global energy conservation laws also can be preserved under several conditions.

The paper is organized in the following way. In Section 2, wavelet collocation scheme and modified Crank–Nicolson wavelet method for GPE are introduced. In Section 3, the analogous time-splitting wavelet method is proposed. In Section 4, some theoretical analysis, such as convergence and conservative properties, are presented. The technique of discrete delta function and numerical experiments for attractive and repulsive nonlinearity are presented in Section 5, which show the effectiveness of the proposed algorithms. Finally, conclusions are made in Section 6.