{'fraction_non_alphanumeric': 0.04110530132437693, 'fraction_numerical': 0.01914049082096103, 'mean_word_length': 4.480912038489698, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fortran programs for the time-dependent Gross-Pitaevskii equation in a fully anisotropic trap

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).


Introduction
After a successful experimental detection of Bose-Einstein condensates (BEC) of dilute trapped bosonic alkali-metal atoms 7 Li, 23 Na, and 87 Rb [1,2] at ultra-low temperatures, there have been intense theoretical activities in studying properties of the condensate using the time-dependent mean-field Gross-Pitaevskii (GP) equation under different trap symmetries. Among many possibilities, the following traps have been used in various studies: three-dimensional (3D) spherically-symmetric, axially-symmetric and anisotropic harmonic traps, two-dimensional (2D) circularly-symmetric and anisotropic harmonic traps, and one-dimensional (1D) harmonic trap. The inter-atomic interaction leads to a nonlinear term in the GP equation, which complicates its accurate numerical solution, specially for a large nonlinearity. The nonlinearity is large for a fixed harmonic trap when either the number of atoms in the condensate or the atomic scattering length is large and this is indeed so under many experimental conditions. Special care is needed for the solution of the time-dependent GP equation with large nonlinearity and there has been an extensive literature on this topic [3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47, The time-dependent GP equation is a partial differential equation in space and time variables involving first-order time and second-order space derivatives together with a harmonic and a nonlinear potential term, and has the structure of a nonlinear Schrödinger equation with a harmonic trap. One commonly used procedure for the solution of the time-dependent GP equation makes use of a discretization of this equation in space and time and subsequent integration and time propagation of the discretized equation. From a knowledge of the solution of this equation at a specific time, this procedure finds the solution after a small time step by solving the discretized equation. A commonly used discretization scheme for the GP equation is the semi-implicit Crank-Nicolson discretization scheme [49,50,51] which has certain advantages and will be used in this work.
In the simplest one-space-variable form of the GP equation, the solution algorithm is executed in two steps. In the first step, using a known initial solution, an intermediate solution after a small interval of time ∆ is found neglecting the harmonic and nonlinear potential terms. The effect of the potential terms is then included by a first-order time integration to obtain the final solution after time ∆. In case of two or three spatial variables, the space derivatives are dealt with in two or three steps and the effect of the potential terms are included next. As the time evolution is executed in different steps it is called a split-step real-time propagation method. This method is equally applicable to stationary ground and excited states as well as non-stationary states, although in this paper we do not consider stationary excited states. The virtues of the semi-implicit Crank-Nicolson scheme [49,50,51] are that it is unconditionally stable and preserves the normalization of the solution under real-time propagation. A simpler and efficient variant of the scheme called the split-step imaginary-time propagation method obtained by replacing the time variable by an imaginary time is also considered. (The GP equation involves complex variables. However, after replacing the time variable by an imaginary time the resultant partial differential equation is real, and hence the imaginarytime propagation method involves real variables only. This trick leads to an imaginary-time operator which results in exponential decay of all states relative to the ground state and can then be applied to any initial trial wave function to compute an approximation to the actual ground state rather accurately. We shall use imaginary-time propagation to compute the ground state in this paper.) The split-step imaginary-time propagation method involving real variables yields very precise result at low computational cost (CPU time) and is very appropriate for the solution of stationary problems involving the ground state. The split-step real-time propagation method uses complex quantities and yields less precise results for stationary problems; however, they are appropriate for the study of non-equilibrium dynamics in addition to stationary problems involving excited states also.
The purpose of the present paper is to develop a simple and efficient algorithm for the numerical solution of the GP equation using time propagation together with the semi-implicit Crank-Nicolson discretization scheme [49,50,51] specially useful to newcomers in this field interested in obtaining a numerical solution of the time-dependent GP equation. Easy-to-use Fortran 77 programs for different trap symmetries with adequate explanation are also included. In case of two and three space variables, Fortran 90/95 programs are more compact in nature and these programs are also included. We include programs using both real-and imaginary-time propagation. For stationary ground states the imaginary-time method has a much quicker convergence rate compared to the real-time method and should be used for the calculation of chemical potential, energy and root-mean-square (rms) sizes. We calculate the chemical potential and rms sizes of the condensate for stationary problems and compare these results with those previously obtained by other workers for different trap symmetries. These results can be easily calculated in a decent PC using the Fortran programs provided. In addition to the results for the stationary states, the real-time propagation routines can also be used to study the non-stationary transitions, as in collapse dynamics [52] and non-equilibrium oscillation [30].
This paper is organized as follows. In Sec. 2 we present the GP equations with different traps that we consider in this paper, e.g., the 3D spherically-symmetric, 2D circularly-symmetric and 1D harmonic traps involving one space variable, the anisotropic 2D and axially-symmetric 3D harmonic traps in two space variables and the fully anisotropic 3D harmonic trap in three space variables. In Sec. 3 we elaborate the numerical algorithm for solving the GP equation in one space variable (the 1D, circularly-symmetric 2D, and spherically-symmetric 3D cases) and for calculating the chemical potential, energy and rms sizes employing both the real-and imaginary-time propagation methods. In Sec. 4 we present the same for solving the GP equation in two and three space variables (the anisotropic 2D and axially-symmetric and anisotropic 3D cases). In Sec. 5 we present a description of the Fortran programs, an explanation about how to use them, and some sample outputs. In Sec. 6 we present the numerical results for chemical potential, rms size, value of the wave function at the center, for the ground-state problem using the imaginary-time propagation routines and compare our finding with previous results for different trap symmetries in 1D, 2D, and 3D. We also present a study of non-stationary oscillation in some of these cases using the real-time propagation routines when the nonlinear coefficient in the GP equation with a stationary solution was suddenly reduced to half its value. Finally, in Sec. 7 we present a brief summary of our study.

Nonlinear Gross-Pitaevskii Equation
At zero temperature, the time-dependent Bose-Einstein condensate wave function Ψ ≡ Ψ(r; τ ) at position r and time τ may be described by the following mean-field nonlinear GP equation [1] i Here m is the mass of an atom and N the number of atoms in the condensate, g = 4π 2 a/m the strength of inter-atomic interaction, with a the atomic scattering length. The normalization condition of the wave function is dr|Ψ(r; τ )| 2 = 1.

Spherically-symmetric GP equation in 3D
In this case the trap potential is given by V (r) = 1 2 mω 2r2 , where ω is the angular frequency andr the radial distance. After a partial-wave projection the radial part ψ of the wave function Ψ can be written as Ψ(r; τ ) = ψ(r, τ ). After a transformation of variables to dimensionless quantities defined by r = √ 2r/l, t = τ ω, l ≡ ( /mω) and φ(r; t) ≡ ϕ(r; t)/r = ψ(r, τ )[l 3 /(2 √ 2)] 1/2 , the GP equation (1) in this case becomes where ℵ = 8 √ 2πN a/l. The purpose of changing the wave function from ψ to ϕ = rψ is a matter of taste and it has certain advantages. First, this transformation removes the first derivative ∂/∂r from the differential equation (2) and thus results in a simpler equation [34]. Secondly, at the origin r = 0, ψ is a constant, or ∂ψ/∂r = 0. But the new variable satisfies ϕ(0, t) = 0. Hence, while solving the differential equation (2), we can implement the simple boundary condition that as r → 0 or ∞, ϕ vanishes. The boundary condition for the differential equation in ψ will be a mixed one, e.g., the function ψ should vanish at infinity and its first space derivative should vanish at the origin. The normalization condition for the wave function is However, Eq. (2) is not the unique form of dimensionless GP equation in this case. Other forms of dimensionless equations have been obtained and used by different workers. For example, using the transformations r =r/l, t = τ ω, l ≡ ( /mω) and φ(r; t) ≡ ϕ(r; t)/r = ψ(r, τ )l 3/2 , the GP equation (1) becomes where ℵ = 4πN a/l with normalization (3). Finally, using the transformations r =r/l, t = τ ω/2, l ≡ ( /mω) and φ(r; t) ≡ ϕ(r; t)/r = ψ(r, τ )l 3/2 , the GP equation (1) becomes where ℵ = 8πN a/l with normalization (3). These three sets of dimensionless GP equations have been widely used in the literature and will be considered here. Equations (2), (4), and (5) allow stationary solutions ϕ(r; t) ≡ ϕ(r) exp(−iµt) where µ is the chemical potential. The boundary conditions for the solution of these equations are ϕ(0, t) = 0 and lim r→∞ ϕ(r, t) = 0 [49].

Circularly-symmetric GP equation in 2D
In the special case of circular symmetry the equations of Sec. 2.5 can be written in one-dimensional form. In this case κ = 1, and we introduce the radial variable r ≡ (x, y), and rewrite the wave function as ϕ(r). Then the GP equation (17) become The normalization of the wave function is 2π ∞ 0 drr|ϕ(r; t)| 2 = 1. In the circularly-symmetric case Eq. (19) becomes Finally, Eq. (20) can be written as The convenient boundary condition in this case is lim r→∞ ϕ(r; t) = 0 and dϕ(r; t)/dr| r=0 = 0 [12].
In this section we have exhibited GP equations for different trap symmetries. In the next section we illustrate the Crank-Nicolson method for the GP equation in one space variable, which is then extended to other types of equations in Sec. 4.

The GP Equation in the 1D and radially-symmetric 3D cases
To introduce the Crank-Nicolson Method [49,50,51] for GP equation we consider first the one-dimensional case of Sec. 2.4. The nonlinear GP equation (13) in this case can be expressed in the following form: where the Hamiltonian H contains the different linear and nonlinear terms including the spatial derivative.
(The spherically-symmetric GP equation in 3D has a similar structure and can be treated similarly.) We solve this equation by time iteration [49,50,51]. A given trial input solution is propagated in time over small time steps until a stable final solution is reached. The GP equation is discretized in space and time using the finite difference scheme. This procedure results in a set of algebraic equations which can be solved by time iteration using an input solution consistent with the known boundary condition. In the present split-step method [50] this iteration is conveniently done in few steps by breaking up the full Hamiltonian into different derivative and non-derivative parts.

Real-time propagation
The time iteration is performed by splitting H into two parts: Essentially we split Eq. (24) into We first solve Eq. (27) with a initial value ϕ(x; t 0 ) at t = t 0 to obtain an intermediate solution at t = t 0 + ∆, where ∆ is the time step. Then this intermediate solution is used as initial value to solve Eq. (28) yielding the final solution at t = t 0 + ∆ as ϕ(x; t 0 + ∆). This procedure is repeated n times to get the final solution at a given time t final = t 0 + n∆.
The time variable is discretized as t n = n∆ where ∆ is the time step. The solution is advanced first over the time step ∆ at time t n by solving the GP equation (24) with H = H 1 to produce an intermediate solution ϕ n+1/2 from ϕ n , where ϕ n is the discretized wave function at time t n . As there is no derivative in H 1 , this propagation is performed essentially exactly for small ∆ through the operation where nd (H 1 ) denotes time-evolution operation with H 1 and the suffix 'nd' denotes non-derivative. Next we perform the time propagation corresponding to the operator H 2 numerically by the semi-implicit Crank-Nicolson scheme (described below) [49]: The formal solution to (30) is which combined with Eq. (29) yields where CN denotes time-evolution operation with H 2 and the suffix 'CN' refers to the Crank-Nicolson algorithm. Operation CN is used to propagate the intermediate solution ϕ n+1/2 by time step ∆ to generate the solution ϕ n+1 at the next time step t n+1 = (n + 1)∆. The advantage of the above split-step method with small time step ∆ is due to the following three factors [50,51]. First, all iterations conserve normalization of the wave function. Second, the error involved in splitting the Hamiltonian is proportional to ∆ 2 and can be neglected and the method preserves the symplectic structure of the Hamiltonian formulation. Finally, as a major part of the Hamiltonian including the nonlinear term is treated fairly accurately without mixing with the delicate Crank-Nicolson propagation, the method can deal with an arbitrarily large nonlinear term and lead to stable and accurate converged result. Now we describe explicitly the semi-implicit Crank-Nicolson algorithm. The GP equation is mapped onto N x one-dimensional spatial grid points in x. Equation (24) is discretized with H = H 2 of (26) by the following Crank-Nicolson scheme [49,50,51]: where for the spherically-symmetric Eq.
has the advantage of taking an equal number of space points on both sides of x = 0 in the 1D case setting the point N x /2 at x = 0.) Equation (33) is the explicit form of the formal Eq. (30). This scheme is constructed by approximating ∂/∂t by a two-point formula connecting the present (n + 1/2) to future (n + 1). The spatial partial derivative ∂ 2 /∂x 2 is approximated by a three-point formula averaged over the present and the future time grid points. This procedure results in a series of tridiagonal sets of equations (33) in ϕ n+1 i+1 , ϕ n+1 i , and ϕ n+1 i−1 at time t n+1 , which are solved using the proper boundary conditions. The Crank-Nicolson scheme (33) possesses certain properties worth mentioning [50,51]. The error in this scheme is both second order in space and time steps so that for small ∆ and h the error is negligible. This scheme is also unconditionally stable [51]. The boundary condition at infinity is preserved for small values of ∆/h 2 [51].
The tridiagonal equations emerging from Eq. (33) are written explicitly as [49] A and All quantities in b i refer to time step t n+1/2 and are considered known. The only unknowns in Eq. (34) are the wave forms ϕ n+1 i±1 and ϕ n+1 i at time step t n+1 . To solve Eq. (34), we assume the one-term forward recursion relation where α i and β i are coefficients to be determined. Substituting Eq. (36) in Eq. (34) we obtain which leads to the solution with From Eqs. (36) and (38) we obtain the following backward recursion relations for the coefficients α i and β i We shall use the recursion relations (38), (39) and (40) in a backward sweep of the lattice to determine α i and β i for i running from N x − 2 down to 0. The initial values chosen are α Nx−1 = 0, β Nx−1 = ϕ n+1 Nx . This ensures the correct value of ϕ at the last lattice point. After determining the coefficients α i , β i and γ i , we can use the recursion relation (36) from i = 0 to N x − 1 to determine the solution for the entire space range using the starting value ϕ n+1 0 (=0) known from the boundary conditions. The value at the last lattice point is also taken to be known (= 0). Thus we have determined the solution by using two sets of recursion relations across the lattice each involving about N x operations.
In the numerical implementation of CN real-time propagation the initial state at t = 0 is usually chosen to be the analytically known solution of the harmonic potential with zero nonlinearity: ℵ = 0. In the course of time iteration the nonlinearity is slowly introduced until the desired final nonlinearity is attained. This procedure will lead us to the final solution of the problem.

