ABSTRACT

We theoretically investigate the quantum states of a Hamiltonian model for quasi-one-dimensional ultracold trapped gases. From the ansatz given by the numerical solution of the Schrödinger equation of the system, we develop a scattering potential functional form and an approximate solution for the analytical approach of the model. We obtain the set of approximate eigenstates and eigenenergies that can be used in future improvements on the study of atomic scattering in low dimensional ultracold gases. We also show that there is a parity inversion of the ground state of the model as the interaction strength increases.

Keywords:

Ultracold Trapped Gases, Confined Eigenstates, Resonance

1. Introduction

The current ability in developing experiments in low dimensional systems has recently renewed the interest for both experimental and theoretical studies of such systems. The techniques of laser cooling of gases confined in optical traps, reaching very low temperatures, allow the physical realization of quasi-1D and quasi-2D systems, with several applications in the study of new collective phenomena (such as Bose-Einstein condensation, BCS superfluidity and the BCS-BEC crossover [1] [2] [3] , exotic (and non-exotic) phases such as Mott, Anderson like, Sarma and mixed phases [4] [5] [6] [7] and Confinement Induced Resonances (CIRs) [8] - [13] ). The interaction between atoms trapped on these systems became a key point on the understanding of the underlying physics involved, with special relevance to the study and modeling of atomic collisions in these low dimensional and low temperature states.

Binary collisions in such confined systems are usually studied in the context of s-wave scattering by means of a very simple model, where the atoms are confined in two of the spatial directions by means of harmonic [8] [14] [15] or infinite quantum well [15] [16] potentials and the interaction between two colliding atoms is modeled by the regularized Huang’s pseudopotential [17] . The strength of the Huang’s potential corresponds to the theoretical representation of the tunable interaction between the atoms by Feshbach resonance [18] [19] [20] in actual experiments. In these approaches, the scattering amplitude is determined by means of the expansion of the scattering wave function [8] [15] [21] into the eigenstates of the transverse unperturbed (Huang’s potential independent) Hamiltonian. As a result, it is possible to relate the effective 1D coupling strength to the s-wave scattering length (an experimentally accessible quantity, related to the Feshbach resonance interaction), showing the existence of phenomena such as CIRs [14] [22] - [29] and multiple CIRs [13] [30] [31] [32] [33] [34] and providing the connection between theoretical and experimental quantities in the study of ultracold gases.

In this paper we will turn our attention to the low energy eigenstates of the atoms in these quasi-1D traps including the scattering potential effects―a nontrivial issue, despite the simplicity of the model we will employ. The aim of such procedure is to provide a set of wave functions that can be used as a basis for the scattering wave function expansion, allowing future refinements on the results previously obtained by employing the eigenstates of the unperturbed Hamiltonian. This set of wave functions constitutes the main result of this paper. Here, we will study a transverse quantum well confined system, defined in section 2, by means of two approaches. First, in section 3, we will numerically study a confined wave function interacting with a Kronecker delta function scattering potential. Then, inspired by the results obtained in the numerical approach, we will also construct, in Section 4, an approximate solution to the exact 1D case subject to a Dirac delta function scattering potential and extend this approach, in Section 5, to define the quasi-1D system for which we will present the referred set of eigenstates. Besides, we will show that, depending on the geometry of the system and as a function of the interaction strength, it is possible the occurrence of the parity inversion of the ground state of the system. Finally, in Section 6, we present our concluding remarks.

2. The Quasi-1D Model for a Scattering Atom Confined in an Optical Trap

As discussed earlier, the quantum states of the system representing a bosonic atom confined in an optical trap can be modeled by a simple quantum model of a particle subject to a harmonic or infinite well transverse potential and to a scattering potential located at the center of the system. Here we will consider the boundary of the system as an infinite well with two of its dimensions considerable lower than the (elongated) other one, $L_z\gg L_x\approx L_y$ , where $L_x,L_y$ and $L_z$ are the dimensions of the potential well. Other geometries can alter the obtained results quantitatively, but the qualitative results for the system should not be dependent on the trap geometry. The scattering potential, representing a target at the center of the system, is modeled by a spherical symmetric potential $V\left(r\right)=g\xi \left(r\right)$ , where g stands for the interaction strength between the scattering atoms. In Olshanii approach [8] , $\xi \left(r\right)$ is the Huang’s regularized pseudopotential [17] ,

$\xi \left(r\right)=\delta \left(r\right)\frac{\partial }{\partial r}$ .

The time independent Schrödinger equation for the system is described thus by

$-{\nabla }^2\Psi \left(r\right)+\left(\frac{2m}{{\hslash }^2}g\right)\xi \left(r\right)\Psi \left(r\right)-\left(\frac{2m}{{\hslash }^2}E\right)\Psi \left(r\right)=0$ (1)

