Efficient numerical methods for computing ground states and dynamics of dipolar Bose-Einstein condensates

New efficient and accurate numerical methods are proposed to compute ground states and dynamics of dipolar Bose-Einstein condensates (BECs) described by a three-dimensional (3D) Gross-Pitaevskii equation (GPE) with a dipolar interaction potential. Due to the high singularity in the dipolar interaction potential, it brings significant difficulties in mathematical analysis and numerical simulations of dipolar BECs. In this paper, by decoupling the two-body dipolar interaction potential into short-range (or local) and long-range interactions (or repulsive and attractive interactions), the GPE for dipolar BECs is reformulated as a Gross-Pitaevskii-Poisson type system. Based on this new mathematical formulation, we prove rigorously existence and uniqueness as well as nonexistence of the ground states, and discuss the existence of global weak solution and finite time blowup of the dynamics in different parameter regimes of dipolar BECs. In addition, a backward Euler sine pseudospectral method is presented for computing the ground states and a time-splitting sine pseudospectral method is proposed for computing the dynamics of dipolar BECs. Due to the adaption of new mathematical formulation, our new numerical methods avoid evaluating integrals with high singularity and thus they are more efficient and accurate than those numerical methods currently used in the literatures for solving the problem. Extensive numerical examples in 3D are reported to demonstrate the efficiency and accuracy of our new numerical methods for computing the ground states and dynamics of dipolar BECs

Error estimates of numerical methods for the nonlinear Dirac equation in the nonrelativistic limit regime

We present several numerical methods and establish their error estimates for the discretization of the nonlinear Dirac equation in the nonrelativistic limit regime, involving a small dimensionless parameter $0<\varepsilon\ll 1$ which is inversely proportional to the speed of light. In this limit regime, the solution is highly oscillatory in time, i.e. there are propagating waves with wavelength $O(\varepsilon^2)$ and $O(1)$ in time and space, respectively. We begin with the conservative Crank-Nicolson finite difference (CNFD) method and establish rigorously its error estimate which depends explicitly on the mesh size $h$ and time step $\tau$ as well as the small parameter $0<\varepsilon\le 1$. Based on the error bound, in order to obtain `correct' numerical solutions in the nonrelativistic limit regime, i.e. $0<\varepsilon\ll 1$, the CNFD method requests the $\varepsilon$-scalability: $\tau=O(\varepsilon^3)$ and $h=O(\sqrt{\varepsilon})$. Then we propose and analyze two numerical methods for the discretization of the nonlinear Dirac equation by using the Fourier spectral discretization for spatial derivatives combined with the exponential wave integrator and time-splitting technique for temporal derivatives, respectively. Rigorous error bounds for the two numerical methods show that their $\varepsilon$-scalability is improved to $\tau=O(\varepsilon^2)$ and $h=O(1)$ when $0<\varepsilon\ll 1$ compared with the CNFD method. Extensive numerical results are reported to confirm our error estimates.Comment: 35 pages. 1 figure. arXiv admin note: substantial text overlap with arXiv:1504.0288

A uniformly accurate (UA) multiscale time integrator pseudospectral method for the Dirac equation in the nonrelativistic limit regime

We propose and rigourously analyze a multiscale time integrator Fourier pseudospectral (MTI-FP) method for the Dirac equation with a dimensionless parameter $\varepsilon\in(0,1]$ which is inversely proportional to the speed of light. In the nonrelativistic limit regime, i.e. $0<\varepsilon\ll 1$, the solution exhibits highly oscillatory propagating waves with wavelength $O(\varepsilon^2)$ and $O(1)$ in time and space, respectively. Due to the rapid temporal oscillation, it is quite challenging in designing and analyzing numerical methods with uniform error bounds in $\varepsilon\in(0,1]$. We present the MTI-FP method based on properly adopting a multiscale decomposition of the solution of the Dirac equation and applying the exponential wave integrator with appropriate numerical quadratures. By a careful study of the error propagation and using the energy method, we establish two independent error estimates via two different mathematical approaches as $h^{m_0}+\frac{\tau^2}{\varepsilon^2}$ and $h^{m_0}+\tau^2+\varepsilon^2$, where $h$ is the mesh size, $\tau$ is the time step and $m_0$ depends on the regularity of the solution. These two error bounds immediately imply that the MTI-FP method converges uniformly and optimally in space with exponential convergence rate if the solution is smooth, and uniformly in time with linear convergence rate at $O(\tau)$ for all $\varepsilon\in(0,1]$ and optimally with quadratic convergence rate at $O(\tau^2)$ in the regimes when either $\varepsilon=O(1)$ or $0<\varepsilon\lesssim \tau$. Numerical results are reported to demonstrate that our error estimates are optimal and sharp. Finally, the MTI-FP method is applied to study numerically the convergence rates of the solution of the Dirac equation to those of its limiting models when $\varepsilon\to0^+$.Comment: 25 pages, 1 figur

Numerical methods and comparison for the Dirac equation in the nonrelativistic limit regime

We analyze rigorously error estimates and compare numerically spatial/temporal resolution of various numerical methods for the discretization of the Dirac equation in the nonrelativistic limit regime, involving a small dimensionless parameter $0<\varepsilon\ll 1$ which is inversely proportional to the speed of light. In this limit regime, the solution is highly oscillatory in time, i.e. there are propagating waves with wavelength $O(\varepsilon^2)$ and $O(1)$ in time and space, respectively. We begin with several frequently used finite difference time domain (FDTD) methods and obtain rigorously their error estimates in the nonrelativistic limit regime by paying particular attention to how error bounds depend explicitly on mesh size $h$ and time step $\tau$ as well as the small parameter $\varepsilon$. Based on the error bounds, in order to obtain `correct' numerical solutions in the nonrelativistic limit regime, i.e. $0<\varepsilon\ll 1$, the FDTD methods share the same $\varepsilon$-scalability on time step: $\tau=O(\varepsilon^3)$. Then we propose and analyze two numerical methods for the discretization of the Dirac equation by using the Fourier spectral discretization for spatial derivatives combined with the exponential wave integrator and time-splitting technique for temporal derivatives, respectively. Rigorous error bounds for the two numerical methods show that their $\varepsilon$-scalability on time step is improved to $\tau=O(\varepsilon^2)$ when $0<\varepsilon\ll 1$. Extensive numerical results are reported to support our error estimates.Comment: 34 pages, 2 figures. arXiv admin note: substantial text overlap with arXiv:1511.0119

Vortex patterns and the critical rotational frequency in rotating dipolar Bose-Einstein condensates

Based on the two-dimensional mean-field equations for pancake-shaped dipolar Bose-Einstein condensates in a rotating frame with both attractive and repulsive dipole-dipole interaction (DDI) as well as arbitrary polarization angle, we study the profiles of the single vortex state and show how the critical rotational frequency change with the s-wave contact interaction strengths, DDI strengths and the polarization angles. In addition, we find numerically that at the `magic angle' $\vartheta=\arccos(\sqrt{3}/3)$, the critical rotational frequency is almost independent of the DDI strength. By numerically solving the dipolar GPE at high rotational speed, we identify different patterns of vortex lattices which strongly depend on the polarization direction. As a result, we undergo a study of vortex lattice structures for the whole regime of polarization direction and find evidence that the vortex lattice orientation tends to be aligned with the direction of the dipoles