{'fraction_non_alphanumeric': 0.04393389762192664, 'fraction_numerical': 0.01724535037715207, 'mean_word_length': 4.38944918541505, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical Solution of the Gross-Pitaevskii Equation for Bose-Einstein Condensation

We study the numerical solution of the time-dependent Gross-Pitaevskii equation (GPE) describing a Bose-Einstein condensate (BEC) at zero or very low temperature. In preparation for the numerics we scale the 3d Gross-Pitaevskii equation and obtain a four-parameter model. Identifying `extreme parameter regimes', the model is accessible to analytical perturbation theory, which justifies formal procedures well known in the physical literature: reduction to 2d and 1d GPEs, approximation of ground state solutions of the GPE and geometrical optics approximations. Then we use a time-splitting spectral method to discretize the time-dependent GPE. Again, perturbation theory is used to understand the discretization scheme and to choose the spatial/temporal grid in dependence of the perturbation parameter. Extensive numerical examples in 1d, 2d and 3d for weak/strong interactions, defocusing/focusing nonlinearity, and zero/nonzero initial phase data are presented to demonstrate the power of the numerical method and to discuss the physics of Bose-Einstein condensation.


Introduction
Recent experimental advances in achieving and observing Bose-Einstein condensation (BEC) in trapped neutral atomic vapors [5,22,12] have spurred great excitement in the atomic physics community and renewed the interest in studying the collective dynamics of macroscopic ensembles of atoms occupying the same oneparticle quantum state [46,21,33]. The condensate typically consists of a few thousands to millions of atoms which are confined by a trap potential. In fact, beside the effects of the internal interactions between the atoms, the macroscopic behavior of BEC matter is highly sensitive to the shape of this external trapping potential. Theoretical predictions of the properties of a BEC like the density profile [11], collective excitations [25] and the formation of vortices [50] can now be compared with experimental data [5,37,45]. Needless to say that this dramatic progress on the experimental front has stimulated a wave of activity on both the theoretical and the numerical front.
The properties of a BEC at temperatures T much smaller than the critical condensation temperature T c [39] are usually well modeled by a nonlinear Schrödinger equation (NLSE) for the macroscopic wave function [39,33] known as the Gross-Pitaevskii equation (GPE) [34,48], which incorporates the trap potential as well as the interactions among the atoms. The effect of the interactions is described by a mean field which leads to a nonlinear term in the GPE. The cases of repulsive and attractive interactions -which can both be realized in the experiment -correspond to defocusing and focusing nonlinearities in the GPE, respectively. Note that equations very similar to the GPE also appear in nonlinear optics where an index of refraction, which depends on the light intensity, leads to a nonlinear term like the one encountered in the GPE.
There has been a series of recent studies which deals with the numerical solution of the time-independent GPE for the ground state and the time-dependent GPE for finding the dynamics of a BEC. Bao et al. [9] presented a general method to compute the ground state solution via directly minimizing the energy functional and used it to compute the ground state of the GPE in different cases. Edwards et al. introduced a Runge-Kutta type method and employed it to solve the spherically symmetric timeindependent GPE [26]. Adhikari [1,2] used this approach to obtain the ground state solution of the GPE in 2d with radial symmetry. Other approaches include a finite difference method proposed by Chiofalo et al. [18] and a simple analytical method proposed by Dodd [23]. For the numerical solution of the time-dependent GPE only few methods are available, a particle-inspired scheme proposed by Cerimele et al. in [16] and a finite difference method used by Ruprecht et al. [51].
In this paper, we take the 3d Gross-Pitaevskii equation, make it dimensionless to obtain a four-parameter model, use (singular) perturbation theory to discuss semiclassical asymptotics, to approximately reduce it to a 2d GPE and a 1d GPE in certain limits, and discuss the approximate ground state solution of the GPE in two extreme regimes: (very) weak interactions and strong repulsive interactions, again using perturbation methods. Numerical computations for similar physical set ups can to a large extent be found in the physical literatures cf. [11], however they are included here since perturbation theory gives a systematic way to obtain rigorously (validate) these approximations and since they are used as important preparatory steps for the numerical simulations (approximate ground states usually serve as initial data). Then we use the time-splitting spectral method, which was studied in Bao et al. [7,8] for the Schrödinger equation in the semiclassical regime, to discretize the time-dependent GPE. The merit of the numerical method is that it is explicit, unconditionally stable, time reversible, time-transverse invariant, and conserves the position density. In fact, the spectral method has shown great success in solving problems arising from many areas of science [32,13] due to its spatially spectral accuracy. The split-step procedure was presented for differential equations in [54] and applied to Schrödinger equations [35,55,28] and the KDV equation [56], as well as used with an iterative procedure for optical fibers [4]. In this paper we perform grid control using singular perturbation theory and present extensive numerical examples in 1d, 2d and 3d for weak/strong interactions, defocusing/focusing nonlinearity, and zero/nonzero initial phase data. The paper is organized as follows. In section 2 we start out with the 3d GPE, scale it to get a four-parameter model, show how to reduce it to lower dimensions and give the approximate ground state solution in the two mentioned extreme regimes of (very) weak interactions and strong repulsive interactions and discuss semiclassical asymptotics. In section 3 we present the time-splitting spectral method for the GPE. In section 4 numerical tests of the GPE for different cases including weak/strong interactions, defocusing/focusing nonlinearity, and zero/nonzero initial phase data are presented. In section 5 a summary is given.

Gross-Pitaevskii equation
At temperatures T much smaller than the critical temperature T c [39], a BEC is well described by the macroscopic wave function ψ = ψ(x, t) whose evolution is governed by a self-consistent, mean field nonlinear Schrödinger equation (NLSE) known as the Gross-Pitaevskii equation [34,48]. If a harmonic trap potential is considered, the equation becomes 1) where x = (x, y, z) T is the spatial coordinate vector, m is the atomic mass,h is the Planck constant, N is the number of atoms in the condensate, and ω x , ω y and ω z are the trap frequencies in x, y and z-direction, respectively. For the following we assume (w.r.o.g.) ω x ≤ ω y ≤ ω z . When ω x = ω y = ω z , the trap potential is isotropic. U 0 describes the interaction between atoms in the condensate and has the form where a is the s-wave scattering length (positive for repulsive interaction and negative for attractive interaction). It is necessary to ensure that the wave function is properly normalized. Specifically, we require A typical set of parameters used in current experiments with 87 Rb is given by

Dimensionless GPE
In order to scale the equation (2.1) under the normalization (2.3), we introducẽ with a 0 the length of the harmonic oscillator ground state (in x-direction). The coefficient of the nonlinearity of (2.6) (interaction strength parameter) can also be expressed as where x h is the healing length [11] with Thus, there are two extreme regimes: one is when ε = O(1) (⇔ a 0 = O(x s )) and κ = δε 5/2 = o(1) (⇔ 4π|a|N ≪ a 0 ), then equation (2.6) describes a weakly interacting condensate. The other is when ε = o(1) (⇔ x s ≫ a 0 ) and κ = δε 5/2 = O(1) (implying 4π|a|N ≫ a 0 ) (or ε = 1 and κ = δε 5/2 = δ with |δ| ≫ 1 by the rescaling x → ε 1/2 x, ψ → ψ/ε 3/4 ), then (2.6) corresponds to a strongly interacting condensate (Thomas-Fermi regime [11]) or, equivalently, to the semiclassical regime. We recall that the equation (2.6) is regularly perturbed in the case of weak interactions and singularly perturbed in the semiclassical regime. Analytical techniques of asymptotic analysis are available in both cases providing structural information on the solutions of (2.6) and on numerical discretization schemes (spatial/temporal grid control, control of the computational domain, error estimates in linearized cases).

Approximate ground state solution in 3d
To find a stationary state of (2.6), we write where µ is the chemical potential of the condensate. Inserting into (2.6) gives the following equation for φ(x) under the normalization condition This is a nonlinear eigenvalue problem. The Bose-Einstein condensate ground-state wave function φ g (x) is found by solving this eigenvalue problem under the normalization condition (2.11) with the minimal chemical potential µ g . Usually, the ground state problem is formulated variationally. Define the energy functional It is easy to see that critical points of E are 'eigenfunctions' of the nonlinear Hamiltonian. To compute the ground state φ g we solve In the case of a defocusing (stable) condensate the energy functional E(ψ) is positive, coercive and weakly lower semicontinuous on the unit sphere in L 2 (R 3 ), thus the existence of a minimum follows from standard theory. For understanding the uniqueness question note that E(αφ g ) = E(φ g ) for all α ∈ C with |α| = 1. Thus an additional constraint has to be introduced to show uniqueness, e.g. φ g real valued and φ g (x) > 0 for all x ∈ R 3 (see [41]). For focusing (unstable) 3-dimensional condensates the energy functional E(ψ) is not bounded from below on the unit sphere of L 2 (R 3 ). Thus, an absolute minimum of E(ψ) does not exist on {ψ ∈ L 2 (R 3 ) | ψ 2 L 2 = 1}. The interpretation of critical points (local minimum, saddle points obtained by min-maxtheory) as physically relevant ground states is by no means clear.
Using (simple) perturbation methods, we present here the approximate ground state solution of (2.6) in the two extreme regimes of weak repulsive or attractive interactions and strong repulsive interactions (see [11] for a discussion in physical literature).
These approximate ground state solutions are used in reducing the 3d GPE to a 2d GPE and a 1d GPE -see the next subsection for details -and as initial data for the numerical solution of the time-dependent GPE in section 4 (see the subsequent discussion).
For a weakly interacting condensate, i.e. ε = O(1) and κ = o(1), we drop the nonlinear term (i.e. the last term on the right hand side of (2.6)) and find the harmonic oscillator equation (2.14) The ground state solution of (2.14) is It can be viewed as an approximate ground state solution of (2.6) in the case of a weakly interacting condensate, with an O(κ)-error in approximating the chemical potential and ground state. From (2.15), we can see that the diameter of the ground state solution (computed according to the formula (4.2)) in the weakly interacting condensate is (after the scaling (2.5)).
Also we remark that the condensate widths in y and z-directions of the approximate ground state φ w Clearly, this is important for the control of the computational domain.
For a condensate with strong repulsive interactions, i.e. ε = o(1), κ = O(1) and κ > 0, we drop the diffusion term (i.e. the first term on the right hand side of (2.6)) corresponding to the Thomas-Fermi approximation [11]: The ground state solution of (2.16) is the compactly supported function φ s g : otherwise.
(2.17) This shows that the diameter of the ground state solution in the strongly interacting repulsive condensate is x s (again after the scaling (2.5)).
Approximate widths of φ s g in the y and z-directions are O (κγ z ) 1/5 /γ 4/5 respectively. This analysis suggests to choose the characteristic condensate length x s such that the (dimensionless) ground state width x w g is O(1) after the scaling (2.5), i.e. x s = O(a 0 ) for weak interaction (as in Remark 2.1) and If we use the typical set of parameter values (2.4) in the above identity, we obtain ε : 0.0078 ∼ 1.

Remark 2.2
The approximate ground state φ w g in the weak-interaction regime has finite energy, more precisely

Contrarily the energy of the Thomas-Fermi approximation is infinite
due to the low regularity of φ s g at the free boundary V (x) = µ s g . More precisely, φ s g is locally C 1/2 at the interface but not H 1 loc (R 3 ). This is a typical behavior for solutions of free boundary value problems, which indicates that φ s g does not approximate φ g to the full O(ε 2 ) -order, as indicated by formal consistency. An interface layer correction has to be constructed in order to improve the approximation quality. For a convergence proof of φ s g → φ g (without convergence rate) we refer to [41]. It is of course tempting to use approximate ground states as initial data for the GPE when simulating Bose-Einstein condensation. In the weak interaction case this produces O(κ) -errors in time dependent simulations on O(1) time intervals. In the strong interaction case an initial wave function φ s g produces time -dependent solutions with infinite energy (which usually generates breathing modes, cf. Example 3 III in Sec. 4.2) and the error in the wave function introduced by this is typically significantly larger than O(ε 2 ).

Reduction to lower dimensions
In two important cases, the 3d Gross-Pitaevskii equation (2.6) can approximately be reduced to a lower dimensional PDE. For a disk-shaped condensate with small height, i.e.
the 3d GPE (2.6) can be reduced to a 2d GPE with x = (x, y) T by assuming that the time evolution does not cause excitations along the z-axis since these have a large energy of approximatelyhω z compared to excitations along the x and y-axis with energies of abouthω x . To understand this, consider the total condensate energy E[ψ(t)]: Multiplying (2.6) by ψ t and integrating by parts show the energy conservation where ψ I = ψ(t = 0) is the initial function which may depend on all parameters ε, γ y , γ z and κ. Now assume that ψ I satisfies Take a sequence γ z → ∞ (and keep all other parameters fixed). Since R 3 |ψ(t)| 2 dx = 1 we conclude from weak compactness that there is a positive measure n 0 (t) such that |ψ(t)| 2 ⇀ n 0 (t) weakly as γ z → ∞.

Energy conservation implies
and thus we conclude concentration of the condensate in the plane z = 0: Denote by S the subspace be the projection on S:

Now write the equation (2.6) in the form
where Aψ stands for the linear part and F (ψ) for the nonlinearity. Applying Π to the GP-equation gives The projection approximation of (2.6) is now obtained by dropping the commutator terms. it reads iε(Πσ) t = ΠA(Πσ) + ΠF (Πσ), (2.26) (Πσ)(t = 0) = Πψ I , (2.27) or explicitly, with Since this GPE is time-transverse invariant, we can replace ψ 2 → ψ e −i C 2ε and drop the constant C in the trap potential. The observables are not affected by this.
The 'effective' GP-equation (2.29) is well known in the physical literature [40], where the projection method is often referred to as 'integrating out the z-coordinate'. However, an analysis of the limit process γ z → ∞ has to be based on the derivation as presented above, in particular on studying the commutators ΠA−AΠ, ΠF −F Π. In the case of small interaction ε = O(1), κ = o(1), a good choice for ψ 3 (z) is the ground state of the harmonic oscillator in z-dimension: Note that |ψ 3 (z)| 2 ⇀ δ(z) as γ z → ∞ and that ΠA = AΠ such that the error in approximating Πψ(t) by Πσ is determined by the commutator of the nonlinearity, which is O(κ).
For condensates with other than small interaction the choice of ψ 3 is much less obvious. Often one assumes that the condensate density along the z-axis well described by the (x, y)-trace of the ground state position density |φ g | 2 and (taking a pure-state-approximation) A mathematical analysis of the limit process γ z → ∞ is currently under study.
For a cigar-shaped condensate the 3d GPE (2.6) can be reduced to a 1d GPE by proceeding analogously. Then the 3d GPE (2.6), 2d GPE and 1d GPE can then be written in a unified way (2.36) By using the approximate ground state of Section 2.2, we derive -after simple calculations -for a weakly interacting condensate and for a condensate with strong repulsive interactions (2.38) We call a d-dimensional (d = 1 or 2) condensate 'very weakly interacting', if ε = O(1) and |κ d | = |κ w d | ≪ 1, which implies |δ| √ γ z ≪ 1 in 2d and |δ| √ γ y γ z ≪ 1 after reduction to 1d. Remark 2.3 By using very elongated trapping potentials it is now possible to produce weakly interacting condensates that are 'truly' in the 1D regime and the assumption that the condensate wave function factorizes is fulfilled with high accuracy. In cases where the interactions cannot be neglected low dimensional simulations [52] should be viewed as model calculations where the effective interaction strength in the reduced equation is estimated by Eq. (2.37)

Geometrical optics
where n = |ψ| 2 and S is the phase of the wave-function. Inserting into the GPequation (2.6) and separating real and imaginary parts give n t + div (n ∇S) = 0, (2.39) The equation (2.39) is the transport equation for the atom density and (2.40) the Hamilton-Jacobi equation for the phase. By formally passing to the limit ε → 0 (cf. [30]), we obtain the system It is well known that this limit process is only correct in the defocusing case κ > 0 before caustic onset, i.e. in time-intervals where the solution of the Hamilton-Jacobin equation (2.40) coupled with the atom-number conservation equation (2.39) is smooth. After the breakdown of regularity oscillations occur which make the term ε 2 2 1 √ n ∆ √ n at least O(1) such that the validity of the formal limit process is destroyed. The limiting behavior after caustics appear is not understood yet except in the one-dimensional case without confinement, see [38]. Also, the focusing case κ < 0 is not fully understood yet.

Numerical approximation
In this section we present a time-splitting Fourier spectral method, which was used by Bao et al. to numerically solve the Schrödinger equation in the semiclassical regime [7,8]. We reiterate that neither time splitting discretisations nor Fourier spectral methods are new, both have been applied successfully to many PDE problems [13,32,54]. Here we adapt the combination of both techniques to the GP equation and infer computational domain and mesh size controls from analytical (perturbation) results. The merit of this method is that it is unconditionally stable, time reversible, time-transverse invariant, and conserves the total particle number. Also, it has very favorable properties with respect to efficiently choosing the spatial/temporal grid in dependence of the semiclassical parameter ε. For simplicity of notation we shall introduce the method in one space dimension (d = 1). Generalizations to d > 1 are straightforward for tensor product grids and the results remain valid without modifications. For d = 1, the equation (2.34) with periodic boundary conditions becomes We choose the spatial mesh size h = ∆x > 0 with h = (b − a)/M for M an even positive integer, the time step k = ∆t > 0 and let the grid points and the time step be x j := a + j h, t n := n k, j = 0, 1, · · · , M, n = 0, 1, 2, · · · Let ψ n j be the approximation of ψ(x j , t n ) and ψ n be the solution vector with components ψ n j . Time-splitting spectral method (TSSP). From time t = t n to t = t n+1 , the GPE (3.1) is solved in two splitting steps. One solves first for the time step of length k, followed by solving and thus can be integrated exactly. From time t = t n to t = t n+1 , we combine the splitting steps via the standard Strang splitting: where ψ * l , the Fourier coefficients of ψ * , are defined as The overall time discretization error comes solely from the splitting, which is second order in k for fixed ε > 0. The spatial discretization is of spectral (i.e. 'infinite') order of accuracy for ε > 0 fixed. An error analysis for linear Schrödinger equations taking into account the ε-dependence of the global error can be found in [7] (where it is shown that k = O(1) and h = b−a M = O(ε) give correct observable), numerical tests for nonlinear problems in the semiclassical regimes in [8]. More restrictive meshing strategies are typically necessary in nonlinear cases, cf. Section 4.1.
For comparison purposes we review now alternative numerical methods [2,16,26] which are currently used for solving the Gross-Pitaevskii equation of BEC. One is the Crank-Nicolson finite difference (CNFD) scheme [2]: Another one is the Crank-Nicolson spectral (CNSP) method: where D f xx , a spectral differential operator approximation of ∂ xx , is defined as Both methods are unconditionally stable, time reversible, conserve the total particle number but they are not time transverse-invariant. We do not present comparism tests with fully implicit and fully explicit finite difference methods since they are not at all competitive with the time splitting-spectral method. Generally 1) they require severe stability constraints on the mesh sizes, 2) they do not conserve the total particle number, 3) they are not time transverse invariant. For a mathematical analysis of FD-methods for Schrödinger type equations in semiclassical regimes we refer to [43,44].

Numerical examples
In this section, we first perform a numerical comparison of TSSP, CNFD and CNSP in terms of accuracy and mesh size strategy for a 1d defocusing GPE. Then we apply the TSSP for solving 1d, 2d and 3d GPEs of Bose-Einstein condensation. Furthermore we also give a physical discussion on our numerical results. In our computations, the initial condition for (2.34) is always chosen in WKB form: ψ(x, t = 0) = ψ 0 (x) = A 0 (x) e iS 0 (x)/ε , with A 0 and S 0 real valued, regular and with A 0 (x) decaying to zero sufficiently fast as |x| → ∞. We compute with TSSP on a domain, which is large enough (as controlled by the initial data) such that the periodic boundary conditions do not introduce a significant aliasing error relative to the whole space problem. There are certainly more sophisticated analysis for controlling aliasing errors, however these do not significantly improve the results for exponentially decaying initial densities.
To quantify the numerical results, we define the condensate widths along the x, y and z-axis as where brackets denote space averaging with respect to the position density

Comparisons of different methods
Example 1 1d defocusing condensate, i.e. we choose d = 1 in (2.34) with positive κ 1 . The initial condition is taken as x ∈ R. First we compare the discretization error in space. We choose a very small time step, e.g., k = 0.00002 such that the error from the time discretization is negligible compared to the spatial discretization error, and solve the GPE using different methods and varying spatial mesh sizes h. Table 1 lists the numerical errors ψ(t) − ψ h,k (t) l 2 at t = 2 for varying spatial mesh sizes h. Clearly TSSP and CNSP show roughly the same errors due to the fact that the temporal discretization is almost 'exact'.
Secondly, we test the discretization error in time. Again, we take ε = 0.1 and κ 1 = 1.2649. Table 2 shows the numerical errors ψ(t) − ψ h,k (t) l 2 at t = 2 with a very small mesh size h = 1 32 for different time steps k and different numerical methods. Here, CNSP and CNFD show almost no difference since the spatial discretization now is almost 'exact'.
We also tested numerically the unconditional stability of the time-splitting spectral method, which was already proven rigorously in [7]. Numerical tests showed no significant accumulation of round-off errors and conservation of the discrete l 2 -norm was observed up to 10 significant digits for tests performed with h = 1 32 and k = 0.2, k = 0.05, k = 0.01 computing up to t = 4.
At last, we test the ε-resolution of different methods. Here we shall compare the meshing strategies required in order to get the 'correct' condensate density |ψ| 2 , for different methods when decreasing the semiclassical parameter ε. Figure  1 shows the numerical results with different combinations of ε, h, k for different methods. Furthermore Figure 2 shows the evolution of ρ = |ψ| 2 in space-time and the condensate width as a function of time by using TSSP for ε = 0.1 and κ 1 = 1.2649.
From Tables 1-3 and Figure 1, one can make the following observations: (1). For TSSP, the spatial and temporal discretization errors are of spectral and second order accuracy, respectively. The admissible meshing strategy for obtaining the 'correct' condensate density in the defocusing case is: h = O(ε) and k = O(ε). This method is explicit, unconditional stable and its extension to 2d and 3d cases is straightforward without additional numerical difficulty.
(2). For CNSP, the spatial and temporal discretization errors are also of spectral and first order accuracy, respectively. But the admissible meshing strategy is: h = O(ε) and k = o(ε). Furthermore this method is implicit and its extension to the 2d or 3d case is expensive except when an ADI technique is used, which destroys the spatial spectral accuracy.
(3). For CNFD, the spatial and temporal discretization errors are of second and first order, respectively, and the admissible meshing strategy is: h = o(ε) and k = o(ε) (see [43]). This method is implicit and the remark of (2) applies.
Example 5 2-dim vortices in BEC. We choose d = 2 and simulate the effect of stirring the (stable) condensate by adding a narrow, circularly moving Gaussian potential W (x, t) to the stationary trap potential in (2.34). W (x, t) represents, for example, a far-blue-detuned laser [14]. We set with the center x s (t) = (r 0 cos ω s t, r 0 sin ω s t) moving on a circle with radius r 0 and frequency ω s . We start the simulation with the ground state in 2 dimensions (no stirrer at t = 0) and minimize transient effects by increasing the stirrer amplitude W s (t) linearly from 0 at t = 0 to a final value W s (t = π) =: W f at t = π. The stirrer is then linearly withdrawn from t = 4π to t = 5π (after constant stirring, i.e. W s (t) = W f for π ≤ t ≤ 4π) and the condensate is left to evolve freely after t = 5π. We recall that 2-dim and 3-dim vortices simulations were already performed in [14,15,3], here we present this example in order to put our numerical method to an important physical test. We take the numerical values ε = 1/ √ 50, κ d = 1, γ y = 1, W f = √ 2, V s = 1/50 1/4 , r 0 = 2/50 1/4 and ω s = 1 for our simulation. We remark that a vortex in the condensate is a point x w with ψ(x w ) = 0 and singular or undefined phase.
In Figure 11,

Discussion
In Sec. 4.1 we compared different numerical methods for solving the GPE with the TSSP. Now we complete our investigation on the validity of using the TSSP for solving the GPE by comparing the results obtained in Sec. 4.2 with well known properties of Bose-Einstein condensates at very low temperatures.
In example 1 we present 1d simulations. Initially the condensate is assumed to be in its noninteracting ground state when at t = 0 repulsive interaction is turned on. In current experiments a change in the interaction strength can be achieved by applying external magnetic fields. Close to a Feshbach resonance the interaction strength shows a strong dependence on the magnetic field and even its sign can be changed by an appropriate choice of the magnetic field [49]. At the same time we also change the trap potential by setting ε = 0.1. These sudden changes lead to many rapid oscillations in the condensates (cf. Fig. 2 a)). The dominant excitation caused by the interactions is, as expected [37,25], the oscillation of the condensate width at approximately twice the trap frequency (cf. Fig. 2 b)).
Example 2 presents 2d simulations for various cases. In cases I there are O(1) -interactions and in case II we investigate a weakly interacting condensate. We assume the condensate to be initially in the noninteracting ground state with possibly a nonuniform phase. As in example 1 turning on the interactions causes the condensate to oscillate at approximately twice the corresponding trap frequency. Higher excitations like in Fig. 2 a) are not visible in Figs. 3 a), 4 a) since compared to example 1 the strength of the interactions is much smaller. A nonuniform phase like in II causes the amplitude of the oscillations to be different in x and y direction. Note also that while in case I the condensate immediately starts to expand (cf. Fig. 3 b)) (due to the repulsive interactions) it starts to contract for the initial condition with a nonzero phase in case II. In case III we investigate the evolution of a strongly interacting condensate initially in the (approximate) Thomas Fermi ground state with an additional phase. Again, since this is not the exact ground state (see Remark 2.2) the width of the condensate starts to oscillate at about twice the trap frequency (cf. Fig. 5 a)). In this case, however, due to the strong nonlinearity the oscillations along the x and y direction are coupled with each other as can be seen from Fig. 5 b). In case IV we investigate the effect of changing the trap frequency and turning on repulsive interactions (see Fig. 6). Like in the previous cases we find the dominant effect to be oscillations at about twice the corresponding trap frequency. The condensate initially starts to contract since the trap frequencies are increased at t = 0. For the initial conditions chosen in this case the amplitudes are sufficiently large to swap the widths σ x and σ y (cf. Fig. 7), whereas for small changes in the trap frequency no swapping of the condensate widths would occur. The numerical results for the oscillations of a BEC obtained in example 2 agree very well with experimental and theoretical results [37,25].
In example 3 we show solutions for a focusing nonlinearity in 2d. Case I shows the effect of turning on O(1) attractive interactions which leads to oscillations similar to those discussed in the previous examples (see Fig. 8). In case II a condensate with negative initial energy is shown. We have not discussed the reduction of the GPE to 2d for this case. Also, our simulations do not contain loss terms which become important in condensates at high densities. Thus we do not give a physical interpretation of these results. However, this example shows that the numerical method is applicable to the case of strong focusing nonlinearities in the GPE. Our numerical results confirm that the attractive GPE in 2d with negative initial energy will blow up at finite time (cf. Fig. 9). Furthermore, we point out that the TSSP allows the inclusion of loss terms into the GPE and it is also feasible to solve the GPE in 3d for the attractive case [10]. Therefore the TSSP is a promising candidate for simulating the recent experiments on collapsing and exploding BEC's by Donley et al. [24] which requires full 3d simulations and the inclusion of loss channels.
Example 4 shows the effects of turning on repulsive interactions and changing the trap frequency in a 3d condensate. As expected from our previous simulations we see in Fig. 10 oscillations at twice the trap frequency in directions x, y and z, respectively. The amplitude of the oscillations decreases with increasing frequency, i.e. it becomes more difficult to excite oscillations for larger trap frequencies. This behavior is one of the basic assumptions allowing the reduction of the GPE to 2d and 1d in the cases where one or two of the trap frequencies are much larger than the others (cf. Sec. 2.3).
The last example 5 shows the creation of vortices in a 2d BEC by stirring it with a blue detuned laser beam (see also [14,15,45]) where we identify the creation of three vortices in the BEC as shown in Fig. 11. We note that since we study the GPE without taking into account the interaction between the condensate and a thermal cloud of atoms (additional dissipation) no stationary state showing an Abrikosov lattice of vortices is found as in recent experiments [45] and numerical studies on the effects of a thermal cloud [47] on the vortex formation. However, we point out that full 3d simulations including the effects of thermal particles with a very high precision are feasible based on the TSSP.
Finally, we note that the TSSP is a very powerful versatile numerical method for solving the GPE which can be applied to a large number of different physical situations. The efficiency of this method and the high precision of the solutions make the TSSP a good choice for solving experimental situations that are numerically very demanding. Among these we believe that numerical studies of collapsing condensates with attractive interactions and multi-component condensates taking into account all experimentally relevant extensions of the GPE as well as on extensions of the GPE dealing with dissipation mechanisms [31] will be feasible by using the TSSP [10].

Summary
We studied a numerical method for solving the time-dependent GPE which describes trapped Bose-Einstein condensates at temperatures T much smaller than the critical condensate temperature T c . We started with the 3d GPE, scaled it to obtain a four-parameter model, and showed how to approximately reduce it to a 2d GPE and a 1d GPE in certain limits. We provided the approximate ground state solution of the GPE in two extreme regimes: (very) weakly interacting condensates and condensates with strong repulsive interactions. Then, most importantly, we used the time-splitting spectral method in connection with analytical considerations based on perturbation theory (mesh-size control, dimension reduction) to solve the time-dependent GPE in 1d, 2d and 3d. Extensive numerical examples in 1d, 2d and 3d for weakly/strongly interacting condensates defocusing/focusing nonlinearity, and zero/nonzero initial phase data were presented to demonstrate the power of the time-splitting spectral numerical method. Finally, we want to point out that equations very similar to the GPE are also encountered in nonlinear optics. In the future we plan to apply this powerful numerical method to physically more complex systems like multi component condensates, collapsing condensates with attractive interactions and also to describe coherent atomic samples in wave guides.