where m is the effective mass of the scattering atoms, $\Psi \left(r\right)=\Psi \left(x,y,z\right)$ is the wave function of the atom into the optical trap, E is the eigenenergy of the

eigenstate described by $\Psi$ and $r=\left|r\right|=\sqrt{x^2+y^2+z^2}$ . Equation (1) is subjected

to the boundary conditions representing the optical trap, given by

$\Psi \left(0,y,z\right)=\Psi \left(x,0,z\right)=\Psi \left(x,y,0\right)=0$ (2)

and

$\Psi \left(L_x,y,z\right)=\Psi \left(x,L_y,z\right)=\Psi \left(x,y,L_z\right)=0$ . (3)

For simplicity, from now on, we will use dimensionless variables, i.e., we will replace $\frac{2m}{{\hslash }^2}g\to g$ and $\frac{2m}{{\hslash }^2}E\to E$ , so that the Schrödinger equation is now described by

$-{\nabla }^2\Psi \left(r\right)+g\xi \left(r\right)\Psi \left(r\right)-E\Psi \left(r\right)=0$ . (4)

Despite the simplicity of the model given by (4), its solution is nontrivial, mainly due to the symmetry of the system―while the boundary conditions have cartesian symmetry (or cylindrical, in the harmonic potential case), the scattering potential has spherical symmetry. The exception is in the particular case $g=0$ , where the analytical solutions of (4), under the boundary conditions given by (2) and (3), constitute a simple quantum mechanics exercise, and are given by

$\Psi \left(x,y,z\right)=\sqrt{\frac{8}{L_x{L}_y{L}_z}}\sin \left(n_x\text{π}\frac{x}{L_x}\right)\sin \left(n_y\text{π}\frac{y}{L_y}\right)\sin \left(n_z\text{π}\frac{z}{L_z}\right)$ , (5)

with eigenenergies given by

$E_{n_x,n_y,n_z}={\text{π}}^2\left(\frac{n_x^2}{L_x^2}+\frac{n_y^2}{L_y^2}+\frac{n_z^2}{L_z^2}\right)$ , (6)

where $n_x,n_y$ and $n_z$ are integers greater than or equal to one. Here, we will refer to the z coordinate dependence of the wave functions as the longitudinal wave function, since we are choosing z as the elongated axis, and to the x and y dependent parts as the transverse wave function. For $L_z\gg L_x\approx L_y$ , the unperturbed solution, given by (5) and (6), shows that for low energies (temperatures), the transverse wave function remains in its ground state, since any transition from $n_{x\left(y\right)}=1$ to $n_{x\left(y\right)}=2$ state implies in a larger amount of energy than the transitions between different $n_z$ levels. So, in the transverse isotropic case ( $L_x=L_y=L$ ) and at low energy ( $n_x=n_y=1$ ) we can rewrite (6) as

$E_{n_z}=\frac{{\text{π}}^2}{L^2}\left(2+n_z^2\frac{L^2}{L_z^2}\right)$ . (7)

The energy difference between the first excited state ( $n_z=2$ ) and the ground state ( $n_z=1$ ) is, thus, given by $\Delta E=E_2-E_1=3{\text{π}}^2/L_z^2$ .

3. Numerical Results on a Discrete Lattice

The solution of (4) can be numerically performed by discretizing the region of interest in a lattice. For the computation of its eigenstates and eigenvalues, we employed the relaxation method [35] in which, from the introduction of a (imaginary) time parameter $\tau$ , we redefine (1) as

$\kappa \frac{\partial \Psi \left(r,\tau \right)}{\partial \tau }=-{\nabla }^2\Psi \left(r,\tau \right)+g\xi \left(r\right)\Psi \left(r,\tau \right)-E\left(\tau \right)\Psi \left(r,\tau \right)$ (8)

where $\kappa$ is the relaxation constant and the energy $E\left(\tau \right)$ , associated to a given wave function configuration $\Psi \left(r,\tau \right)$ , being calculated, from (1), by

$E\left(\tau \right)=\frac{{\displaystyle \int \left(-{\Psi }^{*}{\nabla }^2\Psi +g\xi \left(r\right){\Psi }^{*}\Psi \right)\text{d}r}}{{\displaystyle \int {\Psi }^{*}\Psi \text{d}r}}$ .

Beginning from an initial solution $\Psi \left(r,0\right)$ , and assuming the convergence of

the solution at $\tau \to \infty$ , we get, at this limit, ${\partial \Psi /\partial \tau |}_{\tau \to \infty }=0$ , i.e., $\psi \left(r,\infty \right)$ is the required solution of (4) with eigenenergy $E=E\left(\tau =\infty \right)$ .

