Tunneling of two bosonic atoms from a one-dimensional anharmonic trap

TL;DR

This paper models the quantum tunneling behavior of two interacting bosonic atoms in a one-dimensional anharmonic trap, revealing sequential tunneling and spectrum rearrangements depending on interaction strength.

Contribution

It provides a detailed analysis of tunneling rates as a function of interatomic coupling, linking theoretical calculations with experimental observables in a novel way.

Findings

Tunneling occurs sequentially for the considered states.

Spectrum rearrangement occurs in the noninteracting limit.

Tunneling rate depends on the effective coupling constant g.

Abstract

We investigate the quantum dynamics of two interacting bosonic atoms confined in a one-dimensional anharmonic trap. The tunneling rate, an experimentally measurable parameter of the system, was calculated as a function of the effective coupling interatomic constant $g$ from the ground (n=N=0) and first excited atomic states in the trap with respect to relative (n=2,N=0) and center-of-mass (n=0,N=2) atomic motion. This allows to investigate the initial population and pair correlation, as well as the effective coupling constant $g$, of the system by comparing the calculated tunneling rate with the experimental one. We have observed that the only possible tunneling scenario is a sequential particle tunneling in the cases we considered. We have also analyzed a rearrangement $(0,2)\rightleftarrows(2,0)$ of the spectrum in the limit $g\rightarrow \pm 0$ of noninteracting atoms.

Equations

$$H=$$

$$V_{\text{int}}(x_1-x_2)=-V_0\exp \left\{-\frac{(x_1-x_2{)}^2}{2r_0^2}\right\}$$

$$V(x_j)=V_{\text{d}}{\sin }^2\left(\frac{2\pi }{\lambda }{x}_j\right),j=1,2.$$

$$V(x_j)=-\frac{\hslash \omega }{12\alpha }{\sin }^2\left(\sqrt{-6\alpha }\frac{x_j}{\ell }\right),j=1,2,$$

$$V^{(\text{sw})}(x_j)=$$

$$i\hslash \frac{\partial \psi (x_1,x_2,t)}{\partial t}=H(x_1,x_2)\psi (x_1,x_2,t)$$

$$\psi (x_1,x_2,t+\Delta t)=\exp \left\{-i\frac{\Delta t}{2\hslash }{V}_{\text{int}}(x_1-x_2)\right\}$$

$$\times \exp \left\{-\frac{i\Delta t{H}_1(x_1)}{\hslash }\right\}\exp \left\{-\frac{i\Delta t{H}_2(x_2)}{\hslash }\right\}$$

$$\times \exp \left\{-i\frac{\Delta t}{2\hslash }{V}_{\text{int}}(x_1-x_2)\right\}$$

$$\times \psi (x_1,x_2,t)$$

$$H_j(x_j)=-\frac{{\hslash }^2}{2m}\frac{{\partial }^2}{\partial x_j^2}+V^{(\text{sw})}(x_j),j=1,2.$$

$$\exp \left\{-\frac{i\Delta t{H}_j(x_j)}{\hslash }\right\}={\left(1+\frac{1}{2\hslash }i{H}_j\Delta t\right)}^{-1}$$

$$W(x_j)=w_c(|x_j|-x_c{)}^2\theta (|x_j|-x_c),j=1,2,$$

$$P(t)=\underset{-x_m}{\overset{x_m}{\int }}\underset{-x_m}{\overset{x_m}{\int }}d{x}_{1}d{x}_2|\psi (x_1,x_2,t){|}^2\sim \exp \left\{-\gamma t\right\},$$

$$\gamma =-\frac{1}{P(t)}\frac{dP(t)}{dt}.$$

$$g=-\frac{2{\hslash }^2}{m{a}_{\text{1D}}}.$$

$$\left[-\frac{{\hslash }^2}{2\mu }\frac{d^2}{d{x}^2}-V_0\exp \left\{-\frac{x^2}{2r_0^2}\right\}\right]{\psi }_{\text{sc}}(x)=\frac{{\hslash }^2{k}^2}{2\mu }{\psi }_{\text{sc}}(x),$$

$${\psi }_{\text{sc}}(x)\underset{x\to \pm \infty }{\to }\cos (k|x|+\delta (k))$$

$$a_{\text{1D}}=\underset{k\to 0}{\lim }\frac{\cot (\delta (k))}{k}.$$

$$V^{(6)}(x_j)$$

