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

Over the past decade, much progress has been made in preparing isolated ultracold few-atom systems experimentally [1–4]. Moreover, a variety of tools for manipulating and probing such systems have been developed. On the theoretical side, a number of analytical and numerical approaches have been developed [5–13]. A large number of analytical treatments approximate the true alkali atom–alkali atom potential by a zero-range potential [8, 14–16]. This replacement yields reliable results in the low-energy regime where the de Broglie wave length is larger than the van der Waals length. For example, using zero-range contact interactions, the energy spectrum of two harmonically trapped atoms has been determined analytically [8–10]. These two-body solutions are available in 1D, 2D and 3D [8], and have played a vital role in guiding and interpreting experiments [17–19] as well as in theoretical studies of the two-body dynamics [20, 21] and of larger harmonically trapped systems [11–13, 22–24].

Recently, synthetic gauge fields, which allow for the realization of Hamiltonians that contain spin–orbit coupling terms, have been realized experimentally [25–35]. The purpose of this paper is to address how the trapped two-particle spectrum, obtained by modeling the two-body interaction via a zero-range δ-function, changes in the presence of spin–orbit and Raman coupling. While the two-particle system with spin–orbit coupling in free space [36–39] as well as the trapped single-particle system with spin–orbit coupling [40, 41] have received considerable attention, little is known about the trapped two-particle system with spin–orbit coupling and two-body interaction [42, 43]. In going from the trapped single-atom to the trapped two-atom system, a new length scale, i.e., the atom–atom scattering length, comes into play. Thus, an interesting question concerns the interplay between the interaction energy and the energy scales associated with the spin–orbit and Raman coupling strengths.

Our framework applies to the situation where the spin–orbit (or more precisely, spin–momentum) coupling term acts, as in recent experiments [29–35], along one spatial direction, say the x-direction. This corresponds to an equal mixture of Rashba and Dresselhaus spin–orbit coupling [44, 45]. For simplicity, we assume that the harmonic confinement in the other two spatial directions is much larger than that in the direction where the spin–orbit 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 [46–48]. We find analytical solutions to the two-atom 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 expressing 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 with strength ${{k}_{{\rm so}}}$, a Raman coupling term with strength Ω, detuning δ, and an external harmonic potential with angular trapping frequency ω. For ${{k}_{{\rm so}}}=\Omega =\delta =0$, the two-particle Hamiltonian is given by ${{H}_{{\rm sr}}}$,

$$\begin{eqnarray}&&{{H}_{{\rm sr}}}=\mathop{\sum }\limits_{j=1}^{2}\left( \frac{-{{\hbar }^{2}}}{2\;m}\frac{{{\partial }^{2}}}{\partial x_{j}^{2}}+\frac{1}{2}m{{\omega }^{2}}x_{j}^{2} \right)+{{V}_{2{\rm b}}}({{x}_{1}}-{{x}_{2}}),\end{eqnarray} \tag{ 1 }$$

where x_j denotes the position coordinate of the jth particle and ${{V}_{2{\rm b}}}$ the short-range interaction potential. For non-zero ${{k}_{{\rm so}}}$, Ω and δ, the two-particle Hamiltonian is given by H,

$$\begin{eqnarray}&&H={{H}_{{\rm sr}}}\hat{I}+\mathop{\sum }\limits_{j=1}^{2}\left[ \frac{\hbar {{k}_{{\rm so}}}}{m}{{p}_{xj}}\sigma _{y}^{(j)}+\frac{\Omega }{2}\sigma _{x}^{(j)}+\frac{\delta }{2}\sigma _{y}^{(j)} \right],\end{eqnarray} \tag{ 2 }$$

where $\sigma _{x}^{(j)}$ and $\sigma _{y}^{(j)}$ denote Pauli matrices, $\hat{I}$ the identity matrix and p_xj the momentum of the jth particle. In the following, we first derive solutions to the time-independent Schrödinger equation governed by H with $\Omega =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 define $\tilde{H}$ via a unitary transformation of H, $\tilde{H}={{U}^{\dagger }}HU$, where $U={\rm exp} [\imath (\sigma _{x}^{(1)}+\sigma _{x}^{(2)})\pi /4]$. The eigenenergies of the Hamiltonian H and $\tilde{H}$ coincide while the eigenstates Ψ of the Hamiltonian H are related to the eigenstates $\tilde{\Psi }$ of the Hamiltonian $\tilde{H}$ through $\Psi =U\tilde{\Psi }$. A straightforward calculation shows that ${{U}^{\dagger }}\sigma _{q}^{(j)}U=\sigma _{x}^{(j)}$ and $\sigma _{z}^{(j)}$ for q = x and y, respectively. Correspondingly, we have