Imaginary-time propagation
Although, real-time propagation as described above has many advantages, in this approach one has to deal with complex variables for a complex wave function for non-stationary states. For stationary ground state the wave function is essentially real-and the imaginary-time propagation method dealing with real variables seems to be convenient. In this approach time t is replaced by an imaginary quantity t = −it and Eq. (24) now becomes In this equationt is just a mathematical parameter.
. Hence, if we start with as arbitrary initial ϕ(x;t) which can be taken as a linear combination of all eigenfunctions of H, then upon imaginary-time propagation all the eigenfunctions will decay exponentially with time. However, all the excited states with larger E i will decay exponentially faster compared to the ground state with the smallest eigenvalue. Consequently, after some time only the ground state survives. As all the states are decaying with time during imaginary-time propagation, we need to multiply the wave function by a number greater than unity to preserve its normalization so that the solution does not go to zero.
The imaginary-time iteration is performed by splitting H into two parts as before: H = H 1 + H 2 , with H 1 and H 2 given by Eqs. (25) and (26). It is realized that the entire analysis of Sec. 3.1.1 remains valid provided we replace i by 1 in Eq. (29) and by −1 in the remaining equations. However, there appears one trouble. The CN real-time propagation preserves the normalization of the wave function, whereas the CN imaginarytime propagation does not preserve the normalization. This problem can be circumvented by restoring the normalization of the wave function after each operation of Crank-Nicolson propagation. Once this is done the imaginary-time propagation method for stationary ground state problems yields very accurate result at low computational cost.
Compared to the real-time propagation method, the imaginary-time propagation method is very robust. The initial solution in the imaginary-time method could be any reasonable solution and not the analytically known solution of a related problem as in the real-time method. Also, the full nonlinearity can be added in a small number of time steps or even in a single step and not in a large number of steps as in the real-time method. In the programs using the imaginary-time method we include the nonlinearity in a single step. These two added features together with the use of real algorithm make the imaginary-time propagation method very accurate with quick convergence for stationary ground states as we shall see below.