$$j_k(x_1,x_2,t)=\frac{\hslash }{2mi}\left({\psi }^{\ast }\frac{\partial \psi }{\partial x_k}-\psi \frac{\partial {\psi }^{\ast }}{\partial x_k}\right),$$

$$k=1,2$$

$$P_k(t)=\underset{regionk}{\iint }d{x}_{1}d{x}_2|\Psi (x_1,x_2,t){|}^2,$$

$$k=1,2,3,4.$$

$$P_{\text{fit}}=P_0{e}^{-{\gamma }_{1}t}$$

$$P_L(t)=|⟨\psi (x_1,x_2,t)|{\psi }^{(L)}(x_1,x_2)⟩{|}^2$$

$$P_U(t)=|⟨\psi (x_1,x_2,t)|{\psi }^{(U)}(x_1,x_2)⟩{|}^2$$

$$P_G(t)=|⟨\psi (x_1,x_2,t)|{\psi }^{(G)}(x_1,x_2)⟩{|}^2,$$

Full text

III I. Introduction

One of the basic effects of quantum mechanics - a particle tunneling through a repulsive barrier, responsible for such fundamental processes as alpha decay and nuclear fission and fusion, has attracted in recent years a great attention in connection with the cold atom simulation of different phenomena from solid state to nuclear and high energy physics. Particularly, in the recent experimentszuern such aspects of the atomic tunneling trough the walls of confining traps as pairing and BCS-BEC crossover (from a Bardeen-Cooper-Schrieffer pairing to a Bose-Einstein condensate) have been investigated. To this class of problems can also be attributed a tunneling of BECbeinke and ultracold bosonic few-body systemslode ; lode1 , a transport of the repulsive BEC and a modeling of the Josephson effect in a double-well trapnesterenko . In the workszuern ; zuern2 it was shown that the tunneling rate through the walls of atomic traps is an experimentally measurable parameter containing an important information about the atomic dynamics inside the trap as well as the initial state of the quantum system. However, to extract this information one has to perform a corresponding accurate calculation of the tunneling rate for comparison with the experimental one.

So far, a theoretical description of the tunneling dynamics through repulsive barriers of different form is quite non-trivial task. In the modern computations a semi-classical approach of Wentzel-Kramers-Brillouin (WKB) remains the basic analytical method despite the known shortcomingsrontani ; rontani2 . The main disadvantage here is that the WKB method completely neglects interparticle interactions and therefore can produce significant errors in the end resultsgharashi . Therefore, to treat different tunneling dynamics, which depends on the specific peculiarities of each problem, a variety of numerical approaches was developed during the last two decades in atomic, molecular and nuclear physicsmelezhik ; krass ; gusev ; gharashi ; maruyama ; scamps .

In the present work we investigate the tunneling dynamics of two interacting bosonic atoms through the walls of a one-dimensional (1D) anharmonic trap by using an extension of the computational splitting-up technique suggested in melezhik for ionization of hydrogen-like atoms by strong electric fields. With this approach we calculate the dependence of the tunneling rate on the effective coupling constant $g$ from the first three low-lying atomic states in a confining trap. The rates of the transitions between the states are also investigated. The obtained results can be used to recover the physical picture inside the confining trap by comparing the calculated tunneling rates with the experimental ones. Similar tunneling processes were qualitatively investigated for tunneling through a box-shaped potential model from the ground state of a rectangular potential wellhunn .

The paper is organized as follows. In Section II we define the Hamiltonian of two-atomic system confined in 1D anharmonic trap. Key points of the splitting-up method are given in Section III for numerical integration of the 2D time-dependent Schrödinger equation describing two-body quantum dynamics in a 1D anharmonic atomic trap. Special attention is paid here for the stable and accurate procedure of extracting the tunneling rate of the system. The obtained results of calculation of the tunneling rates and transition probabilities are given and discussed in Section IV. Here, we also discuss the rearrangement of the spectrum of the confined two atoms in the limit $g\rightarrow\pm 0$ of noninteracting atoms. Finally, in Section V we draw our conclusions and provide a short outlook.

IV II. Problem formulation

The quantum dynamics of two identical bosonic atoms with masses $m$ in 1D confining trap $\sum\limits_{j=1,2}V(x_{j})$ is described by the following Hamiltonian:

[TABLE]

where the interatomic potential $V_{\textrm{int}}(x_{1}-x_{2})$ is chosen in the Gaussian form

[TABLE]

with the depth $V_{0}$ and $r_{0}$ defining the range of the interaction.

