Mathematical Modeling of Ultracold Few-Body Processes in Atomic Traps

. We discuss computational aspects of the developed mathematical models for ultracold few-body processes in atomic traps. The key element of the elaborated computational schemes is a nondirect product discrete variable representation (npDVR) we have suggested and applied to the time-dependent and stationary Schrödinger equations with a few spatial variables. It turned out that this approach is very e ﬃ cient in quantitative analysis of low-dimensional ultracold few-body systems arising in conﬁned geometry of atomic traps. The e ﬃ ciency of the method is demonstrated here on two examples. A brief review is also given of novel results obtained recently.


Introduction
The impressive progress in the physics of the ultracold quantum gases has demanded the elaboration of quantitative models for resonant and multichannel processes occurring in confined geometries of atomic traps. Actually, the conventional few-body theory is no longer valid here and the development of a low-dimensional theory, including the influence of the confinement is needed. In our works we have developed computational methods [1][2][3] for pair collisions in tight atomic waveguides and have found several novel effects in its applications: the confinement-induced resonances (CIRs) in multimode regimes including effects of transverse excitations and deexcitations [2], the so-called dual CIR yielding a complete suppression of the quantum scattering [1], and resonant molecule formation with transferring energy relies to center-of-mass excitation while forming molecules [4]. The last effect was recently confirmed in a Heidelberg experiment [5]. Our calculations have also been used for planning and interpreting the Innsbruck experiment on the investigation of CIRs in ultracold Cs gas [6]. Recently, we have calculated the Feshbach resonance shifts and widths induced by atomic waveguides [7,8] and predicted dipolar CIRs [9].
Here we discuss computational aspects of the developed schemes for ultracold few-body processes in atomic traps. In section 2 we present the description of a nondirect product discrete variable representation (npDVR) by using as an example the time-dependent Schrödinger equation with three spatial variables. This problem arises in the quantitative description of binary collisions of ultracold atoms in anisotropic waveguides. The application of the npDVR for constructing the efficient splitting-up a e-mail: melezhik@theor.jinr.ru DOI: 10