The GP Equation in the circularly-symmetric 2D case
The Crank-Nicolson discretization for real-and imaginary-time propagation for the circularly-symmetric GP equation (21) is performed in a similar fashion as for Eq. (24), apart from the difference that here we also have a first derivative in space variable in addition to the second derivative. Another difference is that the wave function is not zero at r = 0. In the case of Eq. (24) we took the boundary condition as ϕ(x; t) = 0 at the boundaries. For the circularly-symmetric case, the convenient boundary conditions are lim r→∞ ϕ(r; t) = 0 and dϕ(r; t)/dr| r=0 = 0.
We describe below the Crank-Nicolson discretization and the solution algorithm in this case for the following equation required for the solution of Eq. (21). The remaining procedure is similar to that described in detail above in the one-dimensional case. Equation (42) is discretized by the following Crank-Nicolson scheme as in Eq. (33) [49,50,51]: where again ϕ n i = ϕ(r i ; t n ), r ≡ r i = ih, i = 0, 1, 2, ..., N r and h is the space step. This scheme is constructed by approximating ∂/∂x by a two-point formula averaged over present and future time grid points. The discretization of the first-order time and second-order space derivatives is done as in Eq. (33). This procedure results in the tridiagonal sets of equations (43) , and ϕ n+1 i−1 at time t n+1 , which are solved using the proper boundary conditions.
The tridiagonal equations emerging from Eq. (43) are written explicitly as where and All quantities in b i refer to time step t n+1/2 and are considered known. The only unknowns in Eq. (44) are the wave forms ϕ n+1 i±1 and ϕ n+1 i at time step t n+1 . To solve Eq. (44), we assume the one-term backward recursion relation where α i and β i are coefficients to be determined. Substituting Eq. (46) in Eq. (44) we obtain which leads to the solution with From Eqs. (46) and (48) we obtain the following forward recursion relations for the coefficients α i and β i We shall use the recursion relations (48), (49) and (50) in a forward sweep of the lattice to determine α i and β i for i running from 1 to N r − 1. The initial values chosen are α 1 = 1, β 1 = 0. This ensures the correct value of the space derivative of ϕ(r; t) = 0 at r = 0. After determining the coefficients α i , β i and γ i , we can use the recursion relation (46) from i = N r to 1 to determine the solution for the entire space range using the starting value ϕ n+1 (N r ) (=0) known from the boundary condition. Thus we have determined the solution by using two sets of recursion relations across the lattice each involving about N r operations.