The region of interest can be discretized in a lattice of $\left(N_x+1\right)\times \left(N_y+1\right)\times \left(N_z+1\right)$ sites, in which we can define the discretized wave function ${\Psi }_{i,j,k,\tau }$ , where $0\le i\le N_x$ and with similar limits on the y and z directions. For symmetry reasons, we take $N_x$ , $N_y$ and $N_z$ as even integers, and the central site is thus indexed as $i=\frac{N_x}{2}$ , $j=\frac{N_y}{2}$ and $k=\frac{N_z}{2}$ . The Schrödinger equation can be discretized by applying the finite differences [36] expression for the Laplacian of the wave function ${\nabla }^2\Psi$ , given by

$\begin{array}{l}{\nabla }^2\Psi ={\left(\frac{N_x}{L_x}\right)}^2\left({\Psi }_{i+1,j,k}-2{\Psi }_{i,j,k}+{\Psi }_{i-1,j,k}\right)+{\left(\frac{N_y}{L_y}\right)}^2\left({\Psi }_{i,j+1,k}-2{\Psi }_{i,j,k}+{\Psi }_{i,j-1,k}\right)\\ +{\left(\frac{N_z}{L_z}\right)}^2\left({\Psi }_{i,j,k+1}-2{\Psi }_{i,j,k}+{\Psi }_{i,j,k-1}\right)\end{array}$

and also by discretizing the $\tau$ derivative as $\frac{\partial \Psi }{\partial \tau }=\frac{{\Psi }_{i,j,k,\tau +1}-{\Psi }_{i,j,k,\tau }}{\Delta \tau }$ , i.e., by

applying the forward-time central-space method [36] . For the discrete simulations we will assume $\xi \left(r\right)$ concentrated at the scattering center, and implement the scattering potential by means of the Kronecker delta function, i.e.,