$$\begin{eqnarray}&&\tilde{H}={{H}_{{\rm sr}}}\hat{I}+\mathop{\sum }\limits_{j=1}^{2}\left[ \frac{\hbar {{k}_{{\rm so}}}}{m}{{p}_{xj}}\sigma _{z}^{(j)}+\frac{\Omega }{2}\sigma _{x}^{(j)}+\frac{\delta }{2}\sigma _{z}^{(j)} \right].\end{eqnarray} \tag{ 3 }$$

For $\Omega =0$, $\tilde{H}$ is diagonal in the pseudo-spin basis $|\uparrow {{\rangle }_{1}}|\uparrow {{\rangle }_{2}}$, $|\uparrow {{\rangle }_{1}}|\downarrow {{\rangle }_{2}}$, $|\downarrow {{\rangle }_{1}}|\uparrow {{\rangle }_{2}}$ and $|\downarrow {{\rangle }_{1}}|\downarrow {{\rangle }_{2}}$ with diagonal elements ${{\tilde{H}}^{\uparrow \uparrow }}$, ${{\tilde{H}}^{\uparrow \downarrow }}$, ${{\tilde{H}}^{\downarrow \uparrow }}$ and ${{\tilde{H}}^{\downarrow \downarrow }}$. To find the corresponding eigenstates, we approximate the two-body interaction by a delta-function interaction with coupling constant g, ${{V}_{2{\rm b}}}({{x}_{1}}-{{x}_{2}})=g\delta ({{x}_{1}}-{{x}_{2}})$. For this interaction model, the eigenenergies and eigenstates of ${{H}_{{\rm 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_{nq}^{{\rm sr}}$ of ${{H}_{{\rm sr}}}$ (see solid lines in figure 1 for the n = 0 energies) are given by $(n+2q+1)\hbar \omega$, where the center of mass quantum number n takes the values $n=0,1,\cdots$ and the non-integer quantum number q is determined by the transcendental equation [8]

$$\begin{eqnarray}&&\frac{2\Gamma (-q+1/2)}{\Gamma (-q)}=-\frac{g}{\sqrt{2}\hbar \omega {{a}_{{\rm ho}}}};\end{eqnarray} \tag{ 4 }$$

here, ${{a}_{{\rm ho}}}$ denotes the harmonic oscillator length, ${{a}_{{\rm ho}}}=\sqrt{\hbar /(m\omega )}$. The corresponding eigenfunctions $\psi _{nq}^{{\rm sr}}(x,X)$ are given by ${{\phi }_{q}}(x){{\Phi }_{n}}(X)$, where the relative and center of mass coordinates are defined through $x=({{x}_{1}}-{{x}_{2}})/\sqrt{2}$ and $X=({{x}_{1}}+{{x}_{2}})/\sqrt{2}$, respectively. The relative wave functions ${{\phi }_{q}}(x)$ can be written in terms of the confluent hypergeometric function U [8],

$$\begin{eqnarray*}&&{{\phi }_{q}}(x)={{N}_{q}}U{{\left( -q,1/2,(x/{{a}_{{\rm ho}}} \right)}^{2}}){\rm exp} [-{{x}^{2}}/(2a_{{\rm ho}}^{2})],\end{eqnarray*}$$

where N_q denotes a normalization constant. The center of mass functions ${{\Phi }_{n}}(X)$ are given by the one-dimensional harmonic oscillator functions for a mass m particle, ${{\Phi }_{n}}(X)=N_{n}^{{\rm ni}}{{H}_{n}}(X/{{a}_{{\rm ho}}}){\rm exp} [-{{X}^{2}}/(2a_{{\rm ho}}^{2})]$, where H_n denotes the Hermite polynomial of order n and $N_{n}^{{\rm ni}}={{\left( \sqrt{\pi }{{2}^{n}}n!{{a}_{{\rm ho}}} \right)}^{-1/2}}$. For states that are odd in the relative coordinate, the eigenenergies $E_{nq}^{{\rm sr}}$ of ${{H}_{{\rm sr}}}$ (see dotted lines in figure 1 for the n = 0 energies) are given by $(n+2q+2)\hbar \omega$, where q and n take the values $0,1,\cdots$. In this case, the eigenfunctions $\psi _{nq}^{{\rm sr}}(x,X)$ are simply products of the non-interacting harmonic oscillator functions in x and X.

Zoom In Zoom Out Reset image size

In addition to using the known properties of ${{H}_{{\rm sr}}}$, we take advantage of the fact that the kinetic energy $(p_{x1}^{2}+p_{x2}^{2})/(2\;m)$ of ${{H}_{{\rm sr}}}$ and the ${{k}_{{\rm so}}}$ dependent terms can be combined,

$$\begin{eqnarray}&&\frac{p_{xj}^{2}}{2\;m}\pm \frac{\hbar {{k}_{{\rm so}}}}{m}{{p}_{xj}}=\frac{{{\left( {{p}_{xj}}\pm \hbar {{k}_{{\rm so}}} \right)}^{2}}}{2\;m}-\frac{{{\hbar }^{2}}k_{{\rm so}}^{2}}{2\;m}.\end{eqnarray} \tag{ 5 }$$

This identity suggests that the momentum-dependent spin–orbit coupling terms add a 'momentum boost' to the solutions $\psi _{nq}^{{\rm sr}}(x,X)$ of ${{H}_{{\rm sr}}}$. Indeed, it is readily verified that the eigenstates of ${{\tilde{H}}^{\uparrow \uparrow }}$, ${{\tilde{H}}^{\uparrow \downarrow }}$, ${{\tilde{H}}^{\downarrow \uparrow }}$, and ${{\tilde{H}}^{\downarrow \downarrow }}$ are given by

$$\begin{eqnarray}&&\tilde{\psi }_{nq}^{\uparrow \uparrow }(x,X)={\rm exp} (-\imath \sqrt{2}{{k}_{{\rm so}}}X)\psi _{nq}^{{\rm sr}}(x,X),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\tilde{\psi }_{nq}^{\uparrow \downarrow }(x,X)={\rm exp} (-\imath \sqrt{2}{{k}_{{\rm so}}}x)\psi _{nq}^{{\rm sr}}(x,X),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\tilde{\psi }_{nq}^{\downarrow \uparrow }(x,X)={\rm exp} (\imath \sqrt{2}{{k}_{{\rm so}}}x)\psi _{nq}^{{\rm sr}}(x,X),\end{eqnarray} \tag{ 8 }$$

and

$$\begin{eqnarray}&&\tilde{\psi }_{nq}^{\downarrow \downarrow }(x,X)={\rm exp} (\imath \sqrt{2}{{k}_{{\rm so}}}X)\psi _{nq}^{{\rm sr}}(x,X),\end{eqnarray} \tag{ 9 }$$

respectively. For fixed g and n and vanishing δ, the states given in equations (6)–(9) are degenerate with eigenenergies ${{E}_{nq}}=E_{nq}^{{\rm sr}}-{{\hbar }^{2}}k_{{\rm so}}^{2}/m$. For $|g|=\infty$, the degeneracy doubles (see the crossings of the solid and dotted lines in figure 1) since q takes the values $1/2,3/2,\cdots$ for $\psi _{nq}^{{\rm sr}}$ that are even in x and the values $0,1,2,\cdots$ for $\psi _{nq}^{{\rm sr}}$ that are odd in x, i.e., since each of the $\psi _{nq}^{{\rm sr}}$ odd in x is degenerate with one of the $\psi _{nq}^{{\rm sr}}$ even in x. For non-vanishing δ, the energies are shifted by δ, 0, 0 and $-\delta$, 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}_{{\rm so}}}$ and a term that depends on the detuning δ.

If Ω is non-zero, the Hamiltonian $\tilde{H}$ expressed in the $|\uparrow {{\rangle }_{1}}|\uparrow {{\rangle }_{2}}$, $|\uparrow {{\rangle }_{1}}|\downarrow {{\rangle }_{2}}$, $|\downarrow {{\rangle }_{1}}|\uparrow {{\rangle }_{2}}$, and $|\downarrow {{\rangle }_{1}}|\downarrow {{\rangle }_{2}}$ pseudo-spin basis is no longer diagonal. To determine the eigenenergies and eigenstates for non-zero Ω, we express $\tilde{H}$ in terms of the eigenstates $\tilde{\psi }_{nq}^{{{\sigma }_{1}}{{\sigma }_{2}}}$, where σ₁ and σ₂ take the values $\uparrow$ and $\downarrow$. The off-diagonal matrix elements $H_{n^{\prime} q^{\prime} ,nq}^{{{\sigma }_{1}}^{\prime} {{\sigma }_{2}}^{\prime} ,{{\sigma }_{1}}{{\sigma }_{2}}}$,

$$\begin{eqnarray}&&H_{n^{\prime} q^{\prime} ,nq}^{{{\sigma }_{1}}^{\prime} {{\sigma }_{2}}^{\prime} ,{{\sigma }_{1}}{{\sigma }_{2}}}=\frac{\Omega }{2}\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }{{\left( \tilde{\psi }_{n^{\prime} q^{\prime} }^{{{\sigma }_{1}}^{\prime} {{\sigma }_{2}}^{\prime} } \right)}^{*}}\tilde{\psi }_{nq}^{{{\sigma }_{1}}{{\sigma }_{2}}}{\rm d}x{\rm d}X\end{eqnarray} \tag{ 10 }$$

