Energy spectrum of a harmonically trapped two-atom system with spin-orbit coupling

Ultracold atomic gases provide a novel platform with which to study spin-orbit coupling, a mechanism that plays a central role in the nuclear shell model, atomic fine structure and two-dimensional electron gases. This paper introduces a theoretical framework that allows for the efficient determination of the eigenenergies and eigenstates of a harmonically trapped two-atom system with short-range interaction subject to an equal mixture of Rashba and Dresselhaus spin-orbit coupling created through Raman coupling of atomic hyperfine states. Energy spectra for experimentally relevant parameter combinations are presented and future extensions of the approach are discussed.

much larger than that in the direction where the spinorbit coupling term acts. This assumption reduces the problem to an effective one-dimensional Hamiltonian in the x-coordinates with effective 1D two-body interaction. The relationship between the true 3D atom-atom and effective 1D atom-atom interaction has been derived in Refs. [46][47][48]. We find analytical solutions to the twoatom system for arbitrary spin-orbit coupling strength and scattering length and vanishing Raman coupling strength. The case of non-zero Raman coupling strength is treated by expanding the system Hamiltonian in terms of the eigenstates for vanishing Raman coupling strength. We find that the relevant Hamiltonian matrix elements have closed analytical expressions, leaving the matrix diagonalization as the only numerical step. The developed framework can, as discussed toward the end of our paper, be readily generalized to a spherically-symmetric harmonic trap or an axisymmetric trap. Moreover, the framework developed also lays the groundwork for treating dynamical aspects of trapped two-body systems with non-vanishing spin-orbit and Raman coupling strengths and for treating the corresponding three-body system.
We consider two structureless one-dimensional particles of mass m subject to a single-particle spin-orbit coupling term of strength k so , a Raman coupling term with strength Ω, detuning δ, and an external harmonic potential with angular trapping frequency ω. For k so = Ω = δ = 0, the two-particle Hamiltonian is given by H sr , where x j denotes the position coordinate of the jth particle and V 2b the short-range interaction potential. For non-zero k so , Ω and δ, the two-particle Hamiltonian is given by H, where σ (j) x and σ (j) y denote Pauli matrices,Î the identity matrix and p xj the momentum of the jth particle. In the following, we first derive solutions to the timeindependent Schrödinger equation governed by H with Ω = 0 and then discuss how to obtain the solutions for non-zero Ω.
To determine the eigenstates and eigenenergies of H, we perform a rotation in spin space [49]. Specifically, we defineH via a unitary transformation of H,H = U † HU , where U = exp[ı(σ (1) x + σ (2) x )π/4]. The eigenenergies of the Hamiltonian H andH coincide while the eigenstates Ψ of the Hamiltonian H are related to the eigenstatesΨ of the HamiltonianH through Ψ = UΨ. A straightforward calculation shows that U † σ (j) x and σ (j) z for q = x and y, respectively. Correspondingly, we havẽ For Ω = 0,H is diagonal in the pseudo-spin basis | ↑ 1 | ↑ 2 , | ↑ 1 | ↓ 2 , | ↓ 1 | ↑ 2 and | ↓ 1 | ↓ 2 with diagonal elementsH ↑↑ ,H ↑↓ ,H ↓↑ andH ↓↓ . To find the corresponding eigenstates, we approximate the two-body interaction by a delta-function interaction with coupling constant g, V 2b (x 1 − x 2 ) = gδ(x 1 − x 2 ). For this interaction model, the eigenenergies and eigenstates of H sr are known in compact form [8]. States that are even in the relative coordinate are affected by the coupling constant g while those that are odd in the relative coordinate are not. For states that are even in the relative coordinate, the eigenenergies E sr nq of H sr (see solid lines in Fig. 1 for the n = 0 energies) are given by (n + 2q + 1) ω, where the center of mass quantum number n takes the values n = 0, 1, · · · and the non-integer quantum number q is determined by the transcendental equation [8] 2Γ here, a ho denotes the harmonic oscillator length, a ho = /(mω). The corresponding eigenfunctions ψ sr nq (x, X) are given by φ q (x)Φ n (X), where the relative and center of mass coordinates are defined through x = (x 1 − x 2 )/ √ 2 and X = (x 1 + x 2 )/ √ 2, respectively. The relative wave functions φ q (x) can be written in terms of the confluent hypergeometric function U [8], where N q denotes a normalization constant. The center of mass functions Φ n (X) are given by the one-dimensional harmonic oscillator functions for a mass m particle, where H n denotes the Hermite polynomial of order n and N ni n = ( √ π2 n n!a ho ) −1/2 . For states that are odd in the relative coordinate, the eigenenergies E sr nq of H sr (see dotted lines in Fig. 1 for the n = 0 energies) are given by (n + 2q + 2) ω, where q and n take the values 0, 1, · · · . In this case, the eigenfunctions ψ sr nq (x, X) are simply products of the non-interacting harmonic oscillator functions in x and X.
In addition to using the known properties of H sr , we take advantage of the fact that the kinetic energy (p 2 x1 + p 2 x2 )/(2m) of H sr and the k so dependent terms can be combined, This identity suggests that the momentum-dependent spin-orbit coupling terms add a "momentum boost" to the solutions ψ sr nq (x, X) of H sr . Indeed, it is readily verified that the eigenstates ofH ↑↑ ,H ↑↓ ,H ↓↑ , andH ↓↓ are given bỹ respectively. For fixed g and n and vanishing δ, the states given in Eqs. (6)-(9) are degenerate with eigenenergies E nq = E sr nq − 2 k 2 so /m. For |g| = ∞, the degeneracy doubles (see the crossings of the solid and dotted lines in Fig. 1) since q takes the values 1/2, 3/2, · · · for ψ sr nq that are even in x and the values 0, 1, 2, · · · for ψ sr nq that are odd in x, i.e., since each of the ψ sr nq odd in x is degenerate with one of the ψ sr nq even in x. For non-vanishing δ, the energies are shifted by δ, 0, 0 and −δ, respectively. The eigenenergies are simply the sum of a term that depends on the coupling constant g, a center of mass contribution that is characterized by n, a term that depends on the square of the spin-orbit coupling strength k so and a term that depends on the detuning δ.
The integral I σ ′ 1 σ ′ 2 ,σ1σ2 q ′ q over the relative coordinate can be performed by expanding [φ q ′ (x)] * and φ q (x) in terms of the non-interacting harmonic oscillator functions φ ni l (x), φ q (x) = lim lmax→∞ lmax l=0 c (q) l φ ni l (x), where the expansion coefficients c (q) l can be obtained analytically [8]. The integral I σ ′ 1 σ ′ 2 ,σ1σ2 q ′ q then becomes a double sum over integrals that have the same structure as the center of mass integrals J σ ′ 1 σ ′ 2 ,σ1σ2 n ′ n . In the calculations reported below, we use a finite cutoff l max . The "optimal" cutoff depends on the value of g considered, the number of relative functions φ q (x) included in the basis and the desired accuracy. For |g| = ∞, we find, as in the g = 0 case, a closed analytical expression for the integral I . Having analytical expressions for the matrix elements ofH, the eigenenergies can be obtained through matrix diagonalization.
To obtain basis functions with good quantum numbers, we work with linear combinations of the functions given in Eqs. (6)-(9), i.e., we work with the basis func- By properly combining the parts of ψ sr nq that are even or odd in the relative coordinate and even or odd in the center of mass coordinate, we construct basis functions that are eigenstates of the operators P 12 and Y 12 . The operator P 12 exchanges the coordinates (position and spin) of particles 1 and 2. Basis functions that are unchanged under the operation P 12 are needed to describe states with bosonic symmetry (p 12 = +1) and those that pick up a minus sign under the operation P 12 are needed to describe states with fermionic symmetry (p 12 = −1). The operator Y 12 can be written as σ (1) x σ (2) x P P 12 , where the parity operator P changes x j to −x j (j = 1 and 2). The Y 12 operator determines the "helicity" of the system. We label the eigenstates by (p 12 , y 12 ), where p 12 = ±1 and y 12 = ±1 are defined by their actions on an eigenstate. The basis functions with (+1, +1) symmetry, for example, are given byψ X,+ with φ q (x) even and Φ n (X) even, byψ X,− with φ q (x) even and Φ n (X) odd, byψ x,+ with φ q (x) even and Φ n (X) even, and byψ x,− with φ q (x) odd and Φ n (X) even.
As an example, Figs. 2(a)-2(c) show energy spectra corresponding to eigenstates with (p 12 , y 12 ) = (+1, +1) as a function of the Raman coupling strength Ω for vanishing detuning δ, small coupling constant g, g ≈ 0.1414a ho ω, and three different spin-orbit coupling strengths k so . The energy spectrum in Fig. 2(a) is, to leading order, given by the spectrum for δ = k so = g = 0. In this limiting case, the energies are equal to (2j + 1) ω ± Ω and 2j ω, where j = 0, 1, · · · . Finite g and k so values introduce shifts and avoided crossings. Specifically, the small positive coupling constant g introduces a positive energy shift for the states that are even in the relative coordinate, which-in first-order perturbation theory-is given by 1 √ 2π (2q)! (q!) 2 4 q g/a ho . The spin-orbit coupling term introduces, in the small Ω regime, a small down shift that is proportional to k 2 so . This down shift is negligible in Fig. 2(a) but clearly visible in Figs. 2(b) and 2(c). Moreover, the spin-orbit coupling introduces avoided crossings. The broadest avoided crossings occur around Ω = ω, where states with the same q but n quantum numbers that differ by one are coupled. The reason is that the spin-orbit coupling term connects, via the total momentum operator, states in first-order perturbation theory if the states' n quantum numbers differ by one and in higher-order perturbation theory otherwise.
Figures 2(d)-2(f) show energy spectra for the strong coupling limit, i.e., for |g| → ∞, as a function of the Raman coupling Ω for vanishing detuning δ. To facilitate the comparison between the large and small g limits, the spin-orbit coupling strengths k so in Figs. 2(d)-2(f) are the same as in Figs. 2(a)-2(c). In the regime where Ω ≪ 2 k 2 so /m, the energies change approximately linearly with Ω (with positive, vanishing or negative slope).
When Ω ≫ 2 k 2 so /m, the low-lying portion of the energy spectrum consists of approximately parallel energy levels that can be parameterized as c−Ω, where c is a constant.
As already aluded to in the introduction, the theoretical framework developed can be generalized to higherdimensional trapping geometries. For a spherically symmetric 3D system, e.g., the eigenstates of the 3D Hamiltonian H sr with 3D contact interaction can be expanded in terms of products of 2D and 1D harmonic oscillator states using cylindrical coordinates. As in the 1D case pursued in this work, the matrix elements for the higherdimensional system can be calculated analytically. Axisymmetric harmonic traps with spin-orbit coupling in one direction can be treated analogously. Furthermore, using the eigenstates of the trapped three-particle system in 1D, 2D or 3D with contact interactions [11,22,23] and expressing the three-particle Hamiltonian in terms of the eight pseudo-spin states, a non-zero Ω introduces off-diagonal elements that can be calculated analytically following steps similar to those discussed in this paper.
Summarizing, this work introduced a theoretical framework that allows for the efficient determination of the energy spectrum and eigenstates of the trapped twoparticle system in 1D with contact interaction and spinorbit and Raman coupling terms. The energy spectra show a rich dependence on the interaction, spin-orbit and Raman coupling strengths. The framework presented provides an important stepping stone for treating more complicated systems with spin-orbit coupling, such as higher-dimensional two-body systems and three-body systems.