Chemical potential
For stationary states the wave functions for the 1D case have the trivial time dependence ϕ(x; t) ≡ ϕ(x) exp(−iµt), where µ is the chemical potential. Substituting this condition in Eq. (13) we obtain Assuming that the wave form is normalized to unity ∞ −∞φ 2 (x)dx = 1, the chemical potential can be calculated from the following expression obtained by multiplying Eq. (51) byφ(x) and integrating over all space where the second derivative has been simplified by an integration by parts. All the programs also calculate the many-body energy, which is of interest. The analytical expression for energy is the same as that of the chemical potential but with the nonlinear term multiplied by 1/2, e.g., [ The programs will write the value of energy as output. However, we shall not study or tabulate the results for energy and we shall not write the explicit algebraic expression of energy in the case of other trap symmetries. The GP equation (2) with spherically-symmetric potential is also an one-variable equation quite similar in structure to the one-dimensional equation (13) considered above. Hence the entire analysis of Secs. 3.1.1, 3.1.2, and 3.3 will be applicable in this case with ϕ(x; t) replaced by ϕ(r; t)/r in the nonlinear term. Now if we consider stationary states of the form ϕ(r, t) ≡φ(r) exp(−iµt/ ), the expression for the chemical potential becomes We shall use Eqs. (52) and (54) for the calculation of chemical potential from Eqs. (13) and (2). The energy will be calculated from Eq. (53). The expressions for chemical potential in other cases can be written down in a straight-forward fashion.

Split-
Step Crank-Nicolson method in two and three space variables

Anisotropic GP equation in 2D
In this case the GP equation (17) can be written as Now we adopt a policy quite similar to that elaborated in Sec. 3.1.1 where the Hamiltonian was broken into two parts and where the time propagation over time step ∆ using the two parts were carried out alternatively.
The same procedure will be adopted in the present case where we perform the time propagation using the pieces H 1 , H 2 , and H 3 of the Hamiltonian successively in three independent time sub-steps ∆ to complete a single time evolution over time step ∆ of the entire GP Hamiltonian H. The time propagation over H 1 is performed as in Eq. (29) and those over H 2 and H 3 as in Eqs. (31) and (33). The chemical potential in this case for the stationary state ϕ(x, y; t) ≡φ(x, y) exp(−iµt) can be written as where we have performed integrations by parts to obtain this form for the chemical potential in terms of first derivatives only.

Axially-symmetric GP equation in 3D
In this case the GP equation (10) can be written as The Hamiltonian H can be conveniently broken into three pieces Now we adopt a policy quite similar to that elaborated in Sec. 4.1 where the Hamiltonian was broken into three parts and where the time propagation over time step ∆ using the three parts were carried out alternatively. The time propagation over H 1 is performed as in Eq. (29) and those over H 2 and H 3 as in Eqs. (43) and (33). The chemical potential in this case for the stationary state ϕ(ρ, z; t) ≡φ(ρ, z) exp(−iµt) can be written as where we have again used integrations by parts to simplify the final expression.

Anisotropic GP equation in 3D
In this case the GP equation (6) can be written as Now we adopt a policy quite similar to that elaborated in Sec. 3.1.1 where the Hamiltonian was broken into two parts and where the time propagation over time step ∆ using the two parts were carried out alternatively.
The same procedure will be adopted in the present case where we perform the time propagation using the pieces H 1 , H 2 , H 3 and H 4 of the Hamiltonian successively in four independent time sub-steps ∆ to complete a single time evolution over time step ∆ of the entire GP Hamiltonian H. The time propagation over H 1 is performed as in Eq. (29) and those over H 2 , H 3 and H 4 as in Eqs. (31) and (33). The chemical potential in this case for the stationary state ϕ(x, y, z; t) ≡φ(x, y, z) exp(−iµt) can be written as where we have again used integrations by parts to simplify the final expression. Now we describe the function of the subroutines, which the user should not need to change. The subroutine NORM calculates by Simpson's rule the normalization of the wave function and sets the normalization to unity, the preassigned value. The real-time propagation preserves the normalization of the wave function and hence the subroutine NORM is not used during time propagation. The subroutines RAD and CHEM calculate the rms size (length or radius), and chemical potential (and energy) of the system. The subroutines COEF and LU together implement the time propagation with the spatial and temporal time derivative terms. The subroutine NU performs the time propagation with the nonlinear term and the potential. In the imaginary-time program the action of the subroutine LU does not preserve normalization and hence each time the subroutine LU is called, the subroutine NORM has to be called to set the normalization of the wave function back to unity. (This is not necessary in the real-time programs which preserve the norm.) The subroutine NONLIN calculates the nonlinear term. The subroutine DIFF calculates the space derivative of the wave function by Richardson's extrapolation formula needed for the computation of the chemical potential and energy. The function SIMP does the numerical space integrations with the Simpson's rule.

Description of Numerical Programs
The programs implement the splitting method described in Sec. 3.1.1 and 3.2 and calculate the wave function, chemical potential, size, normalization, etc. The number of points in the one-dimensional space grid represented by the integer variable N has to be chosen consistent with space step DX such that the total space covered N×DX is significantly larger than the size of the condensate so that at the boundaries the wave function attains the asymptotic limits (e.g., the absolute value of the wave function or its space derivative becomes less than 10 −10 or so for imaginary-time propagation and less than 10 −7 for real-time propagation. Note that in the 1D and spherically-symmetric 3D problems we are using the asymptotic condition that the wave function is zero at both boundaries. In case of the circularly-symmetric 2D problem we use a mixed boundary condition, e.g., at origin the derivative of the wave function is zero and at infinity the wave function is zero.) In the imaginary-time routine the total space covered should be about 1.5 times the extension of the condensate. In the real-time routines, to obtain good precision the total space covered should be at least 2 times the extension of the condensate. This is because the imaginary-time routine is more precise and the solution attains its limiting asymptotic value at the boundary very rapidly as the total space covered is increased. The real-time routine is less accurate and one has to go to a larger distance before the solution or its derivative drops to zero. A couple of runs with a sufficiently large N and large DX are recommended to have an idea of the size/extension of the system. A smaller value of DX leads to a more accurate result provided an appropriate DT is chosen.
The Crank-Nicolson method, described, for example, by Eq. (33), is unconditionally stable for all ∆/h 2 . Nevertheless, for a numerical application to a specific problem one has to fix the time step DT (= ∆) and space step DX (= h) for good convergence.  Table 1.) The total number of space points N in calculation has to be fixed in the line, e.g., "PARAMETER (N = 6000, N2 = N/2, NX = N-1)", in the MAIN program, so that the final wave function is within the space range and its value is negligibly small at the boundaries. The total numbers of time iterations are fixed in the line, e.g., "PARAMETER (NSTP = 500000, NPAS = 10000, NRUN = 20000)", in the MAIN program as described below in Sec. 5.3.
The integer parameter OPTION should be set 1, 2 or 3 in the MAIN routine for solution of equations of type (5), (4), or (2), [or for solution of equations of type (16), (15), or (13),] respectively. The difference between these three types of equations is in the values of the coefficients in the first two terms. The nonlinear term calculated in the subroutine NONLIN is of the form |ϕ(r, t)| 2 with coefficient G. If a different type of nonlinear term (with different functional dependence on the wave function) is to be introduced, it should be done in the subroutine NONLIN. Otherwise, this subroutine should not be changed. A different type of nonlinear term is appropriate for a Fermi super-fluid [53,54] or a Tonks-Girardeau gas [55].