Nondirect Product Discrete-Variable Representation
In our papers [10,11] an alternative approach to the conventional partial-wave analysis was suggested for the Schrödinger equation with a few spatial variables. The angular basis was constructed from eigenfunctions of the angular momentum square defined on the angular grid in the spirit of the discrete variable representation (DVR) [12,13] or the Lagrange-mesh method [14,15] (which is in some sense a subset of DVR). The construction of a DVR orthogonal on the angular grid for one angular variable is a solvable problem if the choice of the grid points coincides with the nodes of the Gauss quadrature [15]. Different kinds of one-dimensional DVRs or Lagrange-meshes are broadly applied for quantum computations [13,15] due to the simplicity (DVR leads to a diagonal form of the interaction and to a compact form of the kinetic energy operator) and efficiency (fast convergence and stability) of this approach. However, an extension of this representation to the two-dimensional case (two angles θ and φ on the unit sphere) is not a trivial problem. Actually, the simple idea to construct the two-dimensional DVR as a direct product of two one-dimensional DVRs leads to significant complication of the matrix of the angular part of the kinetic energy operator [13,16]. As a result, the advantages of the one-dimensional DVR, simplicity and efficiency, are lost [13]. An alternative way to construct twodimensional DVR on an unit sphere is to use the spherical harmonics defined on a two-dimensional grid. However, in this case it becomes not possible to satisfy the orthogonality conditions for all the elements of the fixed set of this basis on the chosen grid [16,17]. To overcome this difficulty, we have suggested [18,19] to use the basis of the orthogonalized combinations of the spherical harmonics on the two-dimensional grid over θ and φ variables. It turned out that this idea was very efficient for the time-dependent Schrödinger equation with three [18][19][20][21][22][23][24] and four [1,4,25,26] nonseparable spatial variables. Especially productive it was for the quantitative description of the ultracold pair atomic collisions in confined geometry of harmonic traps (low-dimensional few-body systems) [1,2,4,9,[25][26][27]. Our idea on how to construct the two-dimensional DVR [18,19] was also successfully extended for computing vibrational levels of four-atom molecules [16] where it has got the title non-direct product DVR (npDVR). An alternative variant of 2D npDVR with non-direct product angular grid coinciding on the unit sphere with the nodes of the Lebedev quadratures was suggested in [17].
We discuss here the key elements of the 2D npDVR by using as an example the three-dimensional time-dependent Schrödinger equation describing the pair atomic collision in the harmonic waveguide U(r) is the relative interatomic variable and ω ix ω iy in the general case of the anisotropic waveguide. We seek the solution ψ(r, t) as the expansion [1,[18][19][20] ψ(r, with respect to the two-dimensional basis defined on an angular grid Ω j = (θ j θ , φ j φ ). For the θ variable, the N θ mesh points θ j θ correspond to the zeros of the Legendre polynomial P N θ (cos θ j θ ) = 0. For the φ variable, the N φ mesh points are chosen is equal to the number of basis functions in the expansion (1) and the number of terms in the definition (2), where the symbol ν represents the twofold index {l, m} and the sum over ν is equivalent to the double sum The l and m indexes show the number of zeros over the θ and φ variables of the polynomials Y ν (Ω) which we specify in the next paragraph. Due to the definition (3) only odd values N φ may be chosen. N θ can take on arbitrary values. The coefficients (Y −1 ) ν j in the definition (2) are the elements of the (2) belongs to the class of Lagrange functions [15]) and the coefficients ψ j (r, t) in the expansion (1) define the values of the searching solution ψ(r, t) at the points of the angular grid Ω j : where C l l = δ ll holds in general, and thus Y ν (Ω) coincides with the usual spherical harmonic with a few possible exceptions for large values of ν such that we receive the orthogonality relation for all ν and ν ≤ N. Here the N weights λ j are the standard Gauss-Legendre weights multiplied by 2π/N φ . For most ν and ν the above relation is automatically satisfied because the basis functions Y ν (Ω) are orthogonal and the Gaussian quadrature is exact. For these ν we have C l l = δ ll in Eq. (4). However, in a few cases involving the highest l values, some polynomials P m l (θ) have to be orthogonalized in the sense of the Gaussian quadrature (C l l δ ll for these specific values of l). With this choice, the matrix It provides the unitary transformation from the ndDVR (2) to the partial-wave representation and back This peculiarity (orthogonality on the grid Ω j ) of the basis (2) was used for constructing the efficient splitting-up method for 3D and 4D time-dependent Schrödinger equations in 2D npDVR [18][19][20]. In the next section we present such a scheme for the 3D time-dependent Schródinger equation.

Splitting-up Method for 3D Time-Dependent Schrödinger Equation
In the 2D npDVR (2) the only nondiagonal part of the Hamiltonian H(r) is the angular part of the kinetic energy operator In our scheme we neglect the terms ∼ 1/(l max (l max , which are of the order of the errors induced by the finite difference approximation over the radial variable and the splitting up procedure. It gives us the possibility to perform the diagonalization of the angular part of the kinetic energy operator by using the unitary transformation S jν = λ 1/2 j Y jν defined in the previous section. This has been exploited for developing an efficient algorithm with a computational time scaling proportional to the number N = N θ × N φ of the unknowns in the system of equations ( = 1) [1,19,20,24] i approximating the initial 3D Schrödinger equation at the 2D npDVR (2). The system of equations (6) arises in the description of pair collisions of the ultracold identical atoms (m 1 = m 2 = m) in optical waiveguides when ω 1 = ω 2 = ω and the center-of-mass of the pair can be separated [1]. For the propagation ψ j (r, t n ) → ψ j (r, t n+1 ) in time t n → t n+1 = t n + Δt we have developed a computational scheme [1,18,20,24] based on the component-by-component split-operator method suggested by G.I. Marchuk in 1971 [28]. The HamiltonianĤ(r) permits the splitting into the following two parts H j j (r) = h j j (r) + U j (r)δ j j where h j j (r) = − δ j j 2μ Here the diagonal in the npDVR potentials V(r) and U(r) describe the spherically symmetric interatomic interaction and the interaction of the atoms with the trap, , and μ = m/2. Subsequently, the time-step ψ j (r, t n ) → ψ j (r, t n+1 = t n + Δt) was approximated according to The time evolution is described in detail in [1,20,24] where the fact that the npDVR (functions (2)) gives the diagonal representation forÛ and the Y ν -representation (functions (4)) gives the diagonal representations forĥ was used for constructing an efficient computational algorithm. The time evolution was accomplished as follows. For the first and the last steps according to the relation (8) we write the function ψ and the operators exp(−iΔtÛ/2) in our 2D npDVR (2) on the grid Ω j = (θ j θ , φ j φ ). Since the potentialÛ(r) = U j (r)δ j j is diagonal in this representation the first and the last steps represent simple multiplications of the diagonal matrices exp(− i 2 ΔtU j (r))δ j j . The intermediate step in (8) EPJ Web of Conferences 01008-p.4 depending onĥ is treated in the basis Y ν (4) where the matrix operatorĥ(r) is diagonal with respect to the indices m and l (ν = (l, m)). For that we approximate the exponential operators according to which ensures the desired accuracy of the numerical algorithm (8). Thus, after the discretization of r with the help of finite-differences the matrixĥ possesses a band structure and we arrive at the following boundary-value problems which can be solved rapidly due to the band structure of the matrixĥ. This computational scheme is unconditionally stable [28], preserves unitarity and is very efficient, i.e. the computational time is proportional to the total number N of the grid points and the basis functions in (1): see Fig. 1 where the comparison with the conventional Crank-Nicolson algorithm demonstrates the advantage of our scheme. The efficiency of the computational procedure is based on the fast transformation carried out with help of the unitary matrix S jν = λ 1/2 j Y jν between the two relevant representations: 2D npDVR (2) and the Y ν -representation (4). Another attractive peculiarity of the suggested scheme (7,8) is the supposed transparent procedure for its parallelization. Actually, the most time-consuming part, the integration of the system of N uncoupled differential equations (10), permits a simple procedure of it parallelization: the system is diagonal, therefore every equation of the system can be integrated independently on a separate processor. It makes the method very promising in application to other actual problems of low-dimensional few-body physics which we briefly discuss in the conclusion.

Shifts and Widths of Feshbach Resonances in Atomic Waveguides
In this section we give an example of application of the npDVR for the stationary four channel scattering problem in a few dimensions. Such a problem arises in the attempt to describe the shifts and Mathematical Modeling and Computational Physics 2015 01008-p.5 widths of the magnetic Feshbach resonance on the pair atomic collision in confined geometries of tight traps [7,8].
Already at the initial stage of the investigation of ultracold pair collisions in atomic traps, it was shown that the confining geometry can drastically change the scattering properties of the ultracold atoms and induces resonances in the collisions (confinement-induced resonances (CIRs)) [29]. Further research in this field has shown that a necessary ingredient for the experimental realization of the CIR in a confining trap is the existence of magnetic Feshbach resonances in the free-space (see for example [6]). However, until recently, only single-channel potential models with zero-energy bound states were used for simulating the magnetic Feshbach resonances in the 3D free-space.
In our work [7] the frame of the single-channel theoretical approaches for the CIRs was overcome and effects coming from different rotational structure of the Feschach resonance, the resonance width and the background scattering have been taken into account. This was possible thanks to our multichannel approach for CIRs and atomic transitions in the waveguides in the multimode regime [2]. We have replaced in this scheme the single-channel scalar interatomic interaction by the four-channel tensorial potentialV(r) modeling resonances of broad, narrow and overlapping character according to the two-channel parametrization of A.D. Lange et. al. [30].
In this approach the scattering amplitude f (k, B) describing a pair collision in the free-space (ω ⊥ = 0) is calculated for different magnetic field strengths B and varying parameters of the potentialV by solving the Schrödinger equation with the scattering boundary conditions at kr → ∞ for the fixed E → 0 (k = 2μE/ → 0) [7]. The diagonal matrixB in (11) is defined as B ii = δμ i (B − B i ) (here δμ i is the relative magnetic moment of the channel i, i = 1, 2, 3) and B ee = 0. After separation of the angular part in (11) we come to the system of four coupled radial equations for the radial part φ α (r) of the desired wave function |ψ = α ψ α (r)|α = α φ α (r)Y l α 0 (r)|α , where α = {e, i = 1, 2, 3}, (l e = 0, l 1 = 0, l 2 = 2, l 3 = 4) and the matrix elements V αβ (r) are defined by four diagonal parameters V c,i , V e and three non-diagonal ones Ω i [7]. The centrifugal barrier 2 l α (l α +1) 2μr 2 in (12) models at r → 0 the correct asymptotic behavior φ i (r) ∼ r l i +1 of the molecular bound states |c, i in the closed channels which couple to the entrance s-wave (l e = 0) channel |e by the nondiagonal terms V αβ (r)(α β).
By varying the V c,i , V e and Ω i we have obtained excellent agreement of the calculated s-wave scattering length a(B) = − f (k → 0, B) with the analytical results reported in [30] for Cs atoms in the hyperfine state |F = 3, m F = 3 for the considered magnetic field regime −40G < B < 60G, corresponding to the conditions of Innsbruck experiment [6]. In this regime, we observe three resonances B i in a(B), which correspond to the coupling to the s-, d-, and g-wave molecular states in the closed channels |c, i .
Then, after finding the parameters of the interparticle interactionV(r), we have analyzed the scattering properties of the s-, d-and g-wave magnetic Feshbach resonances in harmonic waveguides by EPJ Web of Conferences 01008-p. 6 integrating the Schrödinger equation (11) for ω ⊥ 0 with the scattering boundary conditions ψ e (r) = cos(k || z) + f e (B) exp{ik || | z |} Φ 0 (ρ) , ψ c,i (r) → 0 (13) at | z |=| r cos θ |→ ∞ adopted for a confining trap [7]. f e (B, E) is the scattering amplitude, corresponding to symmetric states with respect to the exchange z → −z (we consider collisions of identical bosonic Cs atoms), Φ 0 (ρ) is the wave function of the ground-state of the two-dimensional harmonic oscillator and k || = 2μ(E − ω ⊥ )/ = 2μE / . In the presence of the harmonic trap (ω ⊥ 0) the problem (11), (13) becomes non-separable in the plane (ρ, z) = (r, θ), i.e. the azimuthal angular part is separated in the wave function |ψ and Eq. (11) is reduced to a coupled system of four 2D Schrödinger-type equations. To integrate this coupled channel 2D scattering problem in the plane {r, θ} we have applied an extended version of the computational scheme of [2]. The computations have been performed in a range of variation of ω ⊥ close to the experimental values of the trap frequencies ω ⊥ ∼ 2π × 14.5 kHz [6]. We have integrated Eq. (11) for varying B and fixed longitudinal colliding energy E || = E − ω ⊥ → 0. At input, the experimentally known parameters of the Feshbach resonances in the absence of the waveguide [30], namely the resonant energies E c,i (or the corresponding values of the field strengths B i of the external magnetic field), the widths of the resonances Δ i (Γ i ), the spin characteristics and the background scattering length a bg , were used. The output of the model is the calculated scattering amplitude f e (B) in the confining trap ω ⊥ and the transmission T(B) =| 1 + f e (B) | 2 .
We have calculated the shifts and widths of s-, d-and g-wave magnetic Feshbach resonances of Cs atoms emerging in harmonic waveguides as CIRs (the minimums in the transmission T(B)) and resonant enhancement of transmission T(B) at zeros a(B) = 0 of the free-space scattering length. The results are illustrated by the curves T(B) calculated near the d-wave resonance for different ω ⊥ (figure 2) corresponding to the transverse frequencies of the optical trap being used in the experiment [6]. It is shown that the position of B min of the transmission coefficient T(B) minimum depends on the transverse trap width a ⊥ = /(μω ⊥ ) according to the simplified Olshanii model [29]. Actually, the corresponding scattering length a(B min ) at the point B min is accurately described by the formulas a(B min ) = a ⊥ /1.48 obtained in [29] for the position of the CIR in the s-wave case. This property was experimentally confirmed for d-wave Feshbach resonances in a gas of Cs atoms [6]. The

Conclusion
We have considered computational aspects of the mathematical models developed by us for ultracold two-body resonant processes in atomic traps. We have discussed in detail the key element of the elaborated computational schemes, the npDVR, and two examples of its efficient application to the time-dependent and the stationary Schrödinger equations with a few spatial variables.
Within these approaches, several novel effects for pair collisions in tight atomic waveguides were found: the s-and p-wave CIRs in multimode regimes [1,2,25,26] and resonant molecule formation with transferring energy relies to center-of-mass excitation while forming molecules [4]. The last effect was recently been confirmed in the Heidelberg experiment [5]. Our calculations have also been used for planning and the interpretation of the Innsbruck experiment on the investigation of CIRs in the ultracold Cs gas [6]. Recently, we have calculated the Feshbach resonance shifts and widths induced by atomic waveguides [7,8] and predicted dipolar CIRs [9].
The efficiency of our splitting-up method for the time-dependent Schrödinger equation and the supposed transparent procedure for its parallelization make the method very promising in applications to other actual problems of the low-dimensional few-body physics. Note here the problem of the extension to the ultracold atomic collisions in quasi-2D geometry of confining traps [27] needed for resolving the "puzzle" of the 2D CIR [6], the collisional three-body problem (Efimov resonances) in tight traps, and the non-linear time-dependent Schrödinger equation with a few spatial variables arising in the physics of the Bose-Einstein condensates.