Spin-orbit coupled bosons interacting in a two-dimensional harmonic trap

TL;DR

This paper investigates the ground state and low-energy spectra of two bosons with Rashba spin-orbit coupling in a 2D harmonic trap, revealing energy-level crossings and structural changes due to interactions.

Contribution

It provides an analytic and numerical analysis of the two-boson system, including beyond-approximation results and insights into energy-level crossings caused by spin-orbit coupling and interactions.

Findings

Energy-level crossing between states of different $J_z$ in the ground state.

Discontinuities in energy contributions and spatial density at strong interactions.

Relation between spin-orbit expectation value and $raket{\sigma_z raket{L_z}$ for eigenstates.

Abstract

A system of bosons in a two-dimensional harmonic trap in the presence of Rashba-type spin-orbit coupling is investigated. An analytic treatment of the ground state of a single atom in the weak-coupling regime is presented and used as a basis for a perturbation theory in the interacting two-boson system. The numerical diagonalization of both the single-particle and the two-boson Hamiltonian matrices allows us to go beyond those approximations and obtain not only the ground state, but also the low-energy spectra and the different energy contributions separately. We show that the expectation value of the spin-orbit term is related to the expectation value of $\hat{\sigma}_z \hat{L}_z$ for the eigenstates of the system, regardless of the trapping potential. The low-energy states of the repulsively interacting two-boson system are characterized. With the presence of a sufficiently strong interaction and spin-orbit coupling strength, there is a direct energy-level crossing in the ground state of the system between states of different $J_z$, the third component of the total angular momentum, that changes its structure. This is reflected in a discontinuity in the different energy terms and it is signaled in the spatial density of the system.

Equations

$$\begin{array}{rl}{\hat{H}}_{\mathrm{sp}}& =\frac{1}{2}m{\omega }^2\left({\hat{x}}^2+{\hat{y}}^2\right)+\frac{{\hat{p}}_x^2+{\hat{p}}_y^2}{2m}+\frac{{\kappa }^2}{2m}\\ & +\frac{\kappa }{m}\left({\hat{\sigma }}_x{\hat{p}}_x+{\hat{\sigma }}_y{\hat{p}}_y\right),\end{array}$$

$$\begin{array}{rl}{\hat{H}}_{\mathrm{sp}}& =({\hat{n}}_x+{\hat{n}}_y+1)\\ & +\frac{i\kappa }{\sqrt{2}}\left({\hat{\sigma }}_x\left({\hat{a}}_x^{†}-{\hat{a}}_x\right)+{\hat{\sigma }}_y\left({\hat{a}}_y^{†}-{\hat{a}}_y\right)\right)+\frac{{\kappa }^2}{2}.\end{array}$$

$$⟨\alpha |{\hat{H}}_{\mathrm{sp}}|\beta ⟩={ϵ}_{\alpha ,\beta }+\frac{{\kappa }^2}{2}{\delta }_{\alpha ,\beta },$$

$${ϵ}_{\alpha ,\beta }$$

$$E_{0,d=6}^{\mathrm{sp}}$$

$$E_{0,d=12}^{\mathrm{sp}}$$

$$\begin{array}{rl}{E}_{0,\mathrm{pert}}^{\mathrm{sp}}& =1-\frac{1}{2}{\kappa }^2+\frac{1}{2}{\kappa }^4-\frac{2}{3}{\kappa }^6+\frac{79}{72}{\kappa }^8\\ & -\frac{274}{135}{\kappa }^{10}+\frac{130577}{32400}{\kappa }^{12}.\end{array}$$

$$E_{0,\kappa \gg 1}^{\mathrm{sp}}=\frac{1}{2}+\frac{1}{8{\kappa }^2}.$$

$$|{\psi }_{E,B}^{\mathrm{sp}}⟩$$

$$|{\psi }_{E,A}^{\mathrm{sp}}⟩$$

$$\hat{\mathbf{J}}=\hat{\mathbf{S}}+\hat{\mathbf{L}},$$

$$E_{\kappa \gg 1}^{\mathrm{sp}}=n_r+\frac{1}{2}+\frac{j_z^2}{2{\kappa }^2},$$

$$2⟨{\psi }_E^{\mathrm{sp}}|{\hat{V}}_{\mathrm{ho}}|{\psi }_E^{\mathrm{sp}}⟩-2⟨{\psi }_E^{\mathrm{sp}}|\hat{K}|{\psi }_E^{\mathrm{sp}}⟩-⟨{\psi }_E^{\mathrm{sp}}|{\hat{V}}_{\mathrm{so}}|{\psi }_E^{\mathrm{sp}}⟩=0.$$

$$⟨{\psi }_E^{\mathrm{sp}}|\hat{K}|{\psi }_E^{\mathrm{sp}}⟩=-\frac{3}{4}⟨{\psi }_E^{\mathrm{sp}}|{\hat{V}}_{\mathrm{so}}|{\psi }_E^{\mathrm{sp}}⟩+\frac{2E^{\mathrm{sp}}-{\kappa }^2}{4},$$

$$E^{\mathrm{sp}}$$

$$⟨{\psi }_E^{\mathrm{sp}}|{\hat{V}}_{\mathrm{so}}|{\psi }_E^{\mathrm{sp}}⟩=-2{\kappa }^2\left(1+⟨{\psi }_E^{\mathrm{sp}}|{\hat{L}}_z{\hat{\sigma }}_z|{\psi }_E^{\mathrm{sp}}⟩\right).$$

$$\hat{H}={\hat{H}}_0+{\hat{H}}_{\mathrm{int}}+\frac{N{\kappa }^2}{2}.$$

$${\hat{H}}_0={\hat{V}}_{\mathrm{ho}}^T+{\hat{K}}^T+{\hat{V}}_{\mathrm{so}}^T,$$

$$V(|{\vec{x}}_i-{\vec{x}}_j|)=\frac{g}{\pi s^2}{e}^{-\frac{|{\vec{x}}_i-{\vec{x}}_j{|}^2}{s^2}},$$

$$a_{2D}\approx \sqrt{2}s{e}^{-\frac{\gamma }{2}-\frac{2\pi }{g}},$$

$${\hat{H}}_{\mathrm{int}}={\hat{H}}_{\uparrow \uparrow }+{\hat{H}}_{\downarrow \downarrow }+{\hat{H}}_{\uparrow \downarrow },$$

$${\hat{H}}_{\uparrow \uparrow }$$

$${\hat{H}}_{\downarrow \downarrow }$$

$${\hat{H}}_{\uparrow \downarrow }$$

$${\hat{H}}_0=\sum _{i,j=1}^M{\hat{a}}_i^{†}{\hat{a}}_j{ϵ}_{i,j},$$

$${\hat{H}}_{\mathrm{int}}$$

$$|n_1,...,n_M⟩=\frac{({\hat{a}}_1^{†}{)}^{n_1}\dots ({\hat{a}}_M^{†}{)}^{n_M}}{\sqrt{n_1!...n_M!}}|\mathrm{vac}⟩.$$

$$\sum _{i=1}^M{n}_i{ϵ}_{i,i}⩽E^{\max }=N_E+2,$$

$$D(N_E)=\sum _{k=0}^{N_E}\left(3d_{N_E}^b+d_{N_E}^f\right),$$

$$M=(N_E+1)(N_E+2).$$

$${\hat{J}}_z^T={\hat{J}}_z^{(1)}+{\hat{J}}_z^{(2)},$$

$$|{\Psi }_{0,AA}⟩=|{\psi }_{0,A}^{\mathrm{sp}}⟩|{\psi }_{0,A}^{\mathrm{sp}}⟩,$$

Full text

I Introduction

Spin-orbit coupling in ultracold atoms Dalibard ; Zhang ; Galitski ; Goldman ; Zhai ; Machon has been an issue of great interest in the last years in the atomic physics community. Since the first experiment was carried out successfully Lin dressing the atoms with two Raman lasers, additional investigations have been performed. For example, studying temperature effects Ji or engineering the spin-orbit coupling in alternative ways: with a gradient magnetic field Luo ; and within optical lattices Wu ; Grusdt ; Yamamoto . Interesting phenomena have been observed in spin-orbit coupled systems, for instance, a negative effective mass Khamehchi .

In the absence of a confining potential, in a homogeneous system, the single-particle energy dispersion relation is simple and the Hamiltonian is solvable in momentum space in the presence of spin-orbit coupling. In that case, for the many-body system at zero temperature, two phases were predicted in Ref. Wang in a mean-field approximation: the plane wave phase and the standing wave phase. The transition from one phase to the other was characterized depending on the inter- and intra-spin interactions between the atoms. Further studies in exploring the phase diagram of spin-orbit coupled Bose-Einstein condensates have been done within a mean-field description Li2 , studying the stability of the system against quantum and thermal fluctuations Baym ; Ozawa2 ; Ozawa3 ; Ozawa4 ; Kawasaki .

In the presence of a confining harmonic trap, the situation is fairly different, due to the introduction of a new characteristic length and the fact that the momentum is no longer a good quantum number. At the single-particle level, even when the spin-orbit coupling is strong, the spectrum remains discrete forming a Landau-level-like structure Anderson ; Zhou1 ; Sinah ; Li ; Hu , which is altered when the trap is anisotropic Marchukov ; Marchukov2 . At the mean-field level, more phases, like a half-quantum vortex state, are found in the trapped system Sinah ; Hu ; Ramachandhran ; Zhou2 ; Zhou1 ; Li ; Armaitis .

The inclusion of interactions between the atoms adds an additional challenge, specially in the strongly interacting regime Yin ; Guan2 ; Schillaci ; Guan , where quantum correlations are expected to dominate the physics Ramachandhran2 . Then, methods that go beyond mean field are required Juan .

In this work we make use of analytical approaches and numerical diagonalization techniques in order to describe the trapped single-particle and two-boson systems in the presence of Rashba spin-orbit coupling.

In Sec. II, the ground state of the single-particle system and the first low-energy states are computed and analyzed. We relate the different energy contributions and also the expectation values of different kind of spin-orbit coupling terms applying the virial theorem. In Sec. III, the interacting two-boson system is studied. First, we give the second-quantized $N$-boson Hamiltonian and explain the methodology to diagonalize it for the $N=2$ case. In second place, we discuss the degeneracy breaking in the three-fold degenerate ground-state subspace. In Sec. IV, we analyze the combined effects of the spin-orbit coupling and a spin-independent repulsive interaction in the spectrum. In particular, we find a crossover in the ground state characterized by a discontinuity in the energy contributions as a function of the spin-orbit coupling strength and by a change in the density profile of the system. Finally, conclusions and summary are presented in Sec. V.

II The single-particle system

The physics of a particle of mass $m$ in a two-dimensional isotropic harmonic potential of frequency $\omega$ with Rashba type spin-orbit coupling is described by the Hamiltonian

[TABLE]

where $\kappa$ is the spin-orbit coupling strength and $\hat{\sigma}_{x}$ and $\hat{\sigma}_{y}$ are Pauli matrices. In the present paper, as we consider a bosonic system of ultracold atoms, the spin part does not refer to the intrinsic spin but to an internal degree of freedom or pseudospin, for instance, two hyperfine atomic states as in Ref. Lin . The Hamiltonian is composed by the kinetic energy, $\hat{K}=(\hat{p}^{2}_{x}+\hat{p}^{2}_{y})/(2m)$, the harmonic potential, $\hat{V}_{\rm ho}=(m/2)\omega^{2}\left(\hat{x}^{2}+\hat{y}^{2}\right)$, the spin-orbit coupling, $\hat{V}_{\rm so}=(\kappa/m)\left(\hat{\sigma}_{x}\hat{p}_{x}+\hat{\sigma}_{y}\hat{p}_{y}\right)$, and the constant term $\kappa^{2}/(2m)$. As mentioned in Ref. Schillaci , up to a pseudospin rotation, an alternative and equivalent form of the Rashba term would be $\propto\left(\hat{\sigma}_{x}\hat{p}_{y}-\hat{\sigma}_{y}\hat{p}_{x}\right)$.

From now on, we use harmonic oscillator units, i.e., the energy is measured in units of $\hbar\omega$ and the length in units of $x_{\rm ho}\equiv\sqrt{\hbar/(m\omega)}$. The Hamiltonian in Eq. (1) is written in terms of annihilation operators, $\hat{a}_{x}=(\hat{x}+i\hat{p}_{x})/\sqrt{2}$ and $\hat{a}_{y}=(\hat{y}+i\hat{p}_{y})/\sqrt{2}$, and the corresponding creation operators, $\hat{a}^{\dagger}_{x}$ and $\hat{a}^{\dagger}_{y}$, as

[TABLE]

These operators fulfill the commutation relations $[\hat{{a}}_{i},\hat{a}^{\dagger}_{j}]=\delta_{ij}$ and $[\hat{{a}}_{i},\hat{a}_{j}]=[\hat{{a}}^{\dagger}_{i},\hat{a}^{\dagger}_{j}]=0$, with $i,j=x,y$. We have used the number operators $\hat{n}_{x}=\hat{a}^{\dagger}_{x}\hat{a}_{x}$ and $\hat{n}_{y}=\hat{a}^{\dagger}_{y}\hat{a}_{y}$. Notice that $\kappa$ is not a dimensionless parameter in the original Hamiltonian, Eq. (1), and it is written in units of $\sqrt{\hbar m\omega}$ in Eq. (2).

The single-particle basis can be labeled as, $\{\ket{n_{x},n_{y},m_{s}}\}$, with $n_{x},n_{y}=0,1,2,\,...\,$, and $m_{s}=-1,1$, where $n_{x}$, $n_{y}$ and $m_{s}$ are eigenvalues of $\hat{n}_{x}$, $\hat{n}_{y}$ and $\hat{\sigma}_{z}$, respectively.

The matrix elements of the single-particle Hamiltonian written using the basis introduced above read

[TABLE]

with

[TABLE]

and $\ket{\alpha}\equiv\ket{n_{x}(\alpha),n_{y}(\alpha),m_{s}(\alpha)}$. The index $\alpha$ labels each state of the single-particle basis. The Hamiltonian matrix is fully diagonalized using the first $5112$ states in order of increasing energy $\epsilon_{\alpha,\alpha}$, which corresponds to $(n_{x}+n_{y})\leqslant 70$ and $m_{s}=-1,1$. With this truncated Hilbert space the energies obtained are upper bounds to the exact ones. The method is variational, since we diagonalize in a subspace of the full Hilbert space.

II.1 The single-particle ground state

In this section, we explore the transition from the weak spin-orbit coupling regime, $\kappa<1$, to the strong spin-orbit coupling one, $\kappa\gg 1$, at the single-particle level. Our direct diagonalization results are compared with previously derived analytical expressions valid for the $\kappa\gg 1$ limit in Ref. Sinah , with perturbation theory expressions, $\kappa\ll 1$, derived in Ref. Yin , and with our own truncated analytic predictions valid in the $\kappa\lesssim 1$ regime.

In Fig. 1 we report the single-particle ground-state energy as a function of $\kappa$. The ground state is in all cases two-fold degenerated. For $\kappa=0$, we recover the harmonic oscillator result, $E_{0}^{\rm sp}=1$. As $\kappa$ is increased, the ground-state energy decreases towards an almost constant value of $E_{0}^{\rm sp}\simeq 0.5$, which is already reached for $\kappa\simeq 3$.

For $\kappa<1$, we derive analytical approximate expressions for the ground state of the single-particle Hamiltonian and its energy. The variational method consists in truncating the Hilbert space to a small number of modes (see Appendix A.1 for details). Analytic expressions can be obtained truncating to six or twelve modes,

[TABLE]

The goodness of these expressions is shown in Fig. 1, comparing them with the direct diagonalization and also with the perturbative calculations performed in Ref. Yin , that we write in our units as:

[TABLE]

Eq. (6) is seen to provide the best approximation to the direct diagonalization results, providing an accurate description up to $\kappa=1$. The perturbative expression of Ref. Yin , Eq. (7), reproduces well the results up to $\kappa\simeq 0.7$ while the approximation with six modes already fails for $\kappa\simeq 0.5$. The truncated analytical expressions fail to describe the ground state when it has relevant contributions from basis states that are not in the truncated subspace considered.

The large $\kappa$ domain has been studied previously in Refs. Anderson ; Zhou1 ; Sinah ; Li ; Hu . In this regime, approximate expressions for the two-degenerate states that define the ground-state subspace are given in Ref. Sinah , together with an expression for the ground-state energy,

[TABLE]

This approximation is in very good agreement with our numerical results for $\kappa>2$ (see Fig. 1). In particular, they correctly capture the limiting value in the spin-orbit dominated regime, $E_{0}^{\rm sp}\to 1/2$.

II.2 The single-particle energy spectrum

One of the important advantages of direct diagonalization methods is that they also provide, besides the ground-state properties, the low-energy part of the spectrum. The low-energy spectrum of the single-particle Hamiltonian, Eq. (2), is depicted in Fig. 2.

In the limiting case of $\kappa=0$, the eigenstates of the Hamiltonian are the eigenstates of two independent two-dimensional harmonic oscillators, one for each spin component. Therefore, the energies are $E^{\rm sp}_{n}=n+1$ with degeneracy $2(n+1)$ and $n=n_{x}+n_{y}$. The case of $\kappa<1$ is analyzed in Ref. Yin , where the exact numerical values are compared with perturbation theory calculations in $\kappa$.

For any value of $\kappa$, all energy levels are two-fold Kramers-degenerate because the Hamiltonian is time-reversal symmetric Ramachandhran ; Marchukov ; Hu ; Zhou1 ; Li . This degeneracy can be broken introducing a Zeeman term Marchukov . The effect of deforming the trap was considered in Ref. Marchukov2 , which results in a breaking of the cylindrical symmetry of the system. In our case, the time-reversal symmetry is preserved and, in order to distinguish between the pair of degenerate states, we label them with $A$ and $B$, respectively, for a given energy $E^{\rm sp}$. The action of the time reversal operator, $\hat{T}=i\hat{\sigma}_{y}\mathcal{C}$ Ramachandhran ; Marchukov ; Hu ; Li , with $\mathcal{C}$ the complex conjugation operator, on the two-fold degenerate eigenstates reads

[TABLE]

The eigenstates of the single-particle Hamiltonian can be written in a basis with a well defined total angular momentum,

[TABLE]

where $\hat{\bm{S}}=(\hat{\sigma}_{x},\hat{\sigma}_{y},\hat{\sigma}_{z})/2$, and $\hat{\bm{L}}\equiv\hat{\bm{r}}\times\hat{\bm{p}}=(0,0,\hat{L}_{z})$, both written in units of $\hbar$. The single-particle Hamiltonian commutes with $\hat{J}^{2}$ and $\hat{J}_{z}$. Therefore, the eigenstates of the system can be labeled with the corresponding quantum numbers, $j$ and $j_{z}$, respectively, regardless of the value of $\kappa$. In particular, in the limiting case $\kappa\gg 1$, an additional radial quantum number, $n_{r}$, is introduced to describe the eigenstates of the system (see Ref. Sinah ) and also the eigenenergies, approximately,

[TABLE]

with $n_{r}=0,1,\,...\,$, and $j_{z}=m_{l}+1/2$, with $m_{l}=0,\pm 1,\,...\,$. The two-fold degeneracy is reflected in the fact that the energy depends on $j^{2}_{z}$, so it is independent of its sign. The eigenstates with the same radial quantum number, $n_{r}$, tend to become degenerate with increasing $\kappa$, forming an energy manifold. This kind of physics has been studied in two and three dimensions, where the same type of Landau-level-like spectrum is found and described in terms of dimensional reduction Anderson ; Zhou1 ; Sinah ; Li ; Hu .

The approximate expression, Eq. (11), works very well for $\kappa\gg 1$, as seen in Fig. 3. For a given value of $\kappa$, the lowest eigenenergies are well-described and, as expected, the larger is the value of $\kappa$ the better is the approximation for a larger number of energy levels.

II.3 Energy contributions

As seen above, with increasing $\kappa$ the system goes from a harmonic oscillator behavior to a spin-orbit dominated one. The spectral properties are very different in both limits and feature a particularly involved structure in the intermediate region. To better understand the spin-orbit effects, we consider now the different energy contributions to the total energy of the different eigenstates as we vary the value of $\kappa$.

In Fig. 4, we show, for the first eigenstates of the single-particle system, how the total energy is distributed between the different energy contributions. As can be seen, the degeneracy due to the time-reversal symmetry of the system, that makes all eigenstates two-fold degenerate, is also reflected in the energy contributions. Each pair of degenerate states has also the same kinetic, harmonic potential, and spin-orbit coupling energies.

In the $\kappa=0$ limit, the eigenstates obey the equipartition relation valid for the harmonic oscillator, $\langle K\rangle=\langle V_{\rm ho}\rangle$ [see Fig. 4 panel (a)]. For a sufficiently small value of the spin-orbit coupling strength, those two contributions are not equal but of the same order of magnitude [see panels (b) and (c) of Fig. 4 for the cases $\kappa=0.5$ and $\kappa=1$, respectively]. Further increasing the value of $\kappa$, the situation changes, and the largest contributions, in absolute value, to the total energy are clearly the spin-orbit and kinetic parts [see Fig. 4 panel (d)]. There are, however, large cancellations between these two contributions which result in a total energy comparable to the harmonic oscillator part. Further insights into this energy decomposition and a nontrivial test to our numerical method is provided by the virial theorem (see Appendix B),

[TABLE]

For all the states considered, we have checked that the virial theorem energy relation is fulfilled, i.e., the left part of Eq. (12) represents less than $1\%$ of $E^{\rm sp}$. Actually, the cancellation needed comes from $\langle K\rangle$ and $\langle V_{\rm ho}\rangle$ for $\kappa=0$ and from $\langle K\rangle$ and $\langle V_{\rm so}\rangle$ in the large $\kappa$ domain.

In the presence of spin-orbit coupling, in panels (b)-(d) in Fig. 4, we observe a negative correlation between the spin-orbit coupling term and the kinetic energy. The relation between these two contributions is the following:

[TABLE]

and it arises from the virial theorem, Eq. (12), and from writing the energy as:

[TABLE]

II.4 Expectation value of the spin-orbit potential

The term that commonly appears in atomic and nuclear physics as spin-orbit coupling is proportional to $\hat{L}_{z}\hat{\sigma}_{z}$. The main difference between that kind of term and the Rasbha spin-orbit is that in one case the spin is coupled to the angular momentum and in the other to the linear momentum. However, we can relate the expectation values of both types of spin-orbit coupling terms,

[TABLE]

The eigenstates of the single-particle system obtained by exact diagonalization, whose energies are shown in Fig. 4, fulfill the previous relation, within a numerical error of less than a $1\%$ in the difference between both sides of Eq. (15).

The relation between the expectation values of the two kinds of spin-orbit terms is not a particularity of the pure Rashba case, it also works in a more general case, i.e. a mixture of Rashba and Dresselhaus spin-orbit couplings. Moreover, this property does not depend on the external trapping potential. The derivation of the relation in Eq. (15) is written in Appendix B, where we also generalize it and demonstrate its independence of the external trap.

III The two-boson system

In this section, we turn to the interacting few-body case. We first present our formalism which is developed for the general case of $N$ interacting bosons. Afterwards we specialize for the two-boson case.

Let us thus start with a system of $N$ interacting identical bosons trapped by an isotropic harmonic potential with Rashba spin-orbit coupling. The $N$-boson Hamiltonian reads

[TABLE]

The first part contains the total harmonic potential energy, $\hat{V}_{\rm ho}^{T}$, kinetic energy, $\hat{K}^{T}$, and spin-orbit energy, $\hat{V}_{\rm so}^{T}$,

[TABLE]

with $\hat{V}_{\rm ho}^{T}=(1/2)\sum_{i=1}^{N}\hat{\bm{x}}^{2}_{i}$, $\hat{K}^{T}=(1/2)\sum_{i=1}^{N}\hat{\bm{p}}^{\,2}_{i}$, and $\hat{V}_{\rm so}^{T}=\kappa\sum_{i=1}^{N}\left(\hat{\sigma}_{x_{i}}\hat{p}_{x_{i}}+\hat{\sigma}_{y_{i}}\hat{p}_{y_{i}}\right)$.

We model the atom-atom interaction with a Gaussian potential Mujal ,

[TABLE]

characterized by a finite range $s$ independent of the spin state and an interaction strength, $g$, which can be dependent on the spin Yin . The two interaction parameters are related to the two-dimensional scattering length, $a_{\rm{2D}}$, in Ref. Doganov by comparing the Gaussian potential to zero-range results. The approximate analytical expression that is obtained reads

[TABLE]

where $g$ is written in units of $\hbar^{2}/m$ and $\gamma$ is the Euler-Mascheroni constant. Beyond that approximation, in Ref. Jeszenszki , the authors combine an analytical treatment with numerical calculations to relate the interaction parameters and $a_{\rm{2D}}$.

The interaction part is divided in three contributions,

[TABLE]

where,

[TABLE]

For simplicity, we have introduced the following notation for the spin variable: $\ket{\uparrow}\equiv\ket{m_{s}=1}$, and $\ket{\downarrow}\equiv\ket{m_{s}=-1}$.

III.1 Second-quantized two-boson Hamiltonian

Despite the fact that our approach is in principle valid for a few number of bosons, we concentrate from now on in the two-boson case. The two-boson system provides a nontrivial example where the interplay of interactions and spin-orbit coupling can be studied in detail.

In our approach we solve numerically the time-independent Schrödinger equation for the two-boson Hamiltonian truncating the Hilbert space. We first consider that the particles can populate the first $M$ eigenstates of the harmonic trap, including the spin degree of freedom. In this case, we introduce the creation and annihilation operators, $\hat{a}^{\dagger}_{i}$ and $\hat{a}_{i}$, that create or annihilate bosons in the single-particle state $i=1,\,...\,,M$, respectively. They fulfill the commutation relations $[\hat{a}_{i},\hat{a}^{\dagger}_{j}]=\delta_{i,j}$ and $[\hat{a}_{i},\hat{a}_{j}]=[\hat{a}^{\dagger}_{i},\hat{a}^{\dagger}_{j}]=0$. The index $i$ labels the trio of quantum numbers $n_{x}$, $n_{y}$ and $m_{s}$, and increases with increasing the energy of the harmonic oscillator eigenstate $i$, $\epsilon_{i,i}=n_{x}(i)+n_{y}(i)+1$.

The second-quantized version of the single-particle part of Eq. (16) is

[TABLE]

where the explicit form of $\epsilon_{i,j}$ is given in Eq. (4). The interaction term is written as:

[TABLE]

where $V_{i,j,k,l}$ are computed analytically from the expressions given in Appendix C of Ref. Mujal , being aware that in the present article the indices $i$, $j$, $k$, and $l$ label the single-particle states in a different way and that the integrals depend on the quantum numbers $n_{x}$ and $n_{y}$ corresponding to the previous indices.

The Fock states are built creating particles into the vacuum state, $\ket{\rm vac}\equiv\ket{0,\,...\,,0}$, as follows:

[TABLE]

In the present work, we study the two-boson case, i.e., $N=\sum_{i=1}^{M}n_{i}=2\,.$ The basis that we use is the one formed by all the two-boson Fock states with

[TABLE]

where $N_{E}$ is a non-negative integer number. We truncate the Hilbert-space using this energy criterion Plodzien . In that case, the Hilbert space dimension considered is given by:

[TABLE]

where $d_{N_{E}}^{b}$ and $d_{N_{E}}^{f}$ are the number of spatially symmetric and antisymmetric degenerate two-particle states in a two-dimensional harmonic trap, given in Eqs. (21) and (22) of Ref. Mujal , respectively, and the factors $3$ and $1$ account for the triplet and singlet states of the spin part. The number of modes required to accomplish the energy truncation criterion in Eq. (25) is directly related to $N_{E}$,

[TABLE]

The low-energy eigenstates and eigenenergies of the two-boson Hamiltonian matrix are computed numerically using the ARPACK library. In the following section, we use a Hilbert space of dimension $D=17765$ corresponding to $M=420$ single-particle basis states. In Sec. IV, we need a larger Hilbert space, with $M=812$ and $D=63035$ [see Eqs. (26) and (27)].

III.2 Ground-state energy and degeneracy

In this section, we compute the ground-state energy, concentrating in understanding the way the interaction lifts the degeneracy of the ground-state manifold. To this aim, we compare our direct diagonalization results with approximate expressions for the energy of the ground-state manifold. In all cases discussed below, we set the spin-orbit coupling to a non-zero but small value, $\kappa=0.3$. Larger values of $\kappa$ are discussed in Sec. IV.

In absence of interactions, the ground state is three-fold degenerated. We obtain approximate analytic expressions for the energies of the three states using the six-mode truncation presented in Sec. II.1. The energies of the three states are denoted, $E_{0}^{AA}$, $E_{0}^{AB}$ and $E_{0}^{BB}$. Their explicit expressions are provided in Appendix A.2.

The simplest case we consider is when $g_{\uparrow\uparrow}=g_{\downarrow\downarrow}=g_{\uparrow\downarrow}=g$. In this case, the three orthogonal states that define the ground-state subspace remain quasidegenerate [see Fig. 5 panel (a)]. As we consider a small finite range, $s=0.5$, the $AB$ state, approximated by Eq. (45) at $g\approx 0$, has a slightly different energy within our approximation, and would be truly degenerate with the other two in the limit of $s\rightarrow 0$. The three-fold degeneracy of the ground-state manifold is lifted whenever the interaction strengths are not equal. For instance, fixing $g_{\downarrow\downarrow}=g_{\uparrow\downarrow}=0$, and increasing $g_{\uparrow\uparrow}$ we completely break the degeneracy, since the spin-orbit part of the Hamiltonian induces a nonzero, but different, spin-up spin-up component in all three orthogonal two-boson states. Our perturbative calculations are used to identify which energy level corresponds to each kind of state, as we show in Fig. 5 panel (b). For the case of the state of kind $AA$, the one with a larger spin-up spin-up component, we observe that the prediction of perturbation theory fails for $g_{\uparrow\uparrow}>1$. In contrast, for the state of kind $BB$, with a small spin-up spin-up component, its energy is well-approximated perturbatively up to $g_{\uparrow\uparrow}=20$.

The ground state remains degenerate, although only two-fold, if we set to zero the intraspin interactions, $g_{\downarrow\downarrow}=g_{\uparrow\uparrow}=0$, and vary the inter-spin one, $g_{\uparrow\downarrow}$. Since the effect on the states of kind $AA$ and $BB$ is the same, they remain degenerate and define the ground-state subspace [see Fig. 5 panel (c)]. However, the state $AB$ is very sensitive to changes in $g_{\uparrow\downarrow}$, compared to the two previous ones, and its energy increases more rapidly.

The last case we consider is fixing at finite values two of the interaction strengths, e.g. $g_{\downarrow\downarrow}$ and $g_{\uparrow\downarrow}$, and varying the other one, $g_{\uparrow\uparrow}$ [see Fig. 5 panel (d)]. In this case, we find crossings between the energy levels. The perturbative calculations are useful to predict the value of $g_{\uparrow\uparrow}$ where the crossing occurs, by equating Eqs. (46), (47) and (48), properly, once $g_{\downarrow\downarrow}$ and $g_{\uparrow\downarrow}$ are fixed. In particular, in Fig. 5 panel (d), we see that it happens when $g_{\uparrow\uparrow}=g_{\downarrow\downarrow}$, and also when $g_{\uparrow\uparrow}=g_{\uparrow\downarrow}$.

Finally, we observe that when we further increase the interaction strength, regardless of the spin components, the energy levels tend to saturate. This behavior is not captured by the perturbative expressions discussed. This is an indicator that the system becomes correlated in the proper way in order to reduce the total energy by avoiding the atom-atom interaction. This kind of behavior was found previously in a harmonically trapped system of interacting bosons in two dimensions Mujal ; Mujal2 .

IV Interaction induced crossover in the $g_{\uparrow\uparrow}=g_{\downarrow\downarrow}=g_{\uparrow\downarrow}$ case

Now let us broaden our scope and study not only the ground-state manifold but also the lower part of the energy spectrum. The goal is to discuss the combined effects of the spin-orbit term and the atom-atom interaction. For simplicity, we consider the case $g=g_{\uparrow\uparrow}=g_{\downarrow\downarrow}=g_{\uparrow\downarrow}$, with $g\geqslant 0$.

In this situation, the third component of the total angular momentum of the two-particle system,

[TABLE]

commutes with the Hamiltonian in Eq. (16), so the corresponding quantum number, $J_{z}$, is a good quantum number to label the eigenstates of $\hat{H}$. Moreover, due to the time-reversal symmetry of the Hamiltonian, the states with $J_{z}\neq 0$ are at least two-fold degenerate, i.e., with $\pm J_{z}$. The value of $J_{z}$ is obtained from the single-particle spectrum in the noninteracting two-particle case and by the numerical diagonalization of $\hat{J}_{z}^{T}$ in the degenerate subspaces corresponding to each eigenenergy.

The interaction has three main effects, as seen in Fig. 6, where we compare the low energy spectrum for $g=0$, panels (a) and (b), with the corresponding ones for $g=3$, panel (c) and (d). In Fig. 6 panels (a) and (c), we vary $\kappa\in[0,1]$, while in panels (b) and (d) we consider a wider region, i.e., $\kappa\in[0,3]$. Due to the repulsive character of the interaction, the energies are shifted to higher values, see for instance the case of the three-fold degenerate ground-state energy level. A second effect, is the breaking of degeneracies. For instance, already at $\kappa=0$, the first excited state, with degeneracy 8, breaks in two levels with degeneracy 2 for the lowest level and 6 for the highest one. These degeneracies are further broken when increasing $\kappa$ [see Fig.  6 panels (a) and (c)]. Finally, the breaking of degeneracies is accompanied by the presence of more energy-level crossings, which are avoided between energy levels with the same value of $J_{z}$ and direct otherwise.

As seen in Fig. 6 panel (d), we find a crossing at the ground-state level which appears at $\kappa\approx 2.65$ for $g=3$. The numerical calculation of $J_{z}$ indicates a direct crossing from a nondegenerate ground state with $J_{z}=0$ to a two-fold degenerate ground-state subspace with two states that can be labeled with $J_{z}=-3$ and $J_{z}=3$. Regarding the first-excited energy level, which is two-fold degenerate with $J_{z}=-1$ and $J_{z}=1$, it presents an avoided crossing at $\kappa\approx 2.25$.

As we will discuss later in Fig. 7, we find a discontinuity in the different energy contributions to the ground-state energy as expected for a direct crossing. In the following paragraphs, we concentrate in characterizing this level crossing which corresponds to a change in structure of the ground state induced by the spin-orbit term in the presence of interactions.

Starting from $\kappa=0$ and $g=0$, panel (b) of Fig. 6, the ground state is three-fold degenerate. In this case, one could use as a basis of that subspace the two-boson states formed by putting the two bosons in the ground state of the two-dimensional harmonic trap with parallel spins, both pointing up or both pointing down, and with anti-parallel spins.

For $\kappa>0$ the previous three states are no longer eigenstates, since the spin-orbit imposes a different form for the eigenstates at the single-particle level, that was discussed in Sec. II.1. However, the ground-state degeneracy remains unchanged with increasing $\kappa$ in the noninteracting case. The three states that define the ground-state subspace are

[TABLE]

[TABLE]

and

[TABLE]

constructed with the two-degenerate single-particle eigenstates, $\ket{\psi^{\rm sp}_{0,A}}$ and $\ket{\psi^{\rm sp}_{0,B}}$, of the Hamiltonian in Eq. (1).

In the interacting case the three-fold degenerate ground-state subspace splits in two energy levels: the ground state becomes nondegenerate, with $J_{z}=0$, and the first excitation becomes two-fold degenerate, corresponding to one state with $J_{z}=-1$ and the other with $J_{z}=1$. This effect is more notorious for larger $\kappa$, for instance for $\kappa=1.5$ in Fig. 6 panel (d), where we observe the gap opening. For larger $\kappa$ we observe the previously mentioned crossing. From $\kappa\approx 2.65$ up to $3$, the ground state becomes two-fold degenerate with $J_{z}=-3$ and $J_{z}=3$, respectively. The level which crosses at $\kappa\approx 2.65$ corresponds to the evolution with $\kappa$ of an excited level with $E=4$ at $\kappa=0$ that directly crosses multiple levels [see Fig. 6 panels (c) and (d)]. Let us emphasize that this transition is a joint effect of the spin-orbit coupling and the interaction, since it is only observed when both effects are present.

To characterize the crossing in the ground-state energy we have computed its energy contributions in the cases of Fig. 6 panels (b) and (d). These results are shown in Fig. 7, where we have also tested the fulfillment of the virial theorem energy relation (see Appendix B).

Before the crossing, the dependence on $\kappa$ of the kinetic, the harmonic potential and the spin-orbit energies is qualitatively similar to the noninteracting case (see Fig. 7). In the interacting case, the atoms are farther from the center of the trap resulting in a shift in the harmonic potential energy between the $g=0$ and $g=3$ cases depicted in Fig. 7. The kinetic energy is reduced in the interacting case. The interaction energy and the term coming from the interaction present in the virial relation, $\langle\hat{W}^{T}\rangle$, are mostly independent of $\kappa$. This fact explains that the correlations between the spin-orbit coupling energy and the kinetic energy, given in Eq. (B.1), are similar to the noninteracting case. At the crossing, except from the total energy that remains continuous, all other energy terms feature a discontinuity. After the crossing, the ground state has a different structure. The harmonic potential and the kinetic energies are larger than before. Again, this positive terms are compensated by the negative spin-orbit term that is larger in absolute value. In particular, the state is less sensitive to the presence of the repulsive interaction, since the interaction energy is smaller and closer to zero compared to the other energy terms.

This fact is better observed in Fig. 8, where the crossing is found varying $g$ with fixing $\kappa$. The energy level that finally becomes the ground state is one of the levels forming the third-excited manifold at $g=0$. This level directly crosses with multiple levels when $g$ is increased since its energy is less sensitive to the increase of the interaction strength value. In parallel, there are avoided crossings between the states that have the same quantum number $J_{z}$, as in the cases of $|J_{z}|=1$ and $|J_{z}|=2$, at $g\approx 2$ and $g\approx 3$, respectively.

The effects of the crossover become also apparent in the density of the cloud (see Appendix C for the explicit expressions). To illustrate this phenomenology we compare the densities for the $g=0$ and $g=3$ cases, for two values before and after the level crossing, $\kappa=2$ and $\kappa=3$, respectively. For $\kappa=2$ we observe that the total density of the cloud is similar in both cases (see Fig. 9). The main difference is that the interacting cloud is already larger than the noninteracting one, as expected from the repulsive nature of the atom-atom interactions considered. The densities of the two spin components are different for $g=0$ and $g=3$. In the interacting case, both densities are very similar, while in the noninteracting one $\rho^{\downarrow}$ is much smaller and more peaked at the center of the trap. An important effect of the crossing is that the cloud becomes larger after the level crossing, i.e. going from $\kappa=2$ to $\kappa=3$ for $g=3$ (see the total density in Fig. 9). This is in contrast with the behavior observed in absence of interactions, where the cloud size gets reduced when going from $\kappa=2$ to $\kappa=3$, as seen in Fig. 9. This effect is observed also for the densities of each component separately. Another relevant feature is that, after the crossing, the total density has a dip in the center of the trap, while in the noninteracting case it has a maximum.

V Summary and Conclusions

We have considered one and two bosons trapped in a harmonic potential with the presence of spin-orbit coupling. For the single-particle case, the diagonalization of the Hamiltonian matrix has allowed us to study the properties of the low-energy eigenstates of the system, going from the weak spin-orbit coupling regime to the strong one. We have computed the expectation values of each energy term in the Hamiltonian for the eigenstates, separately, and have derived and tested the virial energy relation between them. In particular, we have found a relation between the expectation value of different kind of spin-orbit coupling terms, which is independent of the trapping potential. For the ground state of the single-particle system, we have derived approximate analytical expressions that are able to reproduce the ground-state energy in the weak spin-orbit coupling regime and that, for the interacting two-boson system, are used to obtain perturbative expressions that explain the breaking of the degeneracy of the ground-state subspace when changing the values of the spin-dependent interaction strengths. In all cases, we have found that the ground-state energy tends to saturate with increasing the strength of the interaction, departing from the perturbation-theory prediction. This signals the formation of repulsive correlations in the system. In addition, in the spin-independent interaction case, for the repulsively interacting two-boson system, we have found a direct crossing in the ground state corresponding to states with different values of $J_{z}$, from $J_{z}=0$ to $J_{z}=\pm 3$. Some avoided crossings between excited states with the same $J_{z}$ are observed when the spin-orbit coupling parameter is sufficiently large. The change in the ground state has been characterized by computing the energy contributions, that present a discontinuity at the point where there is a direct energy-level crossing in the ground-state energy. Moreover, this phenomenon has been observed to be apparent in the density profile of the system, which could be experimentally measured Pyzh ; McDonald ; Subhankar .

Acknowledgements.

Useful discussion with Gordon Baym at the early stage of the work are gratefully acknowledged. We also thank Ferran Mazzanti and Juan Sánchez-Baena for their comments and for sharing their results with us. P.M. wants to thank Doerte Blume for giving him the opportunity to visit her group and thank her and Qingze Guan, Jianwen Jie and Jugal Talukdar for their warm hospitality and useful discussions. We acknowledge financial support from the Spanish Ministerio de Economia y Competitividad Grant No FIS2017-87534-P, from Generalitat de Catalunya Grant No. 2017SGR533 and from the European Union Regional Development Fund within the ERDF Operational Program of Catalunya (project QUASICAT/QuantumCat). P.M. was supported by a FI grant from Generalitat de Catalunya.

Appendix A Analytical approximations in the weak spin-orbit coupling regime

A.1 Single-particle case

In a first approximation, we consider a Hilbert space of dimension $6$, where the particle can populate the ground state of the harmonic oscillator or one of the two first-excited states of the trap, considering also the two possible spin orientations. Therefore, we consider the basis $\{\ket{n_{x},n_{y},m_{s}}\}=\{\ket{0,0,1},\ket{0,0,-1},\ket{1,0,1},\ket{1,0,-1},\ket{0,1,1},\ket{0,1,-1}\}$. In this Hilbert space, we construct the Hamiltonian matrix and diagonalize it analytically with Mathematica. In this way, we find approximate expressions for the ground state and its energy depending on the spin-orbit coupling strength, $\kappa$. The single-particle ground-state energy is approximately given by,

[TABLE]

The ground state is two-fold degenerate, and we label with $A$ and $B$ the orthogonal states,

[TABLE]

and

[TABLE]

where $C_{0}$ and $C_{1}$ are given by

[TABLE]

and

[TABLE]

Repeating the previous procedure with a Hilbert space of dimension $12$, we obtain more accurate expressions for the ground-state energy, given by,

[TABLE]

and also for the coefficients of the two degenerate states

[TABLE]

and

[TABLE]

where $D_{0}$, $D_{1}$ and $D_{2}$ are given by

[TABLE]

[TABLE]

and

[TABLE]

A.2 Two-boson case

Within the first single-particle approximation for small $\kappa$, discussed in Sec. II.1, we compute the energy of the following two-boson states:

[TABLE]

[TABLE]

and

[TABLE]

up to first order in perturbation theory for the interaction strength parameters $g_{\uparrow\uparrow}$, $g_{\downarrow\downarrow}$, and $g_{\uparrow\downarrow}$. The previous three states describe, approximately, the degenerate two-boson ground-state subspace in the noninteracting limit. The approximation becomes exact in the limit of $\kappa\rightarrow 0$. The first part of the energy for all of them is computed multiplying the single-particle energy given in Eq. (5) by the number of particles, that is $2$. The interaction part arises from computing the expectation values $\bra{\Phi_{0,AA}}\hat{H}_{\rm int}\ket{\Phi_{0,AA}}$, $\bra{\Phi_{0,BB}}\hat{H}_{\rm int}\ket{\Phi_{0,BB}}$, and $\bra{\Phi_{0,AB}}\hat{H}_{\rm int}\ket{\Phi_{0,AB}}$, since $\bra{\Phi_{0,AA}}\hat{H}_{\rm int}\ket{\Phi_{0,BB}}=\bra{\Phi_{0,AA}}\hat{H}_{\rm int}\ket{\Phi_{0,AB}}=\bra{\Phi_{0,BB}}\hat{H}_{\rm int}\ket{\Phi_{0,AB}}=0$. Therefore, the energies are

[TABLE]

[TABLE]

and

[TABLE]

where $C_{0}$ and $C_{1}$ depend on $\kappa$ and are given in Eq. (35) and Eq. (36) of Appendix A.1, respectively.

A particular limit case of interest is the short-range limit, $s\rightarrow 0$. In that case, the previous expressions reduce to

[TABLE]

[TABLE]

and

[TABLE]

Appendix B Virial relations

B.1 Virial theorem energy relation

For the eigenstates, $\ket{\Psi_{E}}$, of the Hamiltonian in Eq. (16), i.e., $\hat{H}\ket{\Psi_{E}}=E\ket{\Psi_{E}}$, the virial theorem establishes that

[TABLE]

with $\hat{\mathcal{O}}^{T}=\sum_{i=1}^{N}\left(\hat{x}_{i}\hat{p}_{x_{i}}+\hat{y_{i}}\hat{p}_{y_{i}}\right)$. The explicit computation of the expectation value of the commutator on the left part of the previous equation results in:

[TABLE]

where the last three terms come from the interaction part of the Hamiltonian (20) and the operators involved read:

[TABLE]

[TABLE]

and

[TABLE]

We also define the operator:

[TABLE]

Additional relations between the energy contributions are derived from the virial theorem in Eq. (53) and the fact that the total energy is given by:

[TABLE]

For instance, the generalization of Eq. (13) in the presence of interactions for the $N$-particle system, that relates the spin-orbit coupling term and the kinetic energy, reads

[TABLE]

In the noninteracting case, with the relation in Eq. (53) we can write the eigenenergies of the Hamiltonian in Eq. (16) as:

[TABLE]

In the single-particle case, the virial theorem energy relation, Eq. (53), reduces to Eq. (12).

B.2 Angular momenta and spin-orbit virial relation

Following the same procedure of previous Sec. B.1, we compute the expectation value of the following commutator:

[TABLE]

with $\hat{\mathcal{O}}^{T}=\sum_{i=1}^{N}\kappa\left(\hat{x}_{i}\hat{\sigma}_{x_{i}}+\eta\hat{y_{i}}\hat{\sigma}_{y_{i}}\right)$. In this case, we have used the general many-body Hamiltonian, that describes a noninteracting system,

[TABLE]

where the external trap is an arbitrary potential of the form

[TABLE]

and the spin-orbit term is a mixture of Rashba and Dresselhaus of the form:

[TABLE]

As a result, we find that

[TABLE]

where now, $\ket{\Psi_{E}}$ are the eigenstates of $\hat{H}^{RD}$. The independence of the external trapping potential arises from the fact that

[TABLE]

In the single-particle case and with a pure Rashba-type spin-orbit coupling the relation of Eq. (65) is equivalent to Eq. (15).

Appendix C Densities

The total density is computed as the expectation value of the operator

[TABLE]

which is decomposed as

[TABLE]

with

[TABLE]

and

[TABLE]