5.3 Numerical Techniques
5.3.7 Initial state generation
∂r −V(r, t)−C|Ψ|2. (5.37) In the interaction picture representation this becomes
[Dr+N]I = [Ncyl]I =ei(t−t0) ˜D[Dr+Ncyl]e−i(t−t0) ˜D. (5.38) With the additional evaluation of theDrterm using a finite differencing method, it is thus possible to use the proven and efficient code developed for cartesian 2D and 3D problems to solve cylindrically symmetric problems in “pseudo-3D”.
5.3.7 Initial state generation
The natural choice for a simulation initial state is a trap eigenstate, in particular the trap ground state. An important feature of the trap ground states is that it has a long simulation lifetime and rests in the trap. Thus it is an ideal state for code debugging and an excellent test case. Furthermore it can serve as a starting point for simulations of manipulation and excitation of condensate clouds.
Analytical methods to determine trap eigenstates using conjugate gradient techniques have been described in [133]. These methods are also suitable to calculate higher order excitation states. However, these methods are relatively difficult to implement numerically and may require further interpolation to project onto numerical grids [132]. For the sake of simplicity and at the expense of CPU computing time, we have used a method for generating eigenstates, which allows the use of the existing and proven GPE simulation code. This method will be outlined below.
Trap ground states Ψ0(r, t) evolve over time according to the equation Ψ0(r, t) = Ψ0(r)e−iµt, (5.39) whereµis the chemical potential, which can be written as
µ=
for a wavefunction normalised to unity (R
Ψ∗Ψ = 1). The energies are calculated as integrals over all spaceV as follows.
Ekin = −
0 2000 4000 6000 8000 10000 12000 14000 16000 0
0.5 1 1.5 2 2.5x 10−5
FWHM x, 15*z [m]
neg. imag. timesteps
initial state generation cylindrical ωr = 40Hz, ωz = 1000Hz
Figure 5.1: Evolution of cloud shape during the initial state calculation process using the negative imaginary time GPE propagation technique. The solid line shows the radial (x,y) FWHM, the dashed line shows thez-FWHM, multiplied by a factor of 15 for clarity. 4e5 (real) timesteps are equivalent to one full trap cycle.
Eint is the interaction energy, often called the self-energy, due to the nonlinear interaction C as defined in section 5.3.1. Ekin is the kinetic energy, which is small and may be neglected in the Thomas-Fermi approximation to obtain approximate ground states. Epotis the potential energy due to the external trap magnetic field. As we can see in eq. (5.39), the only parameter that changes during the trap eigenstate’s temporal evolution is its phase. We can use this fact to find ground states iteratively starting from an approximate solution.
Note that if we were to propagate this state in so called “negative imagi-nary time”, the e−iµt factor would turn into a homogeneous damping factor.
This can be used in a very simple and elegant technique to find the eigenstate iteratively. An initial guess, which can, for example, be the linear (gaussian) solution, is propagated over a small negative imaginary timestep and renor-malised afterwards. Repeated application of this procedure will iterate towards a state with a homogeneous flat µacross the whole cloud. Regions with larger values of µ will experience stronger damping and regions with smaller values of µwill experience weaker damping. Thereby the cloud will slowly change its shape towards the eigenstate solution. Technically this is an algorithm using the method of steepest descends [137], and in its simplicity it can be implemented with minimal changes to the RK4IP algorithm solver.
Figure 5.1 shows the evolution of the radial and axial width of a cloud during iterative damping into an eigenstate. Figure 5.2 shows cloud density profiles along y =z = 0. It is evident that the final values are reached asymptotically with a very small rate of change after long run times. The Thomas-Fermi shape is not perfectly reached after 15000 iterations. The convergence process
5.3. NUMERICAL TECHNIQUES 119
0 20 40 60 80 100 120
0 1 2 3 4 5 6 7
8x 10−3 Initial state generation crosscut y=z=0, ωr=40Hz, ωz=1000Hz
x (in numerical units)
|Ψ|2
Figure 5.2: Cloud density cuts along x-axis at iteration levels 1500, 3000 and 15000 (blue, green, red). From an initial Gaussian (not drawn), the cloud spreads out into typical inverted parabola Thomas-Fermi shape.
could theoretically be speeded up by enlarging the magnitude of the imaginary negative timestep size with increasing number of iterations. In doing this, however, care must be taken not to amplify errors due to finite machine floating point precision.
Since eigenstates generated using this method are only close to perfect, it is important to control their quality in order to judge whether or not it will be good enough for a specific simulation. Quality is usually only an issue for high pre-cision measurements where remaining residual excitations such as quadrupole modes will have disturbing effects and for extremely long running simulations.
Rapid profiling of parameter spaces with cloud eigenstates on coarse discrete grids will not be quite as demanding on eigenstate quality. Therefore, in prac-tice, small deviations are not a problem. Before “production use” eigenstates can be continuously improved running the damping program in the background on spare CPU time.
One useful measure of trap eigenstate quality can be found by mean of the virial theorem. Through the use of scaling transforms in expressions for the trap eigenstate and the fact that the eigenstate energy remains constant for small variations [126], we get the following relations (for individual dimensions iand for total values):
Ekin,i−Epot,i+1
2Eint = 0 (5.44)
2Ekin−2Epot+ 3Eint = 0 (5.45) Figure 5.3 shows how these relations approach zero during the iterative damping process (same situation as in figure 5.1). Eigenstate quality will still
0 2000 4000 6000 8000 10000 12000 14000 16000 100
101 102 103
Initial state generation cylindrical, ωr=40Hz, ωz=1000Hz
neg. imag. timesteps E [ hωr]
2Ekin−2Epot+3Eint Ekin EpotEint Etotal
Figure 5.3: Evolution of cloud energies during the initial state calcula-tion process using the negative imaginary time GPE propagacalcula-tion technique.
2Ekin−2Epot+ 3Eint exponentially approaches zero, while the component en-ergies approach a steady solution. 4e5 (real) timesteps are equivalent to one full trap cycle.
improve in terms of this sensitive measure long after noticeable changes in the cloud shape have stopped.