The Hamiltonian (2.1) can be considered as an effective Hamiltonian describing dynamics of two atoms tightly confined in the transverse direction $(y,z)$ (atomic motion is forbidden in the transverse direction) but with a softer confinement in the longitudinal direction $x$ described by a standing-wave form peng ; bloch :

[TABLE]

The interaction of the atom $j$ with the optical trap (2.3) is defined by the wavelength $\lambda$ of the external laser field and the atomic polarizability (included in the parameter $V_{\textmd{d}}$) peng ; bloch . Here we use the parametrization of the workspeng ; ishmukh

[TABLE]

where the parameter of anhamonicity $\alpha=-\dfrac{8\pi^{2}\hbar}{12\lambda^{2}m\omega}$, and $\omega$ and $\ell$, defined as $\omega=\dfrac{2\pi}{\lambda}\sqrt{\frac{2|V_{\textmd{d}}|}{m}}$ and $\ell=\sqrt{\frac{\hbar}{m\omega}}$, were introduced.

To have a realistic scale for the atom-trap interaction (2.4) we use the parameters $\lambda$ and $\omega$ corresponding to the optical traps from the experimenthaller , where the confined 133Cs atoms were investigated (see Table 1).

[TABLE]

Table 1. The trap parameters from the Innsbruck experimenthaller .

In the present work we restrict ourselves by the consideration of the atomic dynamics in the single-well of the 1D lattice (2.4) by approximating the latter as

[TABLE]

Such approximation neglects the tunneling of the atoms through the neighbour walls as well as the reflection from the walls (see Fig.1)

Fig.1: (Color online) Approximation (2.5) (gray solid line) of the exact trap potential (2.4) (blue dashed line) for $\alpha=-0.0304552$.

V III. Method

To calculate a tunneling rate, $\gamma$, from the bound state of the potential $V_{\textrm{int}}(x_{1}-x_{2})+\sum\limits_{j=1,2}V^{(\textsf{sw})}(x_{j})$ we integrate the 2D time-dependent Schrödinger equation (SE)

[TABLE]

with the Hamiltonian $H(x_{1},x_{2})$ defined by Eqs.(2.1), (2.2) and (2.5). Based on ideas from marchuk , which were developed in the worksmelezhik2 ; melezhik3 in application to confined ultracold atom-atom collisions in waveguide-like traps, we employ the component-by-component split-operator method to integrate Eq.(3.1):

[TABLE]

where

[TABLE]

The computational scheme (3.2) is correct up to terms of the $\mathcal{O}(\Delta t^{3})$ order. Following melezhik2 ; melezhik3 , we approximate the action of the differential operators $\exp\left\{-i\Delta t/\hbar H_{j}(x_{j})\right\}$ by implicit Crank-Nicolson scheme

[TABLE]

which maintains the accuracy of the split-operator method (3.2).

Finiteness of the width of the confining potential wall causes broadening of the energy levels in the potential describing interatomic and atom-trap interactions $V_{\textrm{int}}(x_{1}-x_{2})+\sum\limits_{j=1,2}V^{(\textsf{sw})}(x_{j})$ due to the atom tunneling through the wall of the confining trap. This means that we need to have an outgoing wave away from the region of action of the confining potential, i.e. at $x_{1},x_{2}\rightarrow\pm\infty$ meyer ; maruyama . This kind of boundary condition in the time-dependent scheme can be modeled by introducing at the edge of the radial grid $x_{m}$ some type of an absorber meyer or a mask function kulander . Here we choose the scheme with the additional complex absorbing potential (CAP) near the edge of the radial grid in the form suggested in meyer ; scamps

[TABLE]

where $\theta(x)$ - is the Heaviside step function and the parameter $x_{c}$ defining the region where the CAP switches on and it should be chosen at the point behind the barrier of the confining potential $V^{(\textsf{sw})}$. The choice of the parameter $w_{c}$ is discussed in the paragraph after Eq.(3.6).

A plot of CAP (3.4) with the confining potential $V^{(\textsf{sw})}$ (2.5) is shown in Fig.2.

Fig.2: (Color online) Plot of the absorbing potential $W(x)$ (3.4) (thick red solid line) and the confining potential $V^{(\textsf{sw})}(x)$ (2.5) for $\alpha=-0.0304552$ (blue solid line). We also plot the wave function, $\psi_{n=2}(x)$ (black dashed line) of the first excited state in the confining potential for illustrating the scale of the problem.