GP equation in two and three space variables
In addition to the programs in one space variable we have six programs in two and three space variables. The programs imagtime2d.F and realtime2d.F apply to 2D Cartesian space using imaginaryand real-time propagations, respectively. Similarly, imagtime3d.F and realtime3d.F apply to 3D Cartesian space using imaginary-and real-time propagation, respectively. Finally, programs imagtimeaxial.F and realtimeaxial.F apply to an axially symmetric trap in 3D. The Fortran 90/95 versions of these programs are somewhat more condensed and provide some advantage. Hence we also include these programs as imagtime2d.f90, realtime2d.f90, imagtime3d.f90, realtime3d.f90, imagtimeaxial.f90, and realtimeaxial.f90. The output of the Fortran 90/95 programs are identical with their corresponding Fortran 77 versions. All these programs are written using a very similar logic used in the programs in one space variable. So most of the considerations described earlier also applies to these cases. We describe the principal differences below.
In case of total number of space points N, one now has the variables NX, NY, and NZ for total number of space points in X, Y and Z directions in case of three space variables. In case of two space variables the Z component is absent. These variables should be chosen equal to each other for a nearly symmetric case. However, if the problem is anisotropic in space, these variables can and should be chosen differently. Similarly, the space variable X(I) is now replaced by space variables X(I), Y(J), and Z(K) in three directions. Now there are space steps DX, DY, and DZ in place of DX in one space variable. The variables DX, DY, and DZ can now be chosen differently in three directions in case of an anisotropic problem. However, they should be chosen together with NX, NY, and NZ so that the wave function lies entirely inside the chosen space and becomes negligibly small at the boundaries. If the spatial extension of the wave function is much smaller in one space direction than another, one should take a smaller space step in the direction in which the spatial extension of the wave function is small. The potential V and wave function CP are now functions of 2 or 3 space variables and are represented by matrices in place of a column in the case of one space variable. The variables AL (or KAP) and BL (or LAM) now define the anisotropy of the trap and are used in the subroutine INITIALIZE to define the trap. The subroutine LU is now replaced by LUX, LUY, and LUZ to implement the solution in different space directions. In the imaginary-time propagation each time the subroutines LUX, LUY, and LUZ operate, one has to set the normalization of the wave function to unity by calling the subroutine NORM.