can be separated into two one-dimensional integrals,

$$\begin{eqnarray}&&H_{n^{\prime} q^{\prime} ,nq}^{{{\sigma }_{1}}^{\prime} {{\sigma }_{2}}^{\prime} ,{{\sigma }_{1}}{{\sigma }_{2}}}=\frac{\Omega }{2}I_{q^{\prime} q}^{{{\sigma }_{1}}^{\prime} {{\sigma }_{2}}^{\prime} ,{{\sigma }_{1}}{{\sigma }_{2}}}J_{n^{\prime} n}^{{{\sigma }_{1}}^{\prime} {{\sigma }_{2}}^{\prime} ,{{\sigma }_{1}}{{\sigma }_{2}}},\end{eqnarray} \tag{ 11 }$$

where

$$\begin{eqnarray}&&I_{q^{\prime} q}^{{{\sigma }_{1}}^{\prime} {{\sigma }_{2}}^{\prime} ,{{\sigma }_{1}}{{\sigma }_{2}}}=\int _{-\infty }^{\infty }{\rm exp} (\eta \imath \sqrt{2}{{k}_{{\rm so}}}x){{[{{\phi }_{n^{\prime} }}(x)]}^{*}}{{\phi }_{q}}(x){\rm d}x\end{eqnarray} \tag{ 12 }$$

