Source-linked AI summary
Fortran programs for the time-dependent Gross-Pitaevskii equation in a fully anisotropic trap
P. Muruganandam, S. K. Adhikari
TL;DR
The paper develops accessible numerical methods for stationary and non-stationary solutions of the nonlinear GP equation across harmonic-trap geometries. It combines real- and imaginary-time propagation with split-step Crank-Nicolson discretization, supplies Fortran programs, and reports stationary properties and oscillations. The programs agree with prior stationary calculations, while large nonlinear changes can destabilize 3D oscillations.
Problem
Accurate numerical solution of the nonlinear GP equation is needed across diverse trap symmetries describing dilute trapped Bose-Einstein condensates.
Method
The paper uses split-step Crank-Nicolson discretization with real-time propagation for dynamics and imaginary-time propagation for stationary states, implemented across multiple trap geometries.
Results
Stationary chemical potentials and rms sizes agree with previous calculations across several trap cases, while real-time programs produce non-stationary oscillations for different symmetries.
Takeaways & Limitations
The supplied programs provide accessible working implementations for computing stationary condensate properties and studying time-dependent oscillations in six trap symmetries.
Takeaways & Limitations
The programs target the standard cubic GP equation in harmonic potentials, and large nonlinearity changes can destabilize 3D oscillations; systematic oscillation studies are not presented.
Abstract
from arXiv · showhide
We develop simple numerical algorithms for both stationary and non-stationary solutions of the time-dependent Gross-Pitaevskii (GP) equation describing the properties of Bose-Einstein condensates at ultra low temperatures. In particular, we consider algorithms involving real and imaginary-time propagation based on a split-step Crank-Nicolson method. In a one-space-variable form of the GP equation we consider the one-dimensional linear, two-dimensional circularly symmetric, and the three-dimensional spherically-symmetric traps. In the two-space-variable form we consider the GP equation in two-dimensional anisotropic and three-dimensional axially-symmetric traps. The fully-anisotropic three-dimensional GP equation is also considered. Numerical results for the chemical potential and root-mean-square size of stationary states are reported using imaginary-time propagation programs for all the cases and compared with previously obtained results. Also presented are numerical results of non-stationary oscillation for different trap symmetries using real-time propagation programs. A set of convenient working codes developed in Fortran 77 are also provided for all these cases (twelve programs in all). In the case of two or three space variables, {Fortran 90/95 versions provide some simplification over the Fortran 77 programs}, and these programs are also included (six programs in all).
1. Introduction
The paper addresses numerical solution of the nonlinear, time-dependent GP equation across multiple harmonic-trap symmetries. It develops accessible real- and imaginary-time propagation algorithms with Crank-Nicolson discretization and supplies Fortran implementations.
- The nonlinear interaction term complicates accurate numerical solution of the GP equation for condensates in diverse harmonic traps.
- The GP equation combines first-order time and second-order space derivatives with harmonic and nonlinear potential terms.
- The split-step algorithm separately propagates spatial derivatives and then incorporates harmonic and nonlinear potential effects.
- Earlier numerical studies largely focused on stationary states, while relatively few applied Crank-Nicolson methods to time-dependent GP solutions.
- The paper develops a simple, efficient propagation algorithm, provides Fortran 77 and Fortran 90/95 programs, and compares stationary properties across trap symmetries.
2. Nonlinear Gross-Pitaevskii Equation
The paper formulates dimensionless GP equations for spherical, circular, axial, anisotropic, and reduced-dimensional harmonic traps. Coordinate transformations, normalization, interaction scaling, and boundary conditions are specified for numerical treatment.
- The GP equation describes a zero-temperature condensate wave function using mass, atom number, interaction strength, and scattering length parameters.
- Spherically-symmetric GP equation in 3D: For spherical traps, the radial transformation ϕ = rψ removes the first spatial derivative and yields simple vanishing boundary conditions at the origin and infinity.
- Dimensionless formulations use oscillator-length and time scalings, with interaction coefficients ℵ = 4πNa/l or ℵ = 8πNa/l depending on convention.
- Anisotropic GP equation in 3D: The fully anisotropic 3D trap assigns separate frequency ratios and coordinates to the x, y, and z directions, with independent boundary conditions at spatial infinity.
- Axially-symmetric GP equation in 3D: Axial symmetry reduces the problem to radial coordinate ρ and axial coordinate z, with a derivative boundary condition at ρ = 0.
- Quasi-one-dimensional and two-dimensional reductions assume confinement to transverse or axial ground states in elongated or disk-shaped traps.
- Circular symmetry further reduces the two-dimensional equations to one radial variable with zero derivative at the origin and vanishing behavior at infinity.
- The section presents the GP equations for different trap symmetries before extending the Crank-Nicolson method to them.
3. Split-Step Crank-Nicolson Method for the GP Equation in one Space Variable
The method splits GP time evolution into non-derivative and derivative parts, treating them successively with exact small-step propagation and semi-implicit Crank-Nicolson discretization. Imaginary-time propagation renormalizes iterates to obtain stationary ground states, while real-time propagation preserves normalization for dynamical studies.
- Real-time propagation: The Hamiltonian is split into H1 and H2, with each part propagated successively over a small time step.The H1 operation produces an intermediate wave function, and Crank-Nicolson propagation with H2 produces the next-time solution.
- Real-time propagation: Real-time split-step propagation conserves wave-function normalization, while its splitting error is proportional to ∆2.The method also preserves the symplectic structure and treats much of the Hamiltonian without combining it with the delicate derivative term.
- Crank-Nicolson discretization: The spatial derivative is discretized on a one-dimensional grid using a three-point formula averaged over present and future time levels.The resulting equations are tridiagonal and are solved through forward and backward recursions across the lattice.
- Crank-Nicolson discretization: The Crank-Nicolson scheme is second order in space and time, unconditionally stable, and preserves the boundary condition at infinity for small ∆/h2.These properties support accurate propagation when the space and time steps are chosen sufficiently small.
- Imaginary-time propagation: Imaginary-time propagation exponentially suppresses higher-energy eigenstates, allowing convergence toward the stationary ground state.Because it does not preserve normalization, the wave function must be renormalized during iteration; the resulting real-variable algorithm is robust and rapidly convergent.
- Circularly-symmetric case: The circularly-symmetric case adds a first spatial derivative and uses a nonzero-origin boundary condition, dϕ(r; t)/dr|r=0 = 0.At large radius, the wave function satisfies limr→∞ϕ(r; t) = 0.
4. Split-Step Crank-Nicolson method in two and three space variables
For two- and three-space-variable GP equations, the Hamiltonian is decomposed into three or four pieces and propagated through independent successive substeps. The same split-step strategy supports anisotropic two-dimensional, axially symmetric, and fully anisotropic three-dimensional traps.
- Anisotropic 2D GP equation: In the anisotropic two-dimensional case, the Hamiltonian is decomposed as H = H1 + H2 + H3.The three pieces are propagated successively in independent substeps to complete one evolution over time step ∆.
- Anisotropic 2D GP equation: The H1, H2, and H3 propagations use the previously defined non-derivative and Crank-Nicolson operations.The chemical potential is then expressed for stationary states using the corresponding two-dimensional wave function.
- Stationary states: The chemical-potential expressions for the multidimensional stationary states are simplified using integrations by parts.This is stated for both the two-dimensional and three-dimensional formulations.
- Axially symmetric 3D GP equation: In the axially symmetric three-dimensional case, the GP Hamiltonian is likewise divided into H1, H2, and H3.The same three-part propagation policy is applied, with the radial and axial derivative pieces treated using the corresponding Crank-Nicolson schemes.
- Anisotropic 3D GP equation: For the fully anisotropic three-dimensional equation, the Hamiltonian is decomposed as H = H1 + H2 + H3 + H4.Four successive substeps complete one full evolution step, with H1 treated as before and H2–H4 treated by the derivative propagation schemes.
5. Description of Numerical Programs
The paper supplies Fortran programs for imaginary- and real-time GP propagation across one-, two-, and three-dimensional trap geometries. The codes calculate stationary and dynamical observables, with user-controlled grids, steps, nonlinearities, and boundary conditions.
- One-space-variable programs: Six one-variable programs cover 1D, circularly symmetric 2D, and radially symmetric 3D equations using both imaginary- and real-time propagation.The corresponding programs are imagtime1d.F, imagtimecir.F, imagtimesph.F, realtime1d.F, realtimecir.F, and realtimesph.F.
- Program outputs: The programs calculate wave functions, chemical potentials, energies, rms sizes or radii, normalization, central wave-function values, and nonlinearities.The principal program variables expose mesh, potential, propagation, normalization, and observable parameters to the user.
- Grid and boundaries: The spatial grid must extend far enough that the wave function reaches its asymptotic boundary behavior.For imaginary-time propagation, the supplied criterion is an absolute value or derivative below about 10^-10; for real-time propagation it is below about 10^-7.
- Numerical parameters: Crank-Nicolson is unconditionally stable, but DX and DT still require correlated trial values for good convergence.The programs provide preset values, while smaller space and time steps generally improve accuracy.
- Multidimensional programs: Six additional programs cover Cartesian 2D, Cartesian 3D, and axially symmetric 3D traps with imaginary- and real-time propagation.Fortran 90/95 versions are also included for the multidimensional programs.
- Use cases: Imaginary-time programs are recommended for stationary ground states, whereas real-time programs target non-equilibrium dynamics.The stated dynamical applications include soliton dynamics, expansion, and collapse dynamics.
- Parameter changes: Increasing nonlinearity may require more mesh points because the wave function extends over a larger spatial region.The supplied guidance recommends checking output files to verify that the wave function is adequately contained within the computational domain.
- Precision and cost: Computational time can be reduced by increasing the space and time steps and reducing the number of mesh points, at the cost of precision.When DX is doubled, the mesh count should be halved, while DT should increase by roughly a factor of four or more based on experimentation.
6. Numerical Results
Imaginary-time calculations obtain stationary GP results across one-, two-, and three-dimensional trap symmetries, reporting chemical potentials, rms sizes, and central wave-function values. The profiles show expansion with increasing nonlinearity and directional deformation under anisotropy, while comparisons generally agree with or improve upon earlier calculations.
- Stationary results: Imaginary-time programs calculate chemical potential, rms size or radius, and central wave-function values for stationary states across the considered trap symmetries.The reported cases include 1D, circularly symmetric and anisotropic 2D, radially or spherically symmetric 3D, axially symmetric 3D, and fully anisotropic 3D traps.
- One-variable cases: The 1D and radially symmetric 3D calculations use ℵ=627.42 and 627.4, respectively, with correlated space and time steps selected for convergence.The calculations use DX=0.0025 and DT=0.00002 in the reported table results.
- One-variable cases: As nonlinearity increases, repulsion increases and the wave function extends over a larger spatial domain in the radially symmetric 3D and 1D cases.Figure 2 displays profiles calculated with imagtimesph.F and imagtime1d.F for the nonlinearities listed in Tables 2 and 3.
- Two-dimensional cases: The Cartesian 2D results at ℵ=12.5484 converge with DX=0.02, but are less accurate than finer-mesh results because the larger space step reduces CPU time.A finer mesh is reported to increase accuracy while requiring more CPU time.
- Two-dimensional cases: For anisotropic 2D traps, κ=2 compresses the wave function in the y direction, while Table 5 reports results for both anisotropic and circularly symmetric cases.The calculations use imagtime2d.F and imagtimecir.F and establish more accurate results than Ref..
7. Summary and Conclusion
The paper presents split-step Crank–Nicolson algorithms with real- and imaginary-time propagation for time-dependent GP equations across multiple trap symmetries. It supplies Fortran implementations, evaluates stationary properties against prior results, and studies non-stationary oscillations, while noting straightforward adaptation to other equations and potentials.
- The split-step method treats spatial-derivative parts with Crank–Nicolson and handles derivative and non-derivative parts in independent steps.The paper characterizes the resulting calculations as highly stable and accurate.
- Real-time programs produce non-stationary oscillation results for multiple trap symmetries, including the 1D, spherical, circular, and axially symmetric cases shown in Fig. 5.The figure plots rms size versus time for these four cases.
- Real-time propagation solves non-stationary problems, whereas imaginary-time propagation yields accurate stationary ground states with smaller CPU time.Real-time evolution uses the original complex equation; imaginary-time propagation converts the equation to a real form.
- The study covers six harmonic-trap symmetries with both propagation types, providing twelve Fortran 77 programs and six Fortran 90/95 programs for two- and three-variable cases.The programs calculate stationary properties such as chemical potential and rms size and study oscillations after halving the initial nonlinearity.
- The programs are valid for the standard cubic-nonlinearity GP equation in a harmonic potential but can be adapted to other nonlinearities and potentials.Potential changes are made in INITIALIZE, while nonlinear changes are made in NONLIN.