The numerical integration of the Schrödinger equation (3.1) with the Hamiltonian $H(x_{1},x_{2})=\sum\limits_{j=1,2}(H_{j}(x_{j})+iW(x_{j}))+V_{\textrm{int}}(x_{1}-x_{2})$ defined by (3.3), (2.5), (3.4) and (2.2) permits to extract the desired tunneling rate $\gamma$ (or the energy level width $\Gamma=\hbar\gamma$) from the decay of the total probability

[TABLE]

to find the atoms in the box $|x_{1},x_{2}|\leq x_{m}$, i.e. the total population of the atomic bound states in the box. From Eq.(3.5) one can define the tunneling rate as

[TABLE]

The determination of the tunneling rate from (3.6) is obviously holds only for the exponential decay of the probability (3.5). This condition restricts the choice of the parameters $x_{c}$ and $w_{c}$ of CAP (3.4) as well as the time domain where the decay of the norm (3.5) is stabilized after the beginning of the tunneling melezhik . Here we choose $x_{c}=10\ell$ and $w_{c}=-0.1\hbar\omega\ell^{-2}$ so that the tunneling rate, $\gamma$, remains constant to a good accuracy.

Following the pioneering worksolshanii ; busch , laid the foundation for investigations of the confined two-body systems in quasi-1D geometry of atomic traps, we define here the interatomic interaction through the effective coupling constant $g$ connected with the 1D scattering length $a_{\textmd{1D}}$ as

[TABLE]

The scattering length $a_{\textmd{1D}}$ was calculated by the integration of the 1D Schrödinger equation

[TABLE]

describing the atom-atom collision in 1D free-space, with the boundary condition

[TABLE]

at zero energy limit $k=\sqrt{2\mu E}/\hbar\rightarrow 0$. The calculated scattering phase $\delta(k)$ defines at $k\rightarrow 0$ the 1D scattering length

[TABLE]

Here, $\mu=m/2$ and $x=x_{1}-x_{2}$ are the reduced mass and the relative coordinate of the atomic pair, respectively.

The dependence of the coupling constant $g$ on the depth $V_{0}$ of the interaction potential (2.2), calculated at $r_{0}=0.1\ell$ is shown in Fig.3. The choice of the parameter $r_{0}$, fixed in ishmukh2 ; gharashi , adequately corresponds to current experiments zuern ; zuern2 ; haller where the range of confining potential $\ell$ always essentially exceeds the range of interatomic interaction: $r_{0}\ll\ell$.

Fig.3: (Color online) Dependence of the effective coupling constant $g$ on the depth $V_{\textmd{0}}$ of the interaction potential (2.2) (blue solid line) at $r_{0}=0.1\ell$.

In the calculations the parameter $V_{0}$ varied in the range $-64\lesssim V_{0}/(\hbar\omega)\lesssim 53$ which corresponds to $-8\lesssim g/(\hbar\omega\ell)\lesssim+\infty$. It has permitted to investigate the tunneling processes in a wide range of the coupling strength, $g$, from an attraction $g<0$ to a strong repulsion $g\rightarrow+\infty$.

IV IV. Results

IV.1 4.1. Preparation of the initial state

The modern experimental set up permits preparation of the well defined and practically non-decaying initial atomic states in confining traps with the subsequent “switch on” of the tunneling process by means of narrowing the width of the confining potentialzuern ; zuern2 . To model such processgharashi , first, we prepare the non-decaying initial atomic bound state at $t\leq 0$ by solving the eigenvalue problem for the potential $V_{\textrm{int}}(x_{1}-x_{2})+\sum\limits_{j=1,2}V^{(\textsf{6})}(x_{j})$ with the confining trapgrish ; sala

[TABLE]

having infinite width of the walls and repeating the form of the internal part of the confining potential (2.5) $V^{(\textsf{sw})}(x_{j})$. At $t>0$, the trap $V^{(6)}(x_{j})$ is replaced by $V^{(sw)}(x_{j})$ to allow the atoms to tunnel out of it.

A plot of the potentials $V^{(\textsf{sw})}(x_{j})$ and $V^{(6)}(x_{j})$ is presented in Fig.4.

Fig.4: (Color online) “Initial” confining trap $V^{(6)}(x)$ (red solid curve), that confines the motions of the atoms at $t\leq 0$ and a modified trap $V^{(\textsf{sw})}(x)$ (blue dashed curve), which releases the atoms due to the quantum tunneling at $t>0$.