$g\xi \left(r\right)\to g{\xi }_{i,j,k}=g{\delta }_{i,j,k}=\{\begin{array}{l}g,\text{for}i=N_x/2,j=N_y/2\text{and}k=N_z/2;\\ 0,\text{elsewhere,}\end{array}$

resulting in the following expression for ${\Psi }_{i,j,k,\tau +1}$ computed from an earlier solution ${\Psi }_{i,j,k,\tau }$ :

${\Psi }_{i,j,k,\tau +1}={\Psi }_{i,j,k,\tau }+\frac{1}{\kappa }\left(-\frac{{\hslash }^2}{2m}{\nabla }^2\Psi +g{\delta }_{i,j,k}{\Psi }_{i,j,k,\tau }-E\left(\tau \right){\Psi }_{i,j,k,\tau }\right)$

The discrete versions of the analytical solutions given by (5) can be taken as the initial configurations for the relaxation process, and the evolution in $\tau$ , Equation (8), can be performed until convergence is achieved. We will use, as the convergence condition, $\left|\left(E\left(\tau +1\right)-E\left(\tau \right)\right)/E\left(\tau +1\right)\right|<\epsilon$ , where $\epsilon$ is a given convergence radius.

We computed the numerical solution of (8) in a lattice of $51\times 51\times 1001$ sites,

with isotropic lattice spacing $\frac{L_x}{N_x}=\frac{L_y}{N_y}=\frac{L_z}{N_z}=h$ . At this condition, the energies

(and also the interaction strength g) scales with $h^{-2}$ , as one can see by replacing (7) in (8), so the eigenenergies and interaction strength can be expressed by the scale invariant quantities $h^{2}E$ and $h^{2}g$ , respectively. However, the numerical errors involved on the finite differences Laplacian are of order $h^2$ , i.e., the errors are smaller as h is made smaller. For $\kappa =1$ , Equation (8) converges with radius $\epsilon \le {10}^{-7}$ for any value of g used after 1000 steps ( $\tau =1000$ ). The computed eigenenergy for the unperturbed case ( $g=0$ ) differs from the exact analytical value, given by (6), for $n_x=n_y=n_z=1$ , by only 0.033%.

In figure 1 we plot the transverse wave function at $z=L_z/2$ for $g{h}^2=5$ . The resulting wave function shows a reasonable fit to the unperturbed wave function in the region far from the scattering center location, and an inverse peak (the drop in the top point of the surface) at the center of the box, decreasing the probability to locate the atom at this point, as expected for a repulsive interaction. For other z values far from $z=L_z/2$ the transverse wave function is only slightly modified from the unperturbed solution (5).

In figure 2 we show the transverse wave function for different values of g at fixed $y=L/2$ and $z=L_z/2$ . As we can see, increasing g results in the deepening of the inverse peak and in the enlargement of the region perturbed by the target. The eigenenergy dependence with g is displayed in figure 3. The eigenenergy of the ground state has a non-linear dependence with the coupling strength, increasing from its unperturbed value $h^2{E}_1=7.902956\times {10}^{-3}$ up to the finite energy limit $h^2{E}_1=7.916836\times {10}^{-3}$ in the infinity coupling limit. The first excited state, however, has odd parity relative to the central point, i.e.

$\Psi \left(x,y,L_z/2-\epsilon \right)=-\Psi \left(x,y,L_z/2+\epsilon \right)$ . In consequence, $\Psi \left(x,y,L_z/2\right)=0$ , and thus it decouples from the scattering potential. The wave function and the energy of the first excited state are thus independent of g and equal to $h^2{E}_2=7.932565\times {10}^{-3}$ . One can see clearly, thus, that as g increases the energy gap between the ground and the first excited states reduces. The same effect occurs to any pair of other parity even and subsequent odd states.

An interesting phenomenon occurs for large enough $L_z$ : since the energy of the ground state, $E_1$ , increases as g increases, it is possible that this energy reaches the first excited state value (that is independent of g due to its odd parity, as discussed before) for some finite value of g. As we can see from (7), for $g=0$ , $\Delta E$ , the energy gap between the ground and the first excited state, is equal to ${\text{π}}^2\left(2n_z+1\right)/L_z^2$ . Thus, for large enough $L_z$ (small $\Delta E$ ) one should expect the above described resonance condition ( $E_1=E_2$ ) to be reached. For higher values of g, the eigenenergy of the ground state can become higher than the first odd state energy, i.e., there is an inversion of the ground state from an even parity state to an odd parity one. This occurrence is shown in figure 4, for $N_x=N_y=26$ and $N_z=1000$ , where the transition of the ground state from the even to the odd parity state occurs at $h^{2}g\approx 5.7$ . We will study this phenomenon within the analytical approach in Section 5.

The presence of the sharp peak at the scattering center, with discontinuous first derivative and short range, suggests that the $g=0$ wave function is perturbed by the target as an exponential decay dependent on the distance to the scattering center. In comparison with some possible functional forms of wave

functions with these features, we found, for the x dependent term, for example, that

${\Psi }_x\left(x\right)=N\left(1-B{\text{e}}^{-\gamma \left|x-x_c\right|}\right)\sin \left(n_x\text{π}\frac{x}{L_x}\right)$ (9)

fits the numerical computed wave function with a good agreement, as shown in figure 5. Here, $x_c=L_x/2$ and $B$ and $\gamma$ correspond to the fitting parameters of the wave function. By defining ${\Psi }_{x0}=N\sin \left(n_x\text{π}x/L_x\right)$ , we obtain $\ln \left(1-{\Psi }_x/{\Psi }_{x0}\right)=\ln \left(B\right)-\gamma \left|x-x_c\right|$ , and the linear fit of $\ln \left(1-{\Psi }_x/{\Psi }_{x0}\right)$ as a function of $x/L_x$ , depicted on figure 5 for $h^{2}g=5$ , results in $\ln \left(B\right)=40.499$ and $\gamma =81.061$ , with a Pearson correlation coefficient $r^2=0.99788$ .

As suggested by the numerical results, Equation (9) will be used now as an ansatz to the solution to the one dimensional version of the Schrödinger equation of the model.

4. Analytical Approximation―1D Solution

Let us consider the (non normalized) ansatz given by (9),

$\Psi \left(x\right)=\left(1-B{\text{e}}^{-\gamma \left|x-x_C\right|}\right)\sin \left(kx\right)$ (10)

where $x_C=L_x/2$ and $k=n_x\text{π}/L_x$ . For convenience, let us displace, from now on, the origin of the system to the scattering center ( $x-x_C\to x$ ). Thus, the (ground state) wave function and the Schrödinger equation of the system are written, respectively, as

$\Psi \left(x\right)=\left(1-B{\text{e}}^{-\gamma \left|x\right|}\right)\cos \left(kx\right)$ (11)

and

$-{\Psi }^{″}\left(x\right)+g\delta \left(x\right)\Psi \left(x\right)=E\Psi \left(x\right)$ , (12)

where ${\Psi }^{″}\left(x\right)=\frac{{\text{d}}^2\Psi \left(x\right)}{\text{d}{x}^2}$ stands for the Laplacian of $\Psi \left(x\right)$ in 1D and g and E

are given in units of ( $2m/{\hslash }^2$ ), as before. We are also assuming now that the scattering potential is given by $\xi \left(x\right)=\delta \left(x\right)$ , the Dirac delta function. Computing the second derivative of (11) we obtain

$\begin{array}{l}{\Psi }^{″}\left(x\right)=\left(2\gamma B\delta \left(x\right)-{\gamma }^{2}B\mathrm{sgn}{\left(x\right)}^2\right){\text{e}}^{-\gamma \left|x\right|}\cos \left(kx\right)\\ -2k\gamma B\mathrm{sgn}\left(x\right){\text{e}}^{-\gamma \left|x\right|}\sin \left(kx\right)-k^2\left(1-B{\text{e}}^{-\gamma \left|x\right|}\right)\cos \left(kx\right)\end{array}$ (13)

where we have used the following properties for the absolute value of x:

$\frac{\text{d}\left|x\right|}{\text{d}x}=\mathrm{sgn}\left(x\right)=\{\begin{array}{l}-1\text{for}x<0\\ 1\text{for}x>0\end{array}=2\Theta \left(x\right)-1$ , $\frac{\text{d}\mathrm{sgn}\left(x\right)}{\text{d}x}=2\frac{\text{d}\Theta \left(x\right)}{\text{d}x}=2\delta \left(x\right)$ and

$\mathrm{sgn}{\left(x\right)}^2=1$ , with $\Theta \left(x\right)$ standing for the Heaviside function. Replacing these results on (12) we get

$\begin{array}{l}-\left(2\gamma \delta \left(x\right)-{\gamma }^2+k^2\right)B{\text{e}}^{-\gamma \left|x\right|}\cos \left(kx\right)+k^2\cos \left(kx\right)+2k\gamma B\mathrm{sgn}\left(x\right){\text{e}}^{-\gamma \left|x\right|}\sin \left(kx\right)\\ +g\delta \left(x\right)\left(1-B{\text{e}}^{-\gamma \left|x\right|}\right)\cos \left(kx\right)=E\left(1-B{\text{e}}^{-\gamma \left|x\right|}\right)\cos \left(kx\right).\end{array}$ (14)

It is worth to mention that (14) can only be satisfied due to the presence of the absolute value function $\left|x\right|$ on the ansatz solution. The first delta function appearing on the right hand side of (14) comes from the second derivative of $\left|x\right|$ , and this term is necessary to cancels out the other delta function dependent terms that come from the scattering potential $g\delta \left(x\right)$ . By imposing the delta function dependent terms to vanish and considering that these terms will be non null only at $x=0$ (where ${\text{e}}^{-\gamma \left|x\right|}=1$ and $\cos \left(kx\right)=1$ ), we found

$B=\frac{g}{g+2\gamma }$ . (15)

The remaining terms can be expanded in orders of x and $\left|x\right|$ . At zero order, we have

$E=k^2+\frac{B}{1-B}{\gamma }^2=k^2+\frac{g}{2}\gamma$ , (16)

where we have applied (15) in the last equality of (16). Collecting now the terms at first order in $\left|x\right|$ on (14) and applying the identity $x\mathrm{sgn}\left(x\right)=\left|x\right|$ where necessary, we obtain

$E=3k^2-{\gamma }^2$ . (17)

Equations (15)-(17) constitute a set of three equations with three unknown variables, E, $\gamma$ and B. By replacing (17) in (16) we solve the equation for $\gamma$ , obtaining

$\gamma =\frac{\sqrt{g^2+32k^2}-g}{4}$ . (18)

Equation (18) can finally be replaced on (15) and (16) to give us the values of B and E respectively.

The other relations obtained from the expansion of (14) in higher orders are not satisfied by the ansatz (11). These terms are of orders equal or higher than 2, and since the region perturbed by the scattering center is concentrated at the vicinity of $x=0$ , as we saw on the numerical calculations, the errors are of order $x^2/L_x^2<\left|x\right|/L_x$ , i.e., smaller than the terms considered on the solution (15)-(18).

The analytical approximation of the solution for the 1D version of our model will be now used as a guideline for the analytical solution of the quasi-1D case.

5. Analytical Approximation―Quasi-1D Solution

From the 1D solution one can clearly observe that one of the key steps of the solution is to collect the singular part that arises from the Laplacian of ${\text{e}}^{-\gamma \left|x\right|}$ together with the singular scattering potential, resulting in the determination of B, given by (15). The natural extension of the ansatz solution given by (11), taking into account the spherical symmetry of the scattering potential, is to replace the

${\text{e}}^{-\gamma \left|x\right|}$ decay by ${\text{e}}^{-\gamma r}={\text{e}}^{-\gamma \sqrt{x^2+y^2+z^2}}$ . Let us, thus, consider now the analytical solution

for the 3D wave functions confined by the transverse quantum well potential employing the following ansatz for the non normalized wave functions in the transverse isotropic ( $L_x=L_y=L$ ) case:

$\Psi \left(x,y,z\right)=\left(1-B{\text{e}}^{-\gamma r}\right)\cos \left(k_{x}x\right)\cos \left(k_{y}y\right)\cos \left(k_{z}z\right)$ , (19)

where $k_x=\frac{n_x\text{π}}{L_x}$ , $k_y=\frac{n_y\text{π}}{L_y}$ and $k_z=\frac{n_z\text{π}}{L_z}$ . The first spatial derivatives of $\Psi \left(x,y,z\right)$ are given by

$\frac{\partial \Psi }{\partial x}=\left(\gamma Bsgnx\left(x,y,z\right){\text{e}}^{-\gamma r}\cos \left(k_{x}x\right)-\left(1-B{\text{e}}^{-\gamma r}\right)k_x\sin \left(k_{x}x\right)\right)\cos \left(k_{y}y\right)\cos \left(k_{z}z\right)$ (20)

with

$sgnx\left(x,y,z\right)=\frac{x}{\sqrt{x^2+y^2+z^2}}$ (21)

and with similar expressions for $\frac{\partial \Psi }{\partial y}$ and $\frac{\partial \Psi }{\partial z}$ , also with the correspondent

definitions of the sgny and sgnz functions. The functions sgnx, sgny and signz coincide with the usual sgn function only when computed over the axes, remaining undefined at the origin. Now defining $f=B^{-\gamma r}$ we obtain, for the Laplacian of (19)

$\begin{array}{l}{\nabla }^2\Psi \left(x,y,z\right)\\ =\left(2f_x{k}_x\sin \left(k_{x}x\right)-\left(f_{xx}+\left(1-f\right)k_x^2\right)\cos \left(k_{x}x\right)\right)\cos \left(k_{y}y\right)\cos \left(k_{z}z\right)\\ +\left(2f_y{k}_y\sin \left(k_{y}y\right)-\left(f_{yy}+\left(1-f\right)k_y^2\right)\cos \left(k_{y}y\right)\right)\cos \left(k_{x}x\right)\cos \left(k_{z}z\right)\\ +\left(2f_z{k}_z\sin \left(k_{z}z\right)-\left(f_{zz}+\left(1-f\right)k_z^2\right)\cos \left(k_{z}z\right)\right)\cos \left(k_{x}x\right)\cos \left(k_{y}y\right)\end{array}$ (22)

with

$f_x=-\gamma sgnx\left(x,y,z\right)f,f_y=-\gamma sgny\left(x,y,z\right)f,f_z=-\gamma sgnz\left(x,y,z\right)f$ , (23)

$\begin{array}{l}{f}_{xx}=\left({\gamma }^{2}sgn{x}^2-2\gamma {\delta }_x\right)f,f_{yy}=\left({\gamma }^{2}sgn{y}^2-2\gamma {\delta }_y\right)f\\ \text{and}{f}_{zz}=\left({\gamma }^{2}sgn{z}^2-2\gamma {\delta }_z\right)f\end{array}$ (24)

The symbols ${\delta }_x$ , ${\delta }_y$ and ${\delta }_z$ are defined as

$\begin{array}{l}{\delta }_x=\frac{1}{2}\frac{\partial }{\partial x}sgnx\left(x,y,z\right),{\delta }_y=\frac{1}{2}\frac{\partial }{\partial y}sgny\left(x,y,z\right)\\ \text{and}{\delta }_z=\frac{1}{2}\frac{\partial }{\partial z}sgnz\left(x,y,z\right)\end{array}$ (25)

Care must be taken on calculating the results of (25). If one calculates ${\delta }_x$ , for example, directly from (21) and (25), obtains, except at the origin,

${\delta }_x=\frac{1}{2}\frac{\partial }{\partial x}sgnx\left(x,y,z\right)=\frac{1}{\sqrt{x^2+y^2+z^2}}-\frac{x^2}{{\left(x^2+y^2+z^2\right)}^{3/2}}$ , (26)

and with similar results for ${\delta }_y$ and ${\delta }_z$ , resulting in

${\delta }_x+{\delta }_y+{\delta }_z=\frac{1}{\sqrt{x^2+y^2+z^2}}=\frac{1}{r},r\ne 0.$ (27)

However, at the origin the ${\delta }_{x,y,z}$ functions contain singularities, as we saw explicitly in the 1D case ( $\text{d}\mathrm{sgn}\left(x\right)/\text{d}x=2\delta \left(x\right)$ ) that have to be taken into account. From (21) and (25) we can see that

${\delta }_x+{\delta }_y+{\delta }_z=\frac{1}{2}\nabla \cdot \left(\frac{r}{r}\right)$ (28)

that contains a regular Coulomb like part (1/r) plus a singular part at the origin.

Replacing the Laplacian of $\Psi$ in the Schrödinger Equation (4), we obtain

$\begin{array}{l}(\left({\gamma }^{2}sgn{x}^2-2\gamma {\delta }_x\right)f\cos \left(k_{x}x\right)+2\gamma \mathrm{sgnx}f{k}_x\sin \left(k_{x}x\right)\\ +\left(1-f\right)k_x^2\cos \left(k_{x}x\right))\cos \left(k_{y}y\right)\cos \left(k_{z}z\right)\\ +(\left({\gamma }^{2}sgn{y}^2-2\gamma {\delta }_y\right)f\cos \left(k_{y}y\right)+2\gamma \mathrm{sgny}f{k}_y\sin \left(k_{y}y\right)\\ +\left(1-f\right)k_y^2\cos \left(k_{y}y\right))\cos \left(k_{x}x\right)\cos \left(k_{z}z\right)\end{array}$

$\begin{array}{l}+(\left({\gamma }^{2}sgn{z}^2-2\gamma {\delta }_z\right)f\cos \left(k_{z}z\right)+2\gamma \mathrm{sgnz}f{k}_z\sin \left(k_{z}z\right)\\ +\left(1-f\right)k_z^2\cos \left(k_{z}z\right))\cos \left(k_{x}x\right)\cos \left(k_{y}y\right)\\ +g\xi \left(r\right)\left(1-f\right)\cos \left(k_{x}x\right)\cos \left(k_{y}y\right)\cos \left(k_{z}z\right)\\ =E\left(1-f\right)\cos \left(k_{x}x\right)\cos \left(k_{y}y\right)\cos \left(k_{z}z\right).\end{array}$ (29)

Now, following the steps of the 1D solution presented in section 4, by recognizing that the terms of (29) that contain singularities at the origin are the terms depending on the functions $\xi \left(r\right)$ , ${\delta }_x$ , ${\delta }_y$ and ${\delta }_z$ , we collect such terms, resulting in

$-2\gamma \left({\delta }_x+{\delta }_y+{\delta }_z\right)f+g\xi \left(r\right)\left(1-f\right)=0$ . (30)

Here, the role of (27) is twofold―first, it will be used to define the scattering potential $\xi \left(r\right)$ as

$\xi \left(r\right)={\delta }_x+{\delta }_y+{\delta }_z=\frac{1}{2}\nabla \cdot \left(\frac{r}{r}\right)$ , (31)

i.e., we will consider here the case where the approximate solution can be found by following the 1D solution procedure. Equation (31) defines, thus, the quasi-1D model we are studying here. Also, as we shall see, this solution reproduces qualitatively well the numerical results previously obtained for the Kronecker delta scattering potential. The scattering potential given by (31) contains the singular and spherical symmetric nature at the origin of the original Dirac delta function [8] . Besides, it contains also a regular Coulomb like term, reproducing also an expected feature of the atomic scattering. Second, by replacing $f=B^{-\gamma r}$ and (31) on (30), we obtain, at the origin

$B=\frac{g}{g+2\gamma }$ , (32)

the same result obtained in the 1D system, Equation (15). Collecting next the remaining terms of (29) in $x=y=z=0$ , we obtain, again as in (16),

$E=k^2+\frac{B}{1-B}{\gamma }^2=k^2+\frac{g}{2}\gamma$ . (33)

Finally, by collecting now the terms on the first order on the spatial variables, we obtain

$Er=k^{2}r-{\gamma }^{2}r+2\frac{\left(k_x^2{x}^2+k_y^2{y}^2+k_z^2{z}^2\right)}{\left(x^2+y^2+z^2\right)}r$ . (34)

The third term on the right hand side of (34) is a mathematical indeterminacy. The limit of this term at the origin is dependent on the path we choose to approach the origin. Observe that in the one dimensional case this indeterminacy is not present―by eliminating the terms in y and z in (34) we retrieve (17). To fix this indeterminacy we will assume the symmetric path ( $x=y=z\to 0$ ), and compare the result obtained by using this assumption with the one obtained by the numerical simulation. Thus, we have

${\frac{\left(k_x^2{x}^2+k_y^2{y}^2+k_z^2{z}^2\right)}{\left(x^2+y^2+z^2\right)}|}_{x=y=z}=\frac{k^2}{3}$ . (35)

Two physical arguments can also support the specific choice of the path employed to solve the indeterminacy―1) in the case $k_x=k_y=k_z$ there is no indeterminacy for any set of coordinates (x, y, z), and the result is the same as in (35); 2) the result of (35) ensures the wave function energy E to be dependent of the modulus of the wave-number vector $k$ , and not of its individual components, $k_x$ , $k_y$ and $k_z$ , in an independent way, breaking the symmetry of the result relative to the interchanging of the x, y or z dimensions. Using this assumption, we obtain, for the energy of the ansatz wave function,

