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    , exotic (and non-exotic) phases such as Mott, Anderson like, Sarma and mixed phases     and Confinement Induced Resonances (CIRs)  -  ). 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    or infinite quantum well   potentials and the interaction between two colliding atoms is modeled by the regularized Huang’s pseudopotential  . The strength of the Huang’s potential corresponds to the theoretical representation of the tunable interaction between the atoms by Feshbach resonance    in actual experiments. In these approaches, the scattering amplitude is determined by means of the expansion of the scattering wave function    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   -  and multiple CIRs       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, , where and 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 , where g stands for the interaction strength between the scattering atoms. In Olshanii approach  , is the Huang’s regularized pseudopotential  ,
The time independent Schrödinger equation for the system is described thus by
where m is the effective mass of the scattering atoms, is the wave function of the atom into the optical trap, E is the eigenenergy of the
eigenstate described by and . Equation (1) is subjected
to the boundary conditions representing the optical trap, given by
For simplicity, from now on, we will use dimensionless variables, i.e., we will replace and , so that the Schrödinger equation is now described by
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 , 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
with eigenenergies given by
where and 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 , 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 to state implies in a larger amount of energy than the transitions between different levels. So, in the transverse isotropic case ( ) and at low energy ( ) we can rewrite (6) as
The energy difference between the first excited state ( ) and the ground state ( ) is, thus, given by .
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  in which, from the introduction of a (imaginary) time parameter , we redefine (1) as
where is the relaxation constant and the energy , associated to a given wave function configuration , being calculated, from (1), by
Beginning from an initial solution , and assuming the convergence of
the solution at , we get, at this limit, , i.e., is the required solution of (4) with eigenenergy .
The region of interest can be discretized in a lattice of sites, in which we can define the discretized wave function , where and with similar limits on the y and z directions. For symmetry reasons, we take , and as even integers, and the central site is thus indexed as , and . The Schrödinger equation can be discretized by applying the finite differences  expression for the Laplacian of the wave function , given by
and also by discretizing the derivative as , i.e., by
applying the forward-time central-space method  . For the discrete simulations we will assume concentrated at the scattering center, and implement the scattering potential by means of the Kronecker delta function, i.e.,
resulting in the following expression for computed from an earlier solution :
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 , Equation (8), can be performed until convergence is achieved. We will use, as the convergence condition, , where is a given convergence radius.
We computed the numerical solution of (8) in a lattice of sites,
with isotropic lattice spacing . At this condition, the energies
(and also the interaction strength g) scales with , as one can see by replacing (7) in (8), so the eigenenergies and interaction strength can be expressed by the scale invariant quantities and , respectively. However, the numerical errors involved on the finite differences Laplacian are of order , i.e., the errors are smaller as h is made smaller. For , Equation (8) converges with radius for any value of g used after 1000 steps ( ). The computed eigenenergy for the unperturbed case ( ) differs from the exact analytical value, given by (6), for , by only 0.033%.
In figure 1 we plot the transverse wave function at for . 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 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 and . 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 up to the finite energy limit in the infinity coupling limit. The first excited state, however, has odd parity relative to the central point, i.e.
Figure 1. The transverse wave function at z = Lz/2. The wave function presents an inverse peak at the scattering center. Far from this point, the wave function is only slightly changed from the unperturbed solution.
Figure 2. The transverse wave function at y = Ly/2 and z = Lz/2 for different couplings. As g increases, there is a deepening in the inverse peak in the wave function.
Figure 3. The eigenenergies of the even parity ground state (solid line) and of the odd parity first excited state (dashed line) as a function of the scattering strength g for Nz = 19.6Nx. As g increases, the gap between the two states reduces up to a finite difference in the infinite coupling limit.
. In consequence, , 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 . 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 : since the energy of the ground state, , 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 , , the energy gap between the ground and the first excited state, is equal to . Thus, for large enough (small ) one should expect the above described resonance condition ( ) 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 and , where the transition of the ground state from the even to the odd parity state occurs at . 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 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
Figure 4. The eigenenergies of the even parity ground state (solid line) and of the odd parity first excited state (dashed line) as a function of the scattering strength g for Nz = 38.5Nx. At lowers values of g the ground state has even parity symmetry related to the reflection z↔ -z. At a finite g, (h2g » 5.7 for Nx = 26) the ground state and first excited state become degenerate. For larger values of h2g, there is a parity inversion of the ground state―the odd parity state, uncoupled from the scattering potential, becomes the ground state of the system.
functions with these features, we found, for the x dependent term, for example, that
fits the numerical computed wave function with a good agreement, as shown in figure 5. Here, and and correspond to the fitting parameters of the wave function. By defining , we obtain , and the linear fit of as a function of , depicted on figure 5 for , results in and , with a Pearson correlation coefficient .
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),
Figure 5. Linear fit of the numerical computed wave function to the exponential decay in (9), with Pearson correlation coefficient given by r2 = 0.99788, presenting a good agreement of the numerical results to the functional form (9).
where and . For convenience, let us displace, from now on, the origin of the system to the scattering center ( ). Thus, the (ground state) wave function and the Schrödinger equation of the system are written, respectively, as
where stands for the Laplacian of in 1D and g and E
are given in units of ( ), as before. We are also assuming now that the scattering potential is given by , the Dirac delta function. Computing the second derivative of (11) we obtain
where we have used the following properties for the absolute value of x:
, with standing for the Heaviside function. Replacing these results on (12) we get
It is worth to mention that (14) can only be satisfied due to the presence of the absolute value function on the ansatz solution. The first delta function appearing on the right hand side of (14) comes from the second derivative of , and this term is necessary to cancels out the other delta function dependent terms that come from the scattering potential . By imposing the delta function dependent terms to vanish and considering that these terms will be non null only at (where and ), we found
The remaining terms can be expanded in orders of x and . At zero order, we have
where we have applied (15) in the last equality of (16). Collecting now the terms at first order in on (14) and applying the identity where necessary, we obtain
Equations (15)-(17) constitute a set of three equations with three unknown variables, E, and B. By replacing (17) in (16) we solve the equation for , obtaining
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 , as we saw on the numerical calculations, the errors are of order , 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 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
decay by . 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 ( ) case:
where , and . The first spatial derivatives of are given by
and with similar expressions for and , 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 we obtain, for the Laplacian of (19)
The symbols , and are defined as
Care must be taken on calculating the results of (25). If one calculates , for example, directly from (21) and (25), obtains, except at the origin,
and with similar results for and , resulting in
However, at the origin the functions contain singularities, as we saw explicitly in the 1D case ( ) that have to be taken into account. From (21) and (25) we can see that
that contains a regular Coulomb like part (1/r) plus a singular part at the origin.
Replacing the Laplacian of in the Schrödinger Equation (4), we obtain
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 , , and , we collect such terms, resulting in
Here, the role of (27) is twofold―first, it will be used to define the scattering potential as
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  . Besides, it contains also a regular Coulomb like term, reproducing also an expected feature of the atomic scattering. Second, by replacing and (31) on (30), we obtain, at the origin
the same result obtained in the 1D system, Equation (15). Collecting next the remaining terms of (29) in , we obtain, again as in (16),
Finally, by collecting now the terms on the first order on the spatial variables, we obtain
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 ( ), and compare the result obtained by using this assumption with the one obtained by the numerical simulation. Thus, we have
Two physical arguments can also support the specific choice of the path employed to solve the indeterminacy―1) in the case 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 , and not of its individual components, , and , 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,
and, replacing (36) on (33),
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
For the exact solution we have , so we expect, for the approximate solution, to be very small.
Figure 6(a) presents the eigenfunction computed at , as a function of z with and (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 and (arbitrary units). Again, we obtain the same qualitative behavior in comparison with the numerical calculations. In the infinite coupling limit, , we obtain and .
Figure 6. (a) The approximate solution of the Schrödinger equation in the quasi-1D system, with scattering potential given by (31) (solid line). The wave function has the same qualitative behavior of the numerical solution of the Kronecker delta potential. The dashed line represents the residue of the Schrödinger equation (amplified on (b)).
Figure 7. The coefficient B and the ground state eigenenergy as a function of the scattering strength g in the quasi-1D approximate solution. As g ® ¥, B ® 1 and E ® 5k2/3.
The analytical approximation allows us to determine the conditions for when the parity inversion can be observed. As we already had determined, for the ground state energy is given by (7), i.e., , as we can see also from (33). On the other limit, , an expansion of (37) shows that , and from
(36) we get at this limit, so that . On other side, the energy difference between the first excited state and the ground state is given by , 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
or, solving (39) for with ,
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) and (b) ) are depicted in figure 8, illustrating again the occurrence of the parity inversion. Besides, as in the quasi-1D system we have , we can conclude that the condition (40) is always guaranteed in this case.
In conclusion, equation (19) with and B given by (37) and (32), respectively,
Figure 8. The energy difference between the even parity to the odd parity state for Lz = 1.2L (solid line) and for Lz = 1.75L (dashed line). For Lz = 1.2L the even parity state remains as the ground state for any value of the scattering strength g. For Lz = 1.75L there is a parity inversion of the ground state from the even to the odd parity one at g » 0.32.
constitute the basis of approximated eigenstates, with eigenenergies given by (33), for the ultracold atomic scattering confined in a optical trap model.
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.
The authors would like to thank financial support from CNPq, FAPEMIG and CAPES (Brazilian funding agencies).