Fig.5: (Color online) Three lowest energy levels of two bosonic atoms as a function of the coupling constant $g$ in the confining potential $V^{(6)}(x_{j})$ with $\alpha=-0.0304552$ (blue curves) and in the harmonic trap $V^{(2)}(x_{j})=\frac{1}{2}m\omega^{2}x_{j}^{2}$ (gray curves).

In Fig.5 we present the energy spectrum of two bosonic atoms calculated as a function of the coupling constant $g$ in the confining potential $V^{(6)}(x_{j})$ (4.1) and in the harmonic trap $V^{(2)}(x_{j})$ corresponding the case $\alpha=0$.

For identification of the calculated states of the spectrum we use quantum numbers $n$ and $N$ characterizing the quantization (induced by the trap) of the relative atomic and the center-of-mass motions, respectively, which are good ones in the harmonic limit (i.e. when $\alpha=0$ in (4.1)) due to the separation of the relative atomic and center-of-mass variables.

In Fig. 5 the three lowest states of the calculated spectrum are presented: the ground state ($n=0,N=0$) and the first two excited states ($0,2$) and ($2,0$).

When $\alpha=0$ the energy branches $(0,2)$ and $(2,0)$ cross each other at zero coupling, $g=0$ (non-interacting atoms). This corresponds to a pure two-dimensional harmonic oscillator and these levels become degenerate due to rotational symmetry. This symmetry breaks if $\alpha\neq 0$ and we observe an avoided crossing of these energy levels at $g=0$.

These states rearrange when they cross the point $g=0$, that is, the state $(2,0)$ turns into the state $(0,2)$ and, vice versa, the state $(2,0)$ turns into the $(0,2)$ state (see Section 4.3). The quantum numbers (n,N) of the two-atomic state inside of the trap we define by the nodal structure of the initial wave function $\psi_{n,N}(x,y,t=0)$ with respect to variables of the relative motion $x=x_{1}-x_{2}$ and the center-of-mass $y=(x_{1}+x_{2})/2$. The analysis in Sections 4.3 shows that the first excited state at negative $g$ is (0,2) and the second one is (2,0) (see Fig.5). At positive $g$, the first excited state becomes (2,0) and the second one - (0,2) due to the rearrangement $(0,2)\rightleftarrows(2,0)$ of the spectrum in the limit $g\rightarrow\pm 0$.

IV.2 4.2. Tunneling dynamics from bound states of the two-atomic confined system

By numerical integrating the 2D time-dependent Schrödinger equation (3.1) for $t>0$ we calculate the time evolution of the two-atomic wave-packet (3.2) from the ground ($0,0$) and excited states ($0,2$) and ($2,0$) as a function of the coupling constant $g$. In the calculation the sixth-order finite-difference approximation on the uniform spacial grid over $x_{1}$ and $x_{2}$ was usedmelezhik2 ; melezhik3 . The range, $|x_{j}|\leq x_{m}$, of the space of radial variables was chosen as $x_{m}=20\ell$ and the step of integration over radial variables as well as the step of integration over the time $\Delta t=0.01\omega^{-1}$ were chosen to keep the accuracy of the calculation of the tunneling rates within the order of one percent. In Fig.6 we present the calculated tunneling rates $\gamma$ defined by (3.6) from the first two excited states ($0,2$) and ($2,0$) for a wide range of the strength of the coupling constant $g$.

Fig.6: (Color online) (a) Tunneling rates from the first two excited states as a function of the coupling constant $g$ for $\alpha=-0.0304552$. (b) The binding energies of the first two excited states as a function of the coupling constant $g$ for $\alpha=-0.0304552$.

Upper energy branch of the excited states

Fig.7 demonstrates the calculated time-evolution of the total probabilities $P(t)$ (3.5) decaying from the upper branch of energy levels for excited states presented in Fig.6,b and the corresponding tunneling rates $\gamma(t)$ (3.6). This figure demonstrates quite fast transition of the decay of the total probabilities $P(t)$ to the exponential law as well as fast convergence of the tunneling rates $\gamma(t)$ to the limiting value $\gamma(t)\rightarrow\gamma(\infty)$ for a wide range of the coupling constant $g$.

To understand a mechanism of tunneling from the upper energy branch of the excited states we also calculate the probability flux $j_{k}(x_{1},x_{2},t)$

[TABLE]

