Tunneling of Two Interacting Fermions

. We consider two interacting atoms subject to a one-dimensional anharmonic trap and magnetic ﬁeld gradient. This system has been recently investigated by the Heidelberg group in the experiment on two 6 Li atoms. In the present paper the tunneling of two cold 6 Li atoms, initially prepared in the center-of-mass and relative motion excited state, is explored and full time-dependent simulation of the tunneling dynamics is per-formed. The dynamics is analyzed for the interatomic coupling strength ranging from strong attraction to strong repulsion.


Introduction
Description of quantum tunneling, which, although, was discovered about a century ago, still remains challenging, especially when one considers the tunneling of correlated systems. In the context of ultracold physics, there are only a few works that completely take into account the pair correlation impact in the problem of a decay of a quantum system [1][2][3][4][5][6][7][8]. In this work we present an efficient method [1][2][3][4] which completely takes into account the interatomic interaction and also an external field. This situation corresponds to two cold atoms, subject to an anharmonic optical trap and magnetic field gradient, which have been investigated in the recent experiment of the Heidelberg group on two interacting 6 Li atoms [9,10] and theoretically studied in [5,11,12]. The atoms tunnels through the potential barrier of the trap to open space. Here we focus on the excited state branch which has not been considered previously. Due to the tight confinement in the transverse direction the atomic dynamics could be analyzed by solving the one-dimensional Schrödinger equation [13], which we solve by means of a split-operator technique.

Model
Tunneling dynamics of the system considered in this work is modeled by the time-dependent Schrödinger equation where the two-body Hamiltonian reads Here are the single-particle Hamiltonians, in which V (at) describes the atom-trap interaction and m is the atomic mass. The atom-trap potential is taken from [5,9,10] which represents an optical trap and magnetic field gradient (figure 1): where µ B is the Bohr magneton, z R = 8.548ℓ (ℓ = √ /mω) is the Rayleigh range, V 0 = 56.16 ω is the maximum depth of the optical trap, p is the optical trap depth, C | j⟩ depends on external magnetic-field strength, magnetic-field gradient and the hyperfine state of the atoms, ω = 2π×1234 Hz. These values represents realistic trap parameters used in [5,9,10]. For the interatomic interaction potential V (aa) (z 1 − z 2 ) we choose the Gaussian shape [5,14] where V 0 and r 0 -are the depth and range of the interaction (5). The value r 0 = 0.1ℓ [5,14] -well represents the short-range interaction, which atoms experiences at low energies. It is more common to represent the interatomic interaction in terms of the 1D contact coupling strength g, which incorporates V 0 and r 0 via the 1D scattering length, a, as g = −2 2 /(ma). The 1D scattering length is computed by solving the scattering problem of the two atoms in the absence of the trapping potential [1,14]. In order to integrate (1) we use the split-operator method, which is based on ideas [15] and has been developed in the works [2][3][4] in application to confined ultracold atom-atom collisions in waveguide-like traps: The action of the operators exp } is approximated by the Crank-Nicolson scheme, which maintains the accuracy order of the split-operator scheme (6) to O(∆t 3 ). The partial derivatives in (3) are approximated by the sixth-order finite-differences.
To model the outgoing wave we introduce into the original Hamiltonian (2) a complex absorbing potential (CAP) iW(z j ) [16], which absorbs the wave-packet at the edges of the simulation grid: where θ(z) is the Heaviside step function. We find that the parameter values w c = −1 ωℓ −2 and z c = 25ℓ, lie in the domain where final results do not change significantly under the variation of these parameters.

Results
To prepare the initial state we fix the parameters p = 0.795 and C | j⟩ = 1894.18 G/m (for both atoms).
In this case tunneling is highly suppressed, due to the deep trap depth (figure 1) and the system is in a metastable state. In order to completely suppress tunneling and calculate eigenenergies and eigenfunctions we put a hard wall at z = 11ℓ. This procedure models the experimental realization of the initially trapped atoms [5]. To compute the spectrum we use the method from [17]. The spectrum in such a trap is shown in figure 2. We focus on the lower excited state branch. Both excited state branches correspond to doubly excited states with respect to relative and center-of-mass motions [1]. The nodal patterns of the wave functions of these states are mixed when reaching the noninteracting limit, i.e. g = 0. This mixing occurs due to the rotational symmetry breaking, which is caused by anharmonic terms of the potential (4). Similar effect has been observed in [1,18,19]. Energy spectrum of the "closed" trapping potential (4) using p = 0.795 and putting a hard wall at z = 11ℓ. Upper and lower excited state branches correspond to doubly excited states with respect to relative and center-of-mass motions [1]. Now, with the obtained initial solutions we proceed to calculation of the evolution of ψ(z 1 , z 2 , t) from (6). At t > 0 we take p = 0.68 and assume C |1⟩ = 1896 G/m and C |2⟩ = 1894 G/m. These trap values lower the potential barrier and provide the tunneling time scales to be experimentally detectable [5,9,10].
The plot ( figure 3(a)) shows the mean number (experimentally measurable quantity [9,10]) of trapped particles as a function of time calculated as where P 2 (t) is the probability to find the two particles in the box z j ∈ [−2, 13] [in ℓ units]. The probability to find one particle in the trap, P 1 (t), is calculated by integrating the density |ψ(z 1 , z 2 , t)| 2 over the regions z i ∈ [13,25] and z j ∈ [−2, 12] [in ℓ units] [5] ( figure 3(b)). The simulation grid ends at z j = 40ℓ. Configuration space for the two-atom system: R 2 is the region over which we integrate the density |ψ(z 1 , z 2 , t)| 2 to find P 2 (t); R 1 are the regions over which we integrate the density |ψ(z 1 , z 2 , t)| 2 to find P 1 (t).
Another interesting feature of the system tunneling is the behavior of the probability current |j(z 1 , z 2 , t)| (figure 4) From figure 4 we see that at strong attraction, g = −2 ωℓ, the tunneling direction dominates along the axis z 1 = z 2 and, with increasing g, tunneling goes along the axes z 1 and z 2 .
From figure 4 one can notice how the strong attraction, g = −2 ωℓ, modifies the tunneling direction. At this coupling strength, from figure 3(a), we see that the mean atom number decreases, within the considered time domain, approximately to the value 0.85, which implies that one particle remains in the trap during the tunneling. This suggests that this particle has such a small energy that its tunneling time scale is much slower than that of the first particle. Tunneling in other cases, i.e. g = 0 and g = 5 ωℓ, goes in such a way that the two particles leave the trap independently (figure 4). In this paper tunneling of two interacting 6 Li atoms from the relative and center-of-mass motion excited state is computed. This investigation is related to the recent experiments of the Heidelberg group [9,10] and extends the analysis made in [5], where this excited state branch has not been considered. A decay of the mean atom number for strong attraction (g = −2 ωℓ), no interaction (g = 0) and strong repulsion (g = 5 ωℓ) cases is demonstrated. Probability current shows different tunneling directions when varying g. The present study could be related to the investigation of pairing in BEC (Bardeen, Cooper, and Schrieffer) to Bose-Einstein crossover and to topics of shell structure evolution in the context of nuclear physics [10]. Our method performs full time-dependent simulation of atom tunneling by using split-operator and Crank-Nicolson schemes. The method could be extended to more complicated problems with more degrees of freedom.