Qinglin Tang

Weizhu Bao, Yongyong Cai, Xiaowei Jia et al.

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{τ^2}{\varepsilon^2}$ and $h^{m_0}+τ^2+\varepsilon^2$, where $h$ is the mesh size, $τ$ 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(τ)$ for all $\varepsilon\in(0,1]$ and optimally with quadratic convergence rate at $O(τ^2)$ in the regimes when either $\varepsilon=O(1)$ or $0<\varepsilon\lesssim τ$. 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^+$.

Weizhu Bao, Yongyong Cai, Xiaowei Jiao et al.

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 $τ$ 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: $τ=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 $τ=O(\varepsilon^2)$ when $0<\varepsilon\ll 1$. Extensive numerical results are reported to support our error estimates.

Weizhu Bao, Remi Carles, Chunmei Su et al.

We present a regularized finite difference method for the logarithmic Schrödinger equation (LogSE) and establish its error bound. Due to the blow-up of the logarithmic nonlinearity, i.e. $\ln ρ\to -\infty$ when $ρ\rightarrow 0^+$ with $ρ=|u|^2$ being the density and $u$ being the complex-valued wave function or order parameter, there are significant difficulties in designing numerical methods and establishing their error bounds for the LogSE. In order to suppress the round-off error and to avoid blow-up, a regularized logarithmic Schrödinger equation (RLogSE) is proposed with a small regularization parameter $0<\varepsilon\ll 1$ and linear convergence is established between the solutions of RLogSE and LogSE in term of $\varepsilon$. Then a semi-implicit finite difference method is presented for discretizing the RLogSE and error estimates are established in terms of the mesh size $h$ and time step $τ$ as well as the small regularization parameter $\varepsilon$. Finally numerical results are reported to confirm our error bounds.

Weizhu Bao, Rémi Carles, Chunmei Su et al.

We present and analyze two numerical methods for the logarithmic Schr{ö}dinger equation (LogSE) consisting of a regularized splitting method and a regularized conservative Crank-Nicolson finite difference method (CNFD). In order to avoid numerical blow-up and/or to suppress round-off error due to the logarithmic nonlinearity in the LogSE, a regularized logarithmic Schr{ö}dinger equation (RLogSE) with a small regularized parameter 0 < $ε$ $\ll$ 1 is adopted to approximate the LogSE with linear convergence rate O($ε$). Then we use the Lie-Trotter splitting integrator to solve the RLogSE and establish its error bound O($τ$ 1/2 ln($ε$ --1)) with $τ$ > 0 the time step, which implies an error bound at O($ε$ + $τ$ 1/2 ln($ε$ --1)) for the LogSE by the Lie-Trotter splitting method. In addition, the CNFD is also applied to discretize the RLogSE, which conserves the mass and energy in the discretized level. Numerical results are reported to confirm our error bounds and to demonstrate rich and complicated dynamics of the LogSE.

Xavier Antoine, Antoine Levitt, Qinglin Tang

We propose a preconditioned nonlinear conjugate gradient method coupled with a spectral spatial dis-cretization scheme for computing the ground states (GS) of rotating Bose-Einstein condensates (BEC), modeled by the Gross-Pitaevskii Equation (GPE). We first start by reviewing the classical gradient flow (also known as imaginary time (IMT)) method which considers the problem from the PDE standpoint, leading to numerically solve a dissipative equation. Based on this IMT equation, we analyze the forward Euler (FE), Crank-Nicolson (CN) and the classical backward Euler (BE) schemes for linear problems and recognize classical power iterations, allowing us to derive convergence rates. By considering the alternative point of view of minimization problems, we propose the preconditioned gradient (PG) and conjugate gradient (PCG) methods for the GS computation of the GPE. We investigate the choice of the preconditioner, which plays a key role in the acceleration of the convergence process. The performance of the new algorithms is tested in 1D, 2D and 3D. We conclude that the PCG method outperforms all the previous methods, most particularly for 2D and 3D fast rotating BECs, while being simple to implement.

Xavier Antoine, Qinglin Tang, Yong Zhang

In this paper, we propose some efficient and robust numerical methods to compute the ground states and dynamics of Fractional Schrödinger Equation (FSE) with a rotation term and nonlocal nonlinear interactions. In particular, a newly developed Gaussian-sum (GauSum) solver is used for the nonlocal interaction evaluation \cite{EMZ2015}. To compute the ground states, we integrate the preconditioned Krylov subspace pseudo-spectral method \cite{AD1} and the GauSum solver. For the dynamics simulation, using the rotating Lagrangian coordinates transform \cite{BMTZ2013}, we first reformulate the FSE into a new equation without rotation. Then, a time-splitting pseudo-spectral scheme incorporated with the GauSum solver is proposed to simulate the new FSE.