Instruction to use the programs
The programs, as supplied, solve the GP equations for a specific value of nonlinearity and write the wave function, chemical potential, energy, rms size or radius, wave function at the center, and nonlinearity, for specific values of space and time steps. The real-and imaginary-time programs for a specific equation, for example, spherically-symmetric 3D equation, employ similar set of parameters like space and time steps. The real-time programs use a larger value of N, so that the discretization covers a larger region in space. In all cases the supplied programs solve an equation of type (4) with a factor of 1/2 in front of the space derivative, selected by setting the integer parameter OPTION = 2 in the MAIN routine. Other types of equations can be obtained by setting OPTION = 1 for equations of type (5) or = 3 for equations of type (2).
For solving a stationary ground state problem, the imaginary-time programs are far more accurate and should be used. The real-time programs should be used for studying non-equilibrium problems often using an initial wave function calculated by the imaginary-time program. The non-equilibrium problems include the study of soliton dynamics [56], expansion [57], collapse dynamics [52], and other types of problems.
Each program is preset at a fixed nonlinearity G0 (= G), correlated DX-DT values and NSTP, NPAS, and NRUN. Smaller the steps DX and DT, more accurate will be the result. The correlated values of DX and DT on a data line should be found by trial to obtain good convergence. Each supplied program produces result up-to a desired precision consistent with the parameters employed − G0, DX, DT, N, NSTP, NPAS, and NRUN. If the nonlinearity G0 is increased, one might need to increase N to achieve similar precision. In many cases one may need an approximate solution (with lower accuracy) involving less CPU time or one may need to solve the GP equation for a different value of nonlinearity.
(a) If G0 is reduced, just change the card defining G0. However, if G0 is increased, changing the value of G0 may not be enough. For an increased G0, the wave function extends to a larger region in space. One may need to increase the "Number of space mesh points". For a new G0, just plot the output file fort.3 for all programs except realtime3d.F, realtime3d.f90, imagtime3d.F and imagtime3d.f90, where one should plot output files fort.11, fort.12, and fort.13 and see that the wave function is fully and adequately accommodated in the space domain. If not, one needs to increase the number in input cards "Number of space mesh points" until the wave function is fully and adequately accommodated in the space domain.
(b)The CPU time involved can be reduced by sacrificing the precision. This can be done by increasing the space step(s) and time step and reducing the "Number of space mesh points". If the space step DX is increased by a factor of f = 2, the number of space mesh points should be reduced by the same factor. The time step DT should be increased by a larger factor, more like f**2. The optimum increase in time The integer parameter NSTP refers to the number of time iterations during which the nonlinear term is slowly introduced in the real-time propagation. This number should be large (typically more than 100,000 for small nonlinearity, for larger nonlinearity could be 1000,000 ) for good convergence; this means that the nonlinearity should be introduced in small amounts over a large number of time iterations. In real-time propagation, NPAS refers to certain number of time iterations with the constant nonlinear term already introduced in NSTP and should be small (typically 1000). NRUN refers to time iterations with a modified nonlinearity so as to generate a non-equilibrium dynamics. In the imaginary-time propagation the parameters NPAS and NRUN refer to certain number of time iterations with the constant nonlinear term already introduced in one step and should be large (typically NPAS = 200,000 or more) for good convergence.