at certain time, $t$. In Fig.8 the fluxes $|\textbf{j}(x_{1},x_{2},t)|$, calculated for different values of the coupling constant $g$, are presented.

Fig.7: (Color online) The time-evolution of the total probability $P(t)$ (black dashed curves) from the upper branch of energy levels for excited states presented in Fig.6,b and the corresponding tunneling rates $\gamma$ (blue solid curves). Calculations were performed for a few values of the coupling constant: $g/(\hbar\omega\ell)=-6,-1,0,3,20$.

[TABLE]

Fig.8. Distribution of a modulus of the flux $|\textbf{j}(x_{1},x_{2},t)|$ (in $\omega\ell^{-1}$ units) from the upper branch of the excited states presented in Fig.6 and its direction (in arbitrary units) for different values of the coupling constant $g$ at $t=120\omega^{-1}\approx 1.32$msec.

From Fig.8 we see that the flux in all of the cases shows a complicated behavior near the origin $x_{1}=x_{2}=0$. What all graphs in Fig.8 have in common is that in all considered cases the dominate flux directions are the directions along the axes $x_{1}$ and $x_{2}$ what corresponds to a single-particle tunneling.

To analyze the mechanism of the tunneling more quantitatively we divide (following the work maruyama ) the whole space of radial variables ($x_{1},x_{2}$) into several regions (see Fig.9).

Fig.9: Division of the radial variable space into several regions for the partial probabilities $P_{k}$ analysis.

and calculate partial probabilities, $P_{k}$, to find the atom or atomic pair in each region $k$

[TABLE]

Detection of the atom in the regions $3$ corresponds to a single-particle tunneling from the trap and in the regions $1$ corresponds to a tunneling of two particles as a bound system. The situation in which a particle tunnels to the regions $2$ was observed in maruyama . We define the size of the region $4$ to approximately cover the initial atomic distribution and take $R=5\ell$.

The calculated time evolution of the partial probabilities $P_{k}(t)$ is shown in Fig.10.

Fig.10: (Color online) Partial probabilities $P_{k}(t)$ (4.3) to populate the regions $k$ for three different values of the coupling constant $g$.

Fig.10 demonstrates that the partial probabilities - $P_{3}(t)$ - are a few order of magnitude larger than $P_{1}(t)+P_{2}(t)$, which are practically negligible in all of the considered cases. This clearly shows that the sequential particle tunneling is the dominating mechanism of the tunneling.

Lower energy branch of the excited states

To understand the dynamics of a tunneling process of the atoms bound initially in the lower energy branch of the excited states (Fig.6,b) we analyze in details the case $g/(\hbar\omega\ell)=5$ since it captures all the features of such tunneling (see Fig.11).

Fig.11: (Color online) The tunneling rate $\gamma$ (blue solid curve) from the lower energy branch of the excited states (shown in Fig.6,b) and the corresponding total probability $P(t)$ (black dashed curve) at $g/(\hbar\omega\ell)=5$ and $\alpha=-0.0304552$. The inset shows a close-up view of the $\gamma$ and $P(t)$.

From Fig.11 we clearly see that the total probability $P(t)$ decreases in two stages, i.e. we see that the initial exponential behavior of $P(t)$ decreases as $\exp\left\{-\gamma_{1}t\right\}$ which turns to, approximately at $1.5-2.5$ msec, the same exponential law, but now with different tunneling rate $\gamma_{1}\rightarrow\gamma_{2}$. The tunneling rate $\gamma_{1}$ calculated with Eq.(3.7) highly oscillates during the first stage of the decay and after approximately $\simeq 1.5-2.5$ msec, when the second stage of the decay becomes dominating, the oscillations in $\gamma_{2}(t)$ significantly damp out.

To extract the tunneling rate $\gamma_{1}$ we fit the total probability $P(t)$ to the exponential function at the time-interval $t\leq 1.5-2.5$ msec (Fig.12)

[TABLE]

Fig.12: (Color online) The total probability $P(t)$ (black dashed curve) and a fitting function (4.4) (red dashed line) in a logarithmic scale. The tunneling rate extracted from the fitting function is about $\gamma_{1}\simeq 4520$sec*-1* in the region $t\leq 1$ msec.

By using this fitting procedure we extract the tunneling rate $\gamma_{1}$ of the first decay stage for a wide range of coupling constant $g$. The result of calculation of the tunneling rate $\gamma_{1}$ is plotted in Fig.6,a with open circles.