Zixu Feng, Qinglin Tang

In this article, we propose a unified framework for preconditioned Riemannian gradient (P-RG) methods to minimize Gross-Pitaevskii (GP) energy functionals with rotation on a Riemannian manifold. This framework enables comprehensive analysis of existing projected Sobolev gradient methods and facilitates the construction of highly efficient P-RG algorithms. Under mild assumptions on the preconditioner, we prove energy dissipation and global convergence. Local convergence is more challenging due to phase and rotational invariances. Assuming the GP functional is Morse-Bott, we derive a sharp Polyak-Łojasiewicz (PL) inequality near minimizers. This allows precise characterization of the local convergence rate via the condition number $μ/L$, where $μ$ and $L$ are the lower and upper bounds of the spectrum of a combined operator (preconditioner and Hessian) on a closed subspace. By combining spectral analysis with the PL inequality, we identify a quasi-optimal preconditioner achieving the best possible local convergence rate: $(L-μ)/(L+μ)+\varepsilon$ ($\varepsilon>0$ small). To our knowledge, this is the first rigorous derivation of the local convergence rate for P-RG methods applied to GP functionals with two symmetry structures. Numerical experiments on rapidly rotating Bose-Einstein condensates validate the theoretical results and compare the performance of different preconditioners.

Zixu Feng, Patrick Henning, Qinglin Tang

The structure and degeneracy of ground states of the Gross-Pitaevskii energy functional play a central role in both analysis and computation, yet a characterization of the ground-state manifold in the presence of symmetries remains a fundamental challenge. In this paper, we establish sharp results describing the geometric structure of local minimizers and its implications for optimization algorithms. We show that when local minimizers are non-unique, the Morse-Bott condition provides a natural and sufficient criterion under which the ground-state set partitions into finitely many embedded submanifolds, each coinciding with an orbit generated by the intrinsic symmetries of the energy functional, namely phase shifts and spatial rotations. This yields a structural characterization of the ground-state manifold purely in terms of these natural symmetries. Building on this geometric insight, we analyze the local convergence behavior of the preconditioned Riemannian gradient method (P-RG). Under the Morse-Bott condition, we derive the optimal local $Q$-linear convergence rate and prove that the condition holds if and only if the energy sequence generated by P-RG converges locally $Q$-linearly. In particular, on the ground-state set, the Morse-Bott condition is satisfied if and only if the minimizers decompose into finitely many symmetry orbits and the P-RG exhibits local linear convergence in a neighborhood of this set. When the condition fails, we establish a local sublinear convergence rate. Taken together, these results provide a precise picture: for the Gross-Pitaevskii minimization problem, the Morse-Bott condition acts as the exact threshold separating linear from sublinear convergence, while simultaneously determining the symmetry-induced structure of the ground-state manifold. Our analysis thus connects geometric structure, symmetry, and algorithmic performance in a unified framework.

Weizhu Bao, Shidong Jiang, Qinglin Tang et al.

We present efficient and accurate numerical methods for computing the ground state and dynamics of the nonlinear Schrödinger equation (NLSE) with nonlocal interactions based on a fast and accurate evaluation of the long-range interactions via the nonuniform fast Fourier transform (NUFFT). We begin with a review of the fast and accurate NUFFT based method in \cite{JGB} for nonlocal interactions where the singularity of the Fourier symbol of the interaction kernel at the origin can be canceled by switching to spherical or polar coordinates. We then extend the method to compute other nonlocal interactions whose Fourier symbols have stronger singularity at the origin that cannot be canceled by the coordinate transform. Many of these interactions do not decay at infinity in the physical space, which adds another layer of complexity since it is more difficult to impose the correct artificial boundary conditions for the truncated bounded computational domain. The performance of our method against other existing methods is illustrated numerically, with particular attention on the effect of the size of the computational domain in the physical space. Finally, to study the ground state and dynamics of the NLSE, we propose efficient and accurate numerical methods by combining the NUFFT method for potential evaluation with the normalized gradient flow using backward Euler Fourier pseudospectral discretization and time-splitting Fourier pseudospectral method, respectively. Extensive numerical comparisons are carried out between these methods and other existing methods for computing the ground state and dynamics of the NLSE with various nonlocal interactions. Numerical results show that our scheme performs much better than those existing methods in terms of both accuracy and efficiency.