Numerical simulation of nonequilibrium states in a trapped Bose-Einstein condensate

In this work we present numerical study of a trapped Bose-Einstein condensate perturbed by an alternating potential. The relevant physical situation has been recently realized in experiment, where the trapped condensate of 87Rb, being strongly perturbed, exhibits the set of spatial structures. Firstly, regular vortices are detected. Further, increasing either the excitation amplitude or modulation time results in the transition to quantum vortex turbulence, followed by a granular state. Numerical simulation of the nonequilibrium Bose-condensed system is based on the solution of the time-dependent 3D nonlinear Schrodinger equation within the static and dynamical algorithms. The damped gradient step and time split-step Fourier transform methods are employed. We demonstrate that computer simulations qualitatively reproduce the experimental picture, and describe well the main experimental observables.


Introduction
As is well known, trapped Bose atoms in equilibrium at zero temperature form Bose-Einstein condensate (BEC) in the ground state. This system has been thoroughly studied theoretically as well as experimentally. For details we refer to the review articles and books [1,2,3,4] and references therein. If the condensate is slightly moved from equilibrium it usually exhibits collective excitations. Such dynamics is typical not only for BEC, but for the variety of quantum systems [5,6,7].
Nonequilibrium BEC can be created by injecting into the system energy, which can be realized by applying an external modulation field. Two main ways of moving trapped BEC far from equilibrium are usually considered. The first one is a resonant excitation of coherent modes [8,9]. In this case the applied external perturbation is weak, but the modulation frequency of the alternating field is in resonance with one of the transition frequencies between topological modes. This method is a bit difficult for experimental realization, requiring fine tuning to resonance. The other way is to employ a strong external alternating field superimposed onto the stationary trap [10], where either the modulation time or excitation amplitude determine the energy injected into the system. In both these cases, BEC transfers from its ground to excited states, demonstrating a set of dynamical regimes. The second approach has been recently realized in the experiment [11,12].
The experimental realization of nonequilibrium condensates stimulates interest to its numerical modeling and computer simulation [13]. Dealing with a trapped Bose-Einstein condensate, the nonlinear Schröedinger equation (NLSE), which, in BEC physics sometimes is called the Gross-Pitaevskii equation [1], is the main tool for the theoretical investigation. Nowadays, different numerical schemes [13,14] have been developed and successfully applied for the solution of second order partial differential equations, like NLSE. The methods differ in their stability, accuracy, computational cost, and ability of keeping the NLSE original properties. In general, all such methods could be classified as explicit, implicit and finite Fourier transform or pseudospectral methods. Each numerical scheme has its own merits and disadvantages, determining the accuracy of the NLSE solution. The presence of either random forces or potentials, being a prerequisit of nonequilibrium systems, complicates the numerical procedure, since the relevant NLSE solution is not smooth.
In our recent paper [15] we presented the results of numerical simulation of a trapped condensate, under external perturbation, showing the appearance of a granular state, analogous to that observed in the experiment. The simulations were based on the solution of the NLSE, with the parameters corresponding to the experimental setup. In the present paper we give the details of the accomplished numerical modeling. Also, we demonstrate that the chosen numerical methods qualitatively describe the typical physical picture, observed in the experiment.
The paper is organized as follows. In Sec. 1, we briefly review the experiment on the strong perturbation of a trapped BEC by external modulating field. Details of computer simulation are given in Sec. 2. In Sec. 3, we present the numerically simulated time-amplitude diagram.

Experiment on perturbation of 87 Rb
In the experiment [11,12], the 87 Rb atoms (atomic mass m = 1.444 −22 g, scattering length a s = 0.557 × 10 −6 cm), being cooled down to the temperature much lower than T c = 276 nK, form the ground state Bose-Einstein condensate. The total number of atoms in the BEC fraction is N ≈ 1.5 × 10 5 . Stationary trapping potential has cylindrical symmetry, with the trap frequencies: ω r = 2π × 210 Hz, ω z = 2π × 23 Hz. After the stationary BEC state, obeying the Thomas-Fermi shape, is achieved, the trapping potential is modulated by an additional external alternating field, oscillating with the frequency Ω 0 = 2π × 200 Hz. As a result, the total trap potential V ext can be approximated as [12]: The trapping frequency ω r determines the oscillator length l r that defines the modulation parameters of the potential: If the external perturbation imposes a certain angular momentum onto the system, like in the case of the rotating [16] or laser-stirred condensate [17], the formed vortices are aligned along the rotation axis. The number of vortices depends on the modulation frequency. The aligned vortices can form a completely ordered lattice.
But in our case, as is seen from the expression (1), the modulation field used in the experiment does not impose a certain rotation axis. This feature makes the main difference of our experiment, as compared to the case of rotating BEC. Contrary to the rotated BEC, the external perturbation, we use, just shakes the condensate cloud, injecting the energy into the trap and thus generating the diverse coherent topological modes [10].
Mathematically, coherent topological modes are stationary solutions of NLSE, including the ground state and other excited states. Different modes correspond to rather different spatial distributions of the condensate density, forming a variety of spatial shapes. As a result, different structures can be generated in a trapped condensate [8,9].
First modes, generated in the experiment, are quantum vortices. Since the external field does not inject a rotation moment, the vortices appear as vortex-antivortex pairs with the winding numbers plus/minus one These modes are energetically the most stable. The number of created vortices depends on the injected energy, hence can be controlled by either the pumping amplitude or by the modulation time. At a fixed amplitude, almost a linear dependence of the vortex number from time is observed, until the distance between the vortex centers is less than four coherence lengths, when the interaction energy of two vortices is much smaller than the vortex energy. When the number of vortices becomes sufficiently large, the transition to random vortex turbulence occurs.
The continuous energy pumping into the BEC cloud creates more vortices than the trap can host, as a result of which the vortices start strongly interacting and colliding, which destroys the regime of vortex turbulence. Then the system passes to the next dynamical regime, where the condensate forms high-density grains surrounded by sufficiently rarified gas. It can be demonstrated [15] that these grains are coherent objects.