Output files
The programs write via statements WRITE(1,*), WRITE(2,*) WRITE(3,*) in Files 1, 2, and 3, respectively, the initial stationary wave function and that after NSTP, and NPAS time iterations for real-time propagation and after NPAS and NRUN time iterations for imaginary-time propagation. File 3 gives the final stationary wave function of the calculation. However, in the case of the anisotropic 3D programs real-time3d.F and imagtime3d.F, sections of the wave function as plotted in Fig. 4 (b) are written in Files 1, 2, and 3 before NSTP (realtime3d.F) and after NPAS (imagtime3d.F) iterations, and in Files 11, 12, and 13 after NPAS (realtime3d.F) and NRUN (imagtime3d.F) iterations.
In the real-time program a non-stationary oscillation is initiated by suddenly modifying the nonlinearity from G to G/2 after NPAS time iterations. During NRUN time iterations the non-stationary dynamics is studied. The real-time programs write on File 8 the running time and rms size during non-stationary oscillation.
In addition, these programs write in File 7 the values of nonlinearity G, space steps DX, DY, DZ, time step DT, number of space mesh points N, number of time iterations NSTP, NPAS, and NRUN together with the values of normalization of the wave function, chemical potential, energy, rms size, value of the wave function at center, nonlinearity coefficient G.
Below we provide some sample output listed on File 7 from the different programs using OPTION = 2.  Table 1 Convergence of result in the 1D (X) and radially-symmetric 3D (r) cases for nonlinearities ℵ = 627. 42  In this subsection we present results for the stationary ground state problem calculated with the imaginarytime programs. All numerical results presented in this paper are for the nonlinear equations with a factor Table 2 The chemical potential µ, rms radius rrms, and the wave function φ(0) at the center for various nonlinearities in the radially symmetric 3D case calculated using the program imagtimesph.F [Eq. (4), OPTION 2]. Table completed with  In each case a sufficiently large number of space points N is to be taken, so that the space domain of integration covers the extension of the wave function adequately. We exhibit in Table 1, the results for the wave function at center, rms size, and chemical potential in these two cases for a fixed nonlinearity for different DX and DT. We find that convergence is achieved up to six significant digits after the decimal point with space step DX = 0.0025 and time step DT = 0.00002. In Fig. 1 we plot the relative percentage error in chemical potential µ for various space steps. The percentage error rapidly reduces as space step is reduced.

Stationary Problem
In Table 2 we exhibit the chemical potential, rms radius, and the wave function at the center for various nonlinearities in the spherically symmetric 3D case using space step DX = 0.0025 and time step DT = 0.00002 in the imaginary-time propagation program. From the table we find that the results are in agreement with those calculated in Refs. [3,4]. Table 3 The chemical potential µ, rms size xrms, and the wave function ϕ(0) at the center for various nonlinearities in the 1D case calculated using the program imagtime1d.F [Eq. (15), OPTION 2]. In Table 3 we exhibit the chemical potential, rms size, and the wave function at the center for various nonlinearities in the 1D case using space step DX = 0.0025 and time step DT = 0.00002. From the table we find that the results are in agreement with those calculated in Ref. [4].    In Figs. 2 (a) and (b) we plot the wave function profiles for the radially-symmetric 3D and 1D cases calculated using the programs imagtimesph.F and imagtime1d.F, respectively, for different nonlinearities presented in Tables 2 and 3. As the nonlinearity is increased the system becomes more repulsive and the wave function extends to a larger domain in space.
In Table 4 we present the results for the wave function at center, rms radius, and chemical potential for the Cartesian 2D case with nonlinearity ℵ = 12.5484 using the imaginary-time propagation program imagtime2d.F. We find that desired convergence is achieved with space step DX = 0.02. (Note that the converged result in this case is less accurate than those in Tables 1, 2 and 3 as we have used a larger space step in order to keep the CPU time small. A finer mesh will increase the accuracy requiring a larger CPU time.) In Table 5 we present results for the wave function at the center, rms size, and chemical potential for different nonlinearities in the anisotropic and circularly-symmetric 2D case calculated using the imaginarytime routine imagtime2d.F and imagtimecir.F, respectively. In the anisotropic case we used space step 0.02 and time step 0.0001 and in the circularly-symmetric case we used space step 0.0025 and time step 0.00002. We also compare these results with those of Ref. [4] and establish more accurate results in the present numerical calculation.
In Figs. 3 (a) and (b) we plot the wave function profiles for the circularly-symmetric and anisotropic 2D cases using the programs imagtimecir.F and imagtime2D.F, respectively. In the anisotropic 2D case the anisotropy κ = 2 and nonlinearity ℵ = 12.4584. Because of anisotropy the wave function in Fig. 3 (b) is compressed in the y direction (note the different scales in x and y directions of the plot.) Table 5 The chemical potential µ, rms size rrms, and the wave function ϕ(0) at the center for various nonlinearities in the anisotropic 2D case obtained using the programs imagtime2d.F [Eq. (19), OPTION 2] and imagtimecir.F [Eq. (22), OPTION 2]. The case κ = 1 represents circular symmetry and κ = 1 corresponds to anisotropy.   Table 5.
In Table 6 we present results for the wave function at the center, rms sizes, and chemical potential for different nonlinearities in the axially-symmetric 3D case with κ = 1 and λ = 4 calculated using the imaginary-time routine imagtimeaxial.F and space step Dρ = DZ = 0.02 and time step DT = 0.0004. We compare with the results of Ref. [4] and establish more accurate results.
In Figs. 4 (a) and (b) we plot the wave function profiles for the axially-symmetric and fully anisotropic 3D cases using the programs imagtimeaxial.F and imagtime3d.F, respectively. In the axially-symmetric 3D case the anisotropy κ = 1, λ = 4 and nonlinearity ℵ = 1881 were employed. In the fully anisotropic 3D case the anisotropies ν = 1, κ = √ 2, and λ = 2 and nonlinearities ℵ = 22.454, 359.26 and 11496.3 are used. The effect of the anisotropy is explicit in the 3D case in generating different profiles of the wave function along x, y, and z directions for a fixed nonlinearity.
In Table 7 we exhibit the results for our calculations with the fully anisotropic potential together with those for the spherically-symmetric potential in three dimensions. The results for the spherically-symmetric potential are also calculated using the one-dimensional radially symmetric imaginary-time program imagtimesph.F in addition to the fully anisotropic program imagtime3d.F. In the anisotropic case the present results are consistent with those of Ref. [5]. In the spherical case the two sets of the present results (calculated with the spherically-symmetric and anisotropic programs) as well as those of Ref. [5] are consistent with each other.
The input to our calculation is the nonlinearity coefficient ℵ, which is related to the scattering length a, number of atoms N and harmonic oscillator length l via ℵ = 4πaN/l in Eq. (8). We provide the nonlinearity values of our calculations. Although the present results are in agreement with those of Ref. [5], a very precise comparison of the two calculations is not to the point as Schneider and Feder did not provide the nonlinearity coefficient ℵ used in their calculation. Table 7 The chemical potential µ, rms sizes, and wave function ϕ(0) at the center for various number N of condensate of Na atoms. The constants used are m(Na) = 38.175 × 10 −27 kg, a(Na) = 2.75 nm. In all cases the input to numerical calculation was the nonlinearity coefficient ℵ = 4πN a/l shown below. For the spherically-symmetric case, we solved the radially symmetric program imagtimesph.F [Eq. (4), OPTION 2], (in addition to the 3D anisotropic program imagtime3d.F setting equal frequencies in all three directions with DX = DY =DZ = 0.05 and DT = 0.0004) using [5,59]