The tunneling rate, $\gamma_{2}$, of the second stage of the decay from the lower energy branch of the excited states converges to a constant value with growing the time much better than $\gamma_{1}$ but not as good as the tunneling rates from the upper energy branch given in Fig.7, which can be noticed from the inset of Fig.11. If we fit the total probability $P(t)$ with the exponential function in a window $2-10$ msec we get for the $\gamma_{2}\simeq 580$ sec*-1*. This value approximately coincides with the value of tunneling rate from the upper energy branch which equals $\simeq 592$ sec*-1* at $g/(\hbar\omega\ell)=5$ (see tunneling rates in Fig.6,a labeled with closed circles). To understand this effect we have calculated the time-evolution of the populations of the first three low-lying bound states of the two-atomic confined system (see Fig.13) by formulaes

[TABLE]

where $\psi^{(L)},\psi^{(U)}$, and $\psi^{(G)}$ are the wave functions of the two-atomic bound states corresponding to the lower and upper energy branches of the excited states and the ground-state in confined geometry of the trap (4.1).

Fig.13: (Color online) Time-evolution of the populations $P_{L}(t),P_{U}(t)$, and $P_{G}(t)$ and the total probability $P(t)$ (black solid curve) to find atoms in the trap calculated at $g/(\hbar\omega\ell)$.

From Fig.13 one can notice that after approximately $1.5$ msec, the population $P_{U}(t)$ of the upper energy branch becomes dominating due to transition from the lower energy branch. That is, after about $1.5$ msec the tunneling occurs from the upper energy branch and therefore the value $\gamma_{2}$ of the tunneling rate approaches to the tunneling rate from the upper energy branch calculated in the previous subsection.

With increasing time to about $10$ msec the populations of the upper energy branch and the ground states become comparable (see Fig.13) and with further increasing of time the system passes to the ground state where the tunneling rate naturally defined by the decay of the ground state. It is interesting that in the case of very strong interatomic coupling $g=\infty$ the second stage goes with the tunneling rate which approximately coincides with the tunneling rate from the initial ground state (Fig.14) due to more fast population of the ground state than the upper energy branch. To see that, one has to compare the inset of Fig.14 with the converged result for the tunneling rate $\gamma$ from the ground state given in Fig.16.

Fig.14. (Color online) The tunneling rate, $\gamma$, (blue solid curve) from lower branch of the excited state (shown in Fig.6,b) at $g\rightarrow\infty$ and the corresponding total probability $P(t)$ (black dashed curve). The inset shows the convergence $\gamma(t)\rightarrow\gamma(\infty)$ to the value which approximately coincides with $\gamma$ from the ground state.

Ground state

The tunneling rate from the ground state behaves monotonically and it is significantly smaller in magnitude than those from the excited states considered above (Fig.15.a). Moreover, the tunneling is significantly suppressed even at infinite coupling constant $g=\infty$ (Fig.16). This is understood due to considerable increase of the width of the confining potential for the ground state in comparison with excited states.

Fig.15: (a) The tunneling rate $\gamma$ from the ground state as a function of the coupling constant $g$ and (b) the corresponding binding energy.

Fig.16: The tunneling rate $\gamma$ (blue solid curve) and the total probability $P(t)$ (black dashed curve) from the ground state at $g\rightarrow\infty$.

IV.3 4.3. Spectrum rearrangement

In this Subsection we analyze the spectrum rearrangement which occurs for our anharmonic trap $V^{\textsf{sw}}$ at the transition of the special point $g=0$ of non-interacting atoms with increasing $g$ from small negative to small positive values or vice versa. In Fig.5 we labeled the upper and lower energy branches at the negative side of the coupling constant $g$ with quantum numbers ($2,0$) and ($0,2$) correspondingly. These quantum numbers conserve only for the negative side of the coupling $g$. When these branches cross the points $g=0$ the nodal structures of these states rearrange and the states - ($2,0$) and ($0,2$) - interchange between each other.

With decreasing the anharmonic parameter $\alpha$ (see Fig.17) the effect becomes more pronounced.

Fig.17: The confining potential $V^{(6)}(x)$ for two values of the anharmonic parameter $\alpha$: $\alpha=-0.0304552$ (red dashed line) and $\alpha=-0.01$ (blue solid line).

Fig.18 shows the calculated energy levels of the pair ($2,0$) and ($0,2$) of the first excited states for $\alpha=-0.01$.

Fig.18: Energy levels of the pair (2,0) (closed circles) and (0,2) (open circles) of the first excited states for the confining potential $V^{(6)}(x_{j})$ with $\alpha=-0.01$.