Numerical simulation
Experimental characteristics are found to be in good agreement with analytical estimates [10,15]. Now we aim at accomplishing numerical calculations for the same setup as in the experiment. In this section, we focus on the techniques of numerical modeling as applied to the solution of NLSE describing the condensate dynamics driven by external perturbations. The corresponding numerical simulation is realized in two steps: first, we calculate the parameters of the initial conditions describing the stationary state of the BEC in the trap, and then we consider time propagation. The latter is performed via the solution of the time-dependent NLSE for the desired time interval. In both these cases, the full 3D geometry is considered, and the wave function is determined on the 3D space grid. Below, the main points of the accomplished numerical simulations are illustrated.

Initial conditions
First, the wave function, describing the condensate stationary state, has to be computed. To this end, we consider the eigenproblem: where the index n enumerates quantum states and the corresponding wave functions ϕ n and energies E n . In equation (2), the nonlinear Hamiltonian H(ϕ) is typical for a trapped BEC: where g = 4πh 2 as m is the interaction strength depending on the scattering length a s and atomic mass m, V ext is the trapping potential. The parameters of the Hamiltonian (3) exactly correspond to the experimental setup. The BEC of 87 Rb atoms with the related a s and m are considered. The total number of condensed particles is N = 1.5 × 10 5 . The external field V ext is represented by the total trap potential in form (1). Since we interested in stationary solutions, the external perturbation is switched off, so that V ext ≡ V ext (r, t = 0).
The eigenproblem (2) can be solved by a variety of numerical schemes, but iterative methods are the most convenient. We used the damped gradient step method [18,19], where the initial guess φ n is improved by the iteration steps. Mathematically, this can be represented in the form: where D is a damping operator, the symbol O labels the orthonormalization of the set of functions φ n at each iteration step marked as i. The Schmidt orthonormalization is usually employed. Hamiltonian (3) is recalculated after the each iteration step. The main problem of the numerical scheme (4) is the proper choice of the damping operator D providing the convergence of the iteration procedure. The damped gradient step method operates with the damping operator in the form where δ is a parameter controlling the iteration step size, T = −h 2 2m ∇ 2 is the operator of kinetic energy. Since kinetic energy can happen to be close to zero, parameter E 0 stabilizes the iteration process.
If the procedure converges, the final set of φ n , obtained after a number of iteration steps, corresponds to the eigenfunctions ϕ n (r) with the spectra E n .
After the eigenproblem (2) is solved, the total wave function Ψ(r) can be written as a superposition of the states Ψ(r) = n C n ϕ n (r) , and then normalized to the total number of condensed atoms At t = 0, there is no external perturbation, thus the condensate is in equilibrium state corresponding to the lowest energy E 0 . To take this into account, we apply the condition C 0 >> C n>0 (8) to the notation (6). Representation (6) of the wave function corresponds to its expansion over the basis of coherent topological modes [9], which allows one to treat non-ground-state condensates [8].