$E=\frac{5}{3}{k}^2-{\gamma }^2$ , (36)

and, replacing (36) on (33),

$\gamma =\frac{\sqrt{g^2+\frac{32}{3}{k}^2}-g}{4}$ , (37)

which, replaced in (32) and (33), gives us the expressions for B and E in the quasi-1D approach.

The agreement of the approximate solution given by (19), (32), (33) and (37), with the unknown exact one can be evaluated by computing the Schrödinger equation residue, defined by

$R\left(r\right)=-{\nabla }^2\Psi \left(r\right)+g\xi \left(r\right)\Psi \left(r\right)-E\Psi \left(r\right)$ . (38)

For the exact solution we have $R\left(r\right)=0$ , so we expect, for the approximate solution, $R\left(r\right)$ to be very small.

Figure 6(a) presents the eigenfunction $\Psi$ computed at $x=y=0$ , as a function of z with $g=0.2$ and $L_z=40L$ (solid line). As we can see, the result shows the same qualitative behavior as the previously obtained in numerical calculations with the Kronecker delta scattering potential. The dashed line in figure 6(a) (highlighted in figure 6(b)) shows the small residue of the approximate solution.

In figure 7 we show the dependence of the parameter B (solid line) and the ground state eigenenergy (dashed line) as a function of g for $L_z=40L$ and $L=25$ (arbitrary units). Again, we obtain the same qualitative behavior in comparison with the numerical calculations. In the infinite coupling limit, $g\to \infty$ , we obtain $B=1$ and $E=5k^2/3$ .