Non-stationary Oscillation
The real-time propagation programs calculate the stationary states under different trap symmetries. However, they are less efficient than the imaginary-time propagation programs in this task requiring more CPU time and producing less accurate results. However, unlike the imaginary-time propagation programs, the real-time programs can produce time evolution of non-stationary states also and next we present results of such time evolution using the real-time propagation programs under different trap symmetries.
In this subsection we present results for non-stationary oscillation obtained with the use of the real-time programs. After the calculation of the stationary profile, the nonlinearity is suddenly reduced to half. The wave function is no longer an eigenstate of the new nonlinear equation. This sets the system into nonstationary oscillation which continues for ever. In Fig. 5 we plot the rms size of the wave function vs. time t showing this oscillation using the output from File 8 for (a) 1D case (using program realtime1d.F), (b) radially-symmetric 3D case (using program realtimesph.F), (c) Cartesian 2D case with anisotropy κ = 1 (using program realtime2d.F), and (d) 3D axially-symmetric case (using program realtimeaxial.F) with respective nonlinearities ℵ = 62.742, 125.484, 12.5484, and 18.81. The rms size at t = 0 is the rms size of the stationary wave function obtained after NPAS time iterations.
Because of transverse instability, the real-time program realtime3d.F in 3D does not lead to stable sinusoidal oscillation as in other cases for a large change in nonlinearity (nonlinearity reduced to half of its initial value) as shown in Fig. 5. Only for small perturbation a sinusoidal oscillation is observed. However, we do not present a systematic study of such oscillation.

Summary and Conclusion
In this paper we describe a split-step method for the numerical solution of the time-dependent nonlinear GP equation under the action of a general anisotropic 3D trap using real-and imaginary-time propagation. Similar methods for 1D and anisotropic 2D traps are also described. The time propagation is carried with an initial input. The full Hamiltonian is split into several spatial derivative and a non-derivative parts. The  spatial derivative parts are treated by the Crank-Nicolson method. Different spatial derivative and nonderivative parts are dealt in independent steps. This, so called split-step, method leads to highly stable and accurate results. We considered two types of time iterations − real-time propagation and imaginary-time propagation. In the real-time propagation time evolution is performed with the original complex equation. The numerical algorithm in this case requires the use of complex variable but produces solution of non-stationary problems. In the imaginary-time propagation, the time variable is replaced by i (= √ −1) times a new time variable, consequently the GP equation becomes real. The numerical solution of this equation can no longer yield the solution of non-stationary problems; but yields very accurate solution of stationary ground state problems only, requiring much smaller CPU time.
We provide the numerical algorithm in detail in 1D, 2D, and 3D for real-and imaginary-time propagations. We consider six different harmonic oscillator trap symmetries, e.g., a 1D trap, a circularly-symmetric 2D trap, a radially-symmetric 3D trap, an anisotropic trap in 2D, an axially-symmetric 3D trap, and an anisotropic trap in 3D. Each of these cases are treated with real-and imaginary-time propagation algorithms resulting in twelve different Fortran 77 programs supplied. We use the imaginary-time propagation programs to provide results for different stationary properties of the condensate (chemical potential, rms size, etc) in 1D, 2D, 3D, for different nonlinearities ℵ and compare with previously obtained results [4,5]. In addition we study a non-stationary oscillation initiated by suddenly altering the nonlinearity to half its initial value on these preformed condensates using the real-time propagation programs. In addition six Fortran 90/95 programs are supplied in the case of two and three space variables.
Although the present programs are valid for the standard GP equation with cubic nonlinearity in a harmonic potential, they can be easily adopted for other types of bosonic [55] or fermionic equations [53,54] with different nonlinearities and under different types of potentials. To change the potential one should change the variable V in the subroutine INITIALIZE and the change in the nonlinearity can be performed in the subroutine NONLIN.