and

$$\begin{eqnarray}&&J_{n^{\prime} n}^{{{\sigma }_{1}}^{\prime} {{\sigma }_{2}}^{\prime} ,{{\sigma }_{1}}{{\sigma }_{2}}}=\int _{-\infty }^{\infty }{\rm exp} (\xi \imath \sqrt{2}{{k}_{{\rm so}}}X){{[{{\Phi }_{n^{\prime} }}(X)]}^{*}}{{\Phi }_{n}}(X){\rm d}X.\end{eqnarray} \tag{ 13 }$$

The sign of the exponent is determined by the pseudo-spin combinations: $(\eta ,\xi )=(-,+)$, $(+,+)$, $(-,-)$ and $(+,-)$ for $({{\sigma }_{1}}^{\prime} {{\sigma }_{2}}^{\prime} ,{{\sigma }_{1}}{{\sigma }_{2}})=(\uparrow \uparrow ,\uparrow \downarrow )$, $(\uparrow \uparrow ,\downarrow \uparrow )$, $(\downarrow \downarrow ,\uparrow \downarrow )$ and $(\downarrow \downarrow ,\downarrow \uparrow )$, respectively.

The integral $J_{n^{\prime} n}^{{{\sigma }_{1}}^{\prime} {{\sigma }_{2}}^{\prime} ,{{\sigma }_{1}}{{\sigma }_{2}}}$ over the center of mass coordinate coincides with the Fourier transform of the product of two one-dimensional harmonic oscillator eigenstates. For $\xi =\pm$ $(n^{\prime} \leqslant n)$, we find [50]