The analytical approximation allows us to determine the conditions for when the parity inversion can be observed. As we already had determined, for $g=0$ the ground state energy is given by (7), i.e., $E_1=k^2$ , as we can see also from (33). On the other limit, $g\to \infty$ , an expansion of (37) shows that $\gamma \to 0$ , and from

(36) we get $E_1=\frac{5}{3}{k}^2$ at this limit, so that $\Delta E_1=E_{1\left(g\to \infty \right)}-E_{1\left(g=0\right)}=\frac{2}{3}{k}^2$ . On other side, the energy difference between the first excited state and the ground state is given by $\Delta E=3\frac{{\text{π}}^2}{L_z^2}$ , as we have seen. So, the condition for which the

parity inversion (when the odd parity state energy becomes lower than the energy of the even parity state) within the approximation employed here is given by

$\frac{2}{3}{k}^2>3\frac{{\text{π}}^2}{L_z^2}$ (39)

or, solving (39) for $L_z$ with $n_x=n_y=n_z=1$ ,

$L_z>\frac{\sqrt{7}}{2}L$ . (40).

As we can see, the condition for occurrence of the parity inversion at some strength g is purely geometric. The differences of the energies between the first even to the first odd parity states as a function of the interaction strength g in two different conditions ((a) $L_z/L=1.2$ and (b) $L_z/L=2$ ) are depicted in figure 8, illustrating again the occurrence of the parity inversion. Besides, as in the quasi-1D system we have $L_z\gg L$ , we can conclude that the condition (40) is always guaranteed in this case.

