Numerical Solution of the Time Dependent 3D Schrödinger Equation Describing Tunneling of Atoms from Anharmonic Traps

. We present an e ﬃ cient numerical method for the integration of the 3D Schrödinger equation. A tunneling problem of two interacting bosonic atoms conﬁned in a 1D anharmonic trap has been successfully solved by means of this method. We demon-strate fast convergence of the ﬁnal results with respect to spatial and temporal grid steps. The computational scheme is based on the operator-splitting technique with the implicit Crank-Nicolson algorithm on spatial sixth-order ﬁnite-di ﬀ erences. The computational time is proportional to the number of spatial grid points.


Introduction
Pair correlation plays a major role in quantum tunneling and there are a only few works in the context of ultracold physics which accurately take into account the interatomic interaction [1][2][3][4][5][6][7]. In the present paper, we consider a tunneling of two interacting atoms from an anharmonic optical trap. We solve this problem by means of a numerical method which is based on the split-operator approximation [1-3, 8, 9], the Crank-Nicolson scheme and the sixth-order spatial finite-differences.

Method
We consider the time-dependent Schrödinger equation (throughout = ω = m = 1 units are used) with the Hamiltonian where V (aa) describes the atom-atom interaction and H j (x j ) denote the single-atom Hamiltonians in which V (at) describes the atom-trap interaction. The partial derivatives in (3) are approximated by sixth-order finite-differences.
To integrate (1) we use the split-operator method the basic ideas of which have been proposed in [8] and which has been developed in [2,3,9] in connection with the solution of confined ultracold atom-atom collisions in waveguide-like traps: The actions of the operators exp −i∆tH j (x j ) are approximated by the implicit Crank-Nicolson scheme, which maintains the accuracy order of the split-operator scheme (4) to O(∆t 3 ).
When considering a problem of particle tunneling it is necessary to correctly model the outgoing wave function at x 1 , x 2 → ±∞. This can be done by introducing into the original Hamiltonian (2) a complex absorbing potential (CAP) in the form suggested in [10]: where θ(x) is the Heaviside step function. The values of interest of the parameters w c and x c are fixed after Eq. (8).

Results
The efficiency of our method is demonstrated in the solution of the topical problem of tunneling of cold atoms [1]. The atom-trap interaction is chosen of the standing-wave form: where β = 2π/λ · √ /(mω) (in standard units), λ and ω represent the trap wavelength and frequency correspondingly. For β = 0.427471, this models an actual atom-trap interaction [1,11,12] in modern experiments [13] on ultracold atoms (figure 1). The atom-atom interaction is taken of Gaussian form: where the range of the interaction is fixed to r 0 = 0.1 [4,14] and the interaction strength V 0 is varied. We calculate the total probability P(t) in the box |x 1 , x 2 | ≤ 20. For a wide range of values V 0 , P(t) decays exponentially. This allows us to extract the tunneling rate γ as We find that the parameters of the CAP (5), w c = −0.1 and x c = 10, lie in the domain where the tunneling rate γ does not change significantly under the variation of these parameters.
To solve (1), the value of the initial wave function is to be known. We fix it by replacing the trap (6) with the "closed" sextic (i.e., the power series expansion of (6) up to x 6 j ) trapping potential [1,15] and calculating the eigenstates of the doubly excited relative and center-of-mass motion states of the stationary Schrödinger equation by using the numerical scheme from [12]. Figure 2 shows the fast convergence of the errors and tunneling rates γ with the spatial ∆x and temporal ∆t steps. In figure 3 one can see that the computational time grows linearly with the number of spatial grid points N, according to the sweeping method theory [16].

Conclusion
We have solved a tunneling problem of two cold interacting atoms by using a numerical method. The interatomic coupling strength is varied in a wide range -from strong attraction to strong repulsion. Figure 3. Dependence of the CPU time (in minutes) for the derivation of the tunneling rate γ (open circles) at V 0 = 9 (doubly excited relative motion state at t ≤ 0) on the number of spatial grid points N = N 1 × N 2 (N j defines the grid along x j ). The solid line is a fit function. The computation has been performed until the temporal variable reaches the value t = 10 using the server at the Bogoliubov Laboratory of Theoretical Physics, JINR (i7a: 8-core Intel E5-2690 @2.9-3.8GHz, 128 GB ECC RAM, GPU GTX680).
The method is based on the split-operator approximation of the time-evolution operator and Crank-Nicolson approximation of the single-particle Hamiltonians. We have found that the error of the method reduces quite fast with reducing the step sizes of the spatial and temporal grids. The method can be extended to more complicated problems with more degrees of freedom.