$$\begin{eqnarray}\begin{array}{rcl} J_{n^{\prime} n}^{{{\sigma }_{1}}^{\prime} {{\sigma }_{2}}^{\prime} ,{{\sigma }_{1}}{{\sigma }_{2}}} & = & {{a}_{{\rm ho}}}\sqrt{\pi }N_{n^{\prime} }^{{\rm ni}}N_{n}^{{\rm ni}}n^{\prime} !{{2}^{(n+{{n}^{\prime }})/2}} \\ {} & {} & \times L_{n^{\prime} }^{(n-n^{\prime} )}({{\left( {{a}_{{\rm ho}}}{{k}_{{\rm so}}} \right)}^{2}}){{\left( \pm \imath {{k}_{{\rm so}}}{{a}_{{\rm ho}}} \right)}^{n-n^{\prime} }} \\ {} & {} & \times {\rm exp} [-{{\left( {{k}_{{\rm so}}}{{a}_{{\rm ho}}} \right)}^{2}}/2], \\ \end{array}\end{eqnarray} \tag{ 14 }$$

where $L_{n}^{(n-n^{\prime} )}$ denotes the associated Laguerre polynomial.

The integral $I_{q^{\prime} q}^{{{\sigma }_{1}}^{\prime} {{\sigma }_{2}}^{\prime} ,{{\sigma }_{1}}{{\sigma }_{2}}}$ over the relative coordinate can be performed by expanding ${{[{{\phi }_{q^{\prime} }}(x)]}^{*}}$ and ${{\phi }_{q}}(x)$ in terms of the non-interacting harmonic oscillator functions $\phi _{l}^{{\rm ni}}(x)$, ${{\phi }_{q}}(x)={{{\rm lim} }_{{{l}_{{\rm max} }}\to \infty }}\sum _{l=0}^{{{l}_{{\rm max} }}}c_{l}^{(q)}\phi _{l}^{{\rm ni}}(x)$, where the expansion coefficients $c_{l}^{(q)}$ can be obtained analytically [8]. The integral $I_{q^{\prime} q}^{{{\sigma }_{1}}^{\prime} {{\sigma }_{2}}^{\prime} ,{{\sigma }_{1}}{{\sigma }_{2}}}$ then becomes a double sum over integrals that have the same structure as the center of mass integrals $J_{n^{\prime} n}^{{{\sigma }_{1}}^{\prime} {{\sigma }_{2}}^{\prime} ,{{\sigma }_{1}}{{\sigma }_{2}}}$. In the calculations reported below, we use a finite cutoff ${{l}_{{\rm max} }}$. The 'optimal' cutoff depends on the value of g considered, the number of relative functions ${{\phi }_{q}}(x)$ included in the basis and the desired accuracy. For $|g|=\infty$, we find, as in the g = 0 case, a closed analytical expression for the integral $I_{q^{\prime} q}^{{{\sigma }_{1}}^{\prime} {{\sigma }_{2}}^{\prime} ,{{\sigma }_{1}}{{\sigma }_{2}}}$. Having analytical expressions for the matrix elements of $\tilde{H}$, 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 equations (6)–(9), i.e., we work with the basis functions ${{\tilde{\psi }}_{X,\pm }}={{(\tilde{\psi }_{nq}^{\uparrow \uparrow }|\uparrow \rangle }_{1}}|\uparrow {{\rangle }_{2}}\pm \tilde{\psi }_{nq}^{\downarrow \downarrow }|\downarrow {{\rangle }_{1}}|\downarrow {{\rangle }_{2}})/\sqrt{2}$ and ${{\tilde{\psi }}_{x,\pm }}={{(\tilde{\psi }_{nq}^{\uparrow \downarrow }|\uparrow \rangle }_{1}}|\downarrow {{\rangle }_{2}}\pm \tilde{\psi }_{nq}^{\downarrow \uparrow }|\downarrow {{\rangle }_{1}}|\uparrow {{\rangle }_{2}})/\sqrt{2}$. By properly combining the parts of $\psi _{nq}^{{\rm sr}}$ 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₁₂ and Y₁₂. The operator P₁₂ exchanges the coordinates (position and spin) of particles 1 and 2. Basis functions that are unchanged under the operation P₁₂ are needed to describe states with bosonic symmetry $({{p}_{12}}=+1)$ and those that pick up a minus sign under the operation P₁₂ are needed to describe states with fermionic symmetry $({{p}_{12}}=-1)$. The operator Y₁₂ can be written as $\sigma _{x}^{(1)}\sigma _{x}^{(2)}P{{P}_{12}}$, where the parity operator P changes x_j to $-{{x}_{j}}$ (j = 1 and 2). The Y₁₂ operator determines the 'helicity' of the system. We label the eigenstates by $({{p}_{12}},{{y}_{12}})$, where ${{p}_{12}}=\pm 1$ and ${{y}_{12}}=\pm 1$ are defined by their actions on an eigenstate. The basis functions with $(+1,+1)$ symmetry, for example, are given by ${{\tilde{\psi }}_{X,+}}$ with ${{\phi }_{q}}(x)$ even and ${{\Phi }_{n}}(X)$ even, by ${{\tilde{\psi }}_{X,-}}$ with ${{\phi }_{q}}(x)$ even and ${{\Phi }_{n}}(X)$ odd, by ${{\tilde{\psi }}_{x,+}}$ with ${{\phi }_{q}}(x)$ even and ${{\Phi }_{n}}(X)$ even, and by ${{\tilde{\psi }}_{x,-}}$ with ${{\phi }_{q}}(x)$ odd and ${{\Phi }_{n}}(X)$ even.