In conclusion, equation (19) with $\gamma$ and B given by (37) and (32), respectively,

constitute the basis of approximated eigenstates, with eigenenergies given by (33), for the ultracold atomic scattering confined in a optical trap model.

6. Conclusion

In summary, we present both a numerical and an analytical analysis of the eigenstates and eigenenergies of a simple model for atomic scattering of ultracold atoms confined in quasi-1D optical traps. We show that the eigenstates wave functions are modified from the unperturbed solution only in a region near the scattering center, the wave function profile presenting an inverse sharp peak at the vicinity of the scattering center with an exponential decay to the unperturbed solution as we move away from this point. We used the functional form suggested by the numerical calculations as an ansatz that allowed both the definition of a quasi-1D model and the obtaining of its approximate solutions, providing a set of wave functions that can be used as a basis for the scattering wave function expansion. The definition of a quasi-1D model for scattering atoms in ultracold trapped gases and the obtaining of the set of approximate solutions that include the effects of the scattering potential and can be used as a basis for other applications are the main results of the present contribution. We also show that, as the interaction strength increases, a parity inversion of the ground state can occur, depending on the geometry of the system.

Acknowledgements

The authors would like to thank financial support from CNPq, FAPEMIG and CAPES (Brazilian funding agencies).

Cite this paper

e Miranda, P.G., de Freitas, D.N. and Mota, A.L. (2018) Quantum States of Interacting Atoms in Quasi-One-Dimensional Ultracold Trapped Gases. Journal of Modern Physics, 9, 443-460. https://doi.org/10.4236/jmp.2018.93031

References