The corresponding wave functions for the two cases - $\alpha=-0.01$ and $\alpha=-0.0304552$ of the upper and lower branches, given in Fig.18 and Fig.5 correspondingly, are presented in Fig.19 and Fig.20, which are plotted as functions of relative $x=x_{1}-x_{2}$ and center-of-mass $y=(x_{1}+x_{2})/2$ variables.

From Figs.19 and 20 we see how the wave functions change their nodal structures when they cross the point $g=0$: the nodal structure of the upper branch’s eigenfunction interchange between ($2,0$) and ($0,2$) states; similar effect occurs for the nodal structure of lower branch’s eigenfunction, only vice versa. The similar effect of the spectrum rearrangement for two atoms confined in 3D anharmonic trap was observed in calculations of sala

We have also observed some kind of rearrangement of the nodal structure of the calculated wave function when crossing the point $g/(\hbar\omega\ell)\simeq 3$ at $\alpha=-0.0304552$ (see Figs.19 and 20). However, due to strong interatomic coupling $g$ and considerable anharmonic parameter $\alpha$ we have strong mixing of the states with different quantum numbers here and cannot interpret the effect as a simple transition from one pure quantum state to another one like near the point $g=0$.

The calculated dependence of the initial atomic distribution (the probability density $|\psi(x,y,t=0)|^{2}$) on the coupling constant $g$ clarifies the monotonic increase of the tunneling rate from the lower energy branch of the excited states with the increase of $g$: with the increase of $g$ the maximums of the probability density moves closer to the regions (3) in Fig.9 (see Fig.20). The dependence on $g$ for the tunneling rate from the upper energy branch of the excited states has non-monotonic character due to more complicated dependence of the probability density $|\psi(x,y,t=0)|^{2}$ on $g$: with deviation to the left or right from the point $g=0$ the maximums of the probability density first approach to the regions (3) in Fig.9, but then start to move from the regions (3) to (4) and finally stabilize (see Fig.19).

[TABLE]

Fig.19: Wave function, $\psi(x,y,t=0)$ (in $\ell^{-1}$ units), of the upper energy branch in Fig.17. Coupling strength $g$ is in $\hbar\omega\ell$ units.

[TABLE]

Fig.20: Wave function, $\psi(x,y,t=0)$ (in $\ell^{-1}$ units), of the lower energy branch in Fig.17. Coupling strength $g$ is in $\hbar\omega\ell$ units.

V V. Conclusion

We have investigated the tunneling dynamics of a 1D system of two interacting atoms confined in the anharmonic trap (2.5). We have calculated the tunneling rates $\gamma$ from the three lowest atomic bound states as a function of the coupling constant $g$ for different values $\alpha$ of the trap anharmonicity. It was found that in the tunneling from the upper energy branch of the excited states, $\gamma$ behaves non-monotonically and the sequential particle channel dominates in the tunneling. Note, that the domination of the sequential tunneling was also observed in the box-shaped potential model from the ground state of a rectangular potential wellhunn . When the atoms are initially in the lower energy branch of the excited states, $\gamma$ grows very fast with increasing of the coupling strength $g$ at the beginning of the tunneling. Then, it was found, that the tunneling passes in new regime(s) with more slow tunneling rate(s) due to the competition of the tunneling with the transition to the upper energy branch of the excited states and to the ground state. When the atoms tunnel from the ground state, $\gamma$ grows slowly and monotonically with increasing $g$.

We have also analyzed a rearrangement $(0,2)\rightleftarrows(2,0)$ of the spectrum in the limit $g\rightarrow\pm 0$ of noninteracting atoms with the exchange of the wave function nodal structure. More complicated rearrangement of the nodal structure of the calculated wave function of the confined pair of atoms was found when crossing the point $g/(\hbar\omega\ell)\simeq 3$.

The developed computational scheme can be extended to technically more complicated, but close to current and planning experimentszuern ; zuern2 , a problem about the tunneling from quasi-1D cigar-like and quasi-2D pancake-like traps. Including the spin dynamics into the model for tunneling process is another actual problem which can be investigated with the developed technique.

The authors are very grateful to P. Schmelcher, V. Pupyshev, Yu.V. Popov, and S.I. Vinitsky for very helpful comments and fruitful discussions.

The work was supported by the Grant of the Plenipotentiary Representative of the Republic of Kazakhstan to JINR.