As an example, figures 2(a)–(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\approx 0.1414{{a}_{{\rm ho}}}\hbar \omega$, and three different spin–orbit coupling strengths ${{k}_{{\rm so}}}$. The energy spectrum in figure 2(a) is, to leading order, given by the spectrum for $\delta ={{k}_{{\rm so}}}=g=0$. In this limiting case, the energies are equal to $(2j+1)\hbar \omega \pm \Omega$ and $2j\hbar \omega$, where $j=0,1,\ldots$. Finite g and ${{k}_{{\rm 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 $\frac{1}{\sqrt{2\pi }}\frac{(2q)!}{{{\left( q! \right)}^{2}}{{4}^{q}}}g/{{a}_{{\rm ho}}}$. The spin–orbit coupling term introduces, in the small Ω regime, a small down shift that is proportional to $k_{{\rm so}}^{2}$. This down shift is negligible in figure 2(a) but clearly visible in figures 2(b) and (c). Moreover, the spin–orbit coupling introduces avoided crossings. The broadest avoided crossings occur around $\Omega =\hbar \omega$, 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.

Zoom In Zoom Out Reset image size

Figures 2(d)–(f) show energy spectra for the strong coupling limit, i.e., for $|g|\to \infty$, as a function of the Raman coupling strength Ω for vanishing detuning δ. To facilitate the comparison between the large and small g limits, the spin–orbit coupling strengths ${{k}_{{\rm so}}}$ in figures 2(d)–(f) are the same as in figures 2(a)–(c). In the regime where $\Omega \ll {{\hbar }^{2}}k_{{\rm so}}^{2}/m$, the energies change approximately linearly with Ω (with positive, vanishing or negative slope). When $\Omega \gg {{\hbar }^{2}}k_{{\rm so}}^{2}/m$, the low-lying portion of the energy spectrum consists of approximately parallel energy levels that can be parameterized as $c-\Omega$, where c is a constant.

As already aluded to in the introduction, the theoretical framework developed can be generalized to higher-dimensional trapping geometries. For a spherically symmetric 3D system, e.g., the eigenstates of the 3D Hamiltonian ${{H}_{{\rm 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 higher-dimensional 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 two-particle system in 1D with contact interaction and spin–orbit 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.

Acknowledgements