Time evolution
The static solution is followed by the time propagation of the condensate wave function ψ(r). For this purpose, we solve numerically the nonlinear 3D time-dependent NLSE using the obtained Ψ(r) as the initial condition. Similarly to the previous case, the parameters of NLSE exactly correspond to the experimental setup. To solve NLSE, the time split-step Fourier transform (TSSFT) method [20,21] is chosen. The advantage of this numerical scheme [13] is that it is unconditionally stable and keeps the main properties of NLSE, such as mass conservation, time reversibility, i.e., time-inverse invariance. Also, the method is second-order accurate in both time and space. The main idea of the TSSFT is in splitting the time-propagation between the kinetic and potential components of the total energy. Both these components can be approximated by exponents related to the solutions of the linear Schrödinger equation (see Ref. [21] for details). Mathematically, the mechanism of time propagation has the following form: Ψ(r, t + ∆t) ≈ exp (iT ∆t/2) exp (−iV ∆t) exp (iT ∆t/2)Ψ(r, t) , where T = −h 2 2m ∇ 2 and V = V ext + g|Ψ(r, t)| 2 are the operators of kinetic and potential energy, respectively.
The relation between the exponential operators obeys the Baker-Cambell-Hausdorf theorem, which states that e A e B = e C is valid if C = A + B + [A, B] + .... Since the TSSFT method skips the commutation relations of the kinetic and potential operators, numerical error grows with time. The symmetrical form of (10), with respect to the operators, allows to remove the first commutator term, increasing the accuracy of time propagation [22].
The exponential derivative operator is easily calculated in the momentum space rather than in the coordinate space, because of the well known relation where symbol F means the Fourier transform, k is a variable conjugated to q. Taking into account expressions (10) and (11), the full time propagation of the condensate wave function for the time step ∆t is performed in three stages: (i) The Fourier transform of the initial function Ψ(r, t) is followed by the half time-step ∆t/2, with the operator of kinetic energy in the momentum space. The first stage is completed by the inverse Fourier transform that results in the new function Ψ(r, t) →Ψ(r, t). (ii) Recalculation of the potential operator V with the functionΨ(r, t) and further full-time step propagation ofΨ(r, t), with the recalculated term V . As a result, the transformatioñ Ψ(r, t) →Ψ(r, t) is obtained.
(iii) The final stage equals the first step, but with the functionΨ(r, t). This completes the time propagation, sinceΨ(r, t) → Ψ(r, t + ∆t) .
The time step ∆t should be small enough to keep the accuracy of time propagation. The increment ∆t fulfills the condition [23] < Ψ|H|Ψ > << 2π ∆t .

Other numerical aspects
Dealing with a grid-based technique, the parameters of the computation box are extremely important. The size of the computational box strictly depends on the trapping frequencies, total number of condensed particles and the kind of boundary conditions. The proper size of the 3D grid avoids the unphysical truncation of the operated wave function and its reflection from the grid boundaries during the perturbation. The later is important, if reflecting boundary conditions are employed. In the present simulation, we used the reflecting boundary conditions for both, the static and dynamical simulations. In the experiment, the applied external perturbation heats up the BEC leading to the transition from condensed to thermal phase for a part of atoms. This process is almost regular in time. Consequently, the number of condensed particles in the trap permanently decreases. To follow the experiment, we introduced into the system a small dissipation replacing the term ih ∂ ∂t Ψ(r, t) in the NLSE (9) by (i − γ)h ∂ ∂t Ψ(r, t), where γ is a phenomenological term. The same approach is also used by the other authors [24]. In simulations, during the 60 ms of perturbation, approximately 15 percents of the initial condensed particles are lost, which is in good agreement with the experiment. Of course, introducing the dissipation term γ does not conserve the initial norm of the wave function (7), which implies the diminishing number of condensed atoms.

Results and discussion
Computer simulations reproduce well the experimental results.
The applied external perturbation, modelled by (1) creates, first, dipole oscillations of the condensate cloud, followed by vortex structures corresponding to vortex-antivortex pairs. The vortex nature of the observed structures is confirmed by quantized circulation having the values plus or minus one. The circulation is calculated by the method proposed in Ref [25]. Increasing either the amplitude of pumping or the excitation time produces more and more vortices, creating the random vortex tangle, representing, according to Feynman [26], quantum turbulence. The next regime observed in the simulations is grain turbulence, where continuous trap modulation destroys all the vortices, forming high-density grains surrounded by rarified gas. The typical physical picture, obtained in numerical simulations, is represented by the time-amplitude diagram shown in Fig. 1. The experimental diagram, published in Ref [12], is in qualitative agreement with that obtained in the numerical simulation. The area inside the rectangle reproduces the experiment in the best way. Numerical simulations give very good quantitative values for the main characteristics observed in the experiment and also found from theoretical arguments. Thus, the critical number of vortices N cr ≈ 25, required for creating the regime of quantum turbulence, is the same in the experiment, numerical simulations, as well as in theoretical estimates. The grains can be considered as coherent objects, since their typical linear size practically coincides with the coherence length, as is found in the experiment and in the simulation.
In numerical simulations, we also observe one more regime [15], arising after the granular state, the regime of wave turbulence, when all BEC is destroyed, and turbulence is formed by small-amplitude waves. Here we compare the experimental and numerical results, because of which the regime of wave turbulence is not included in the diagram, since it has not yet been reached in experiment.