Nematic-orbit coupling and nematic density waves in spin-1 condensates

TL;DR

This paper introduces a method to create artificial nematic-orbit coupling in spin-1 Bose-Einstein condensates using a microwave chip, revealing new quantum phases and density modulations without the need for Raman lasers.

Contribution

It demonstrates a novel way to induce nematic-orbit coupling in spin-1 condensates via a designed microwave field, exploring resulting phases and excitations.

Findings

Identification of three nematic quantum phases.

Prediction of periodic nematic density modulations.

Low energy excitation spectra for each phase.

Abstract

We propose the creation of artificial nematic-orbit coupling in spin-1 Bose-Einstein condensates, in analogy to spin-orbit coupling. Using a suitably designed microwave chip, the quadratic Zeeman shift, normally uniform in space, can be made to be spatio-temporally varying, leading to a coupling between spatial and nematic degrees of freedom. A phase diagram is explored where three quantum phases with the nematic order emerge: easy-axis, easy-plane with single-well and easy-plane with double well structure in momentum space. By including spin-dependent and spin-independent interactions, we also obtain the low energy excitation spectra in these three phases. Lastly, we show that the nematic-orbit coupling leads to a periodic nematic density modulation in relation to the period $\lambda_T$ of the cosinusoidal quadratic Zeeman term. Our results point to the rich possibilities for manipulation of tensorial degrees of freedom in ultracold gases without requiring Raman lasers, and therefore, obviating light-scattering induced heating.

Equations

$${\hat{H}}_{{\rm IP}}=\int\mathrm{d}{\bf r}\sum_{a}\psi_{a}^{\dagger}({\bf r})\bigg{[}\frac{{\bf p}^{2}}{2m}{\hat{\bf 1}}+V({\bf r}){\hat{\bf 1}}+q({\bf r},t){\hat{\bf F}}_{z}^{2}\bigg{]}\psi_{a}({\bf r}),$$

$${\hat{H}}_{\mathrm{IP}}=\sum _{{\mathbf{k}}_{\perp }n}{\hat{\varphi }}_{{\mathbf{k}}_{\perp }n}^{†}{\mathbf{H}}_{\mathrm{D}}{\hat{\varphi }}_{{\mathbf{k}}_{\perp }n}+\left[\Omega {\hat{\varphi }}_{{\mathbf{k}}_{-},1}^{†}{\hat{\mathbf{F}}}_z^2{\hat{\varphi }}_{{\mathbf{k}}_{+},2}+H.c.\right].$$

$$E_{\alpha,\beta}({\bf k}_{\perp})=q+\frac{\hbar^{2}}{2m}\bigg{[}k_{\perp}^{2}+\frac{1}{4}k^{2}_{T}\bigg{]}\pm\sqrt{\bigg{[}\frac{\hbar^{2}}{2m}k_{x}k_{T}\bigg{]}^{2}+\Omega^{2}}.$$

$$\left(\begin{array}{c}{\chi }_{a\alpha }({\mathbf{k}}_{\perp })\\ {\chi }_{a\beta }({\mathbf{k}}_{\perp })\end{array}\right)=\left(\begin{array}{cc}{u}_{+\alpha }({\mathbf{k}}_{\perp })& u_{-\alpha }({\mathbf{k}}_{\perp })\\ u_{+\beta }({\mathbf{k}}_{\perp })& u_{-\beta }({\mathbf{k}}_{\perp })\end{array}\right)\left(\begin{array}{c}{\varphi }_{1,a}({\mathbf{k}}_{-})\\ {\varphi }_{2,a}({\mathbf{k}}_{+})\end{array}\right),$$

$${\tilde{H}}_0=\sum _{\begin{array}{c}{\mathbf{k}}_{\perp }{{\mathbf{k}}^{\prime }}_{\perp }{\mathbf{p}}_{\perp }\\ a{a}^{\prime }\{n_i\}\end{array}}{C}_{\{n_i\}}{\Lambda }_{n_1{n}_2}^{†a{a}^{\prime }}({\mathbf{k}}_{p-},{\mathbf{k}}_{p+}^{\prime }){\Lambda }_{n_3{n}_4}^{a^{\prime }a}({\mathbf{k}}_{p-}^{\prime },{\mathbf{k}}_{p+}),$$

$$\begin{array}{rl}& {\Lambda }_{n_1{n}_2}^{†a{a}^{\prime }}({\mathbf{k}}_{p-},{\mathbf{k}}_{p+}^{\prime })={\varphi }_{n_1,a}^{†}({\mathbf{k}}_{p-}){\varphi }_{n_2,a^{\prime }}^{†}({\mathbf{k}}_{p+}^{\prime }),\\ & {\Lambda }_{n_3{n}_4}^{a^{\prime }a}({\mathbf{k}}_{p-}^{\prime },{\mathbf{k}}_{p+})={\varphi }_{n_3,a^{\prime }}({\mathbf{k}}_{p-}^{\prime }){\varphi }_{n_4,a}({\mathbf{k}}_{p+}).\end{array}$$

$${\tilde{H}}_2=\sum _{\begin{array}{c}{\mathbf{k}}_{\perp }{\mathbf{k}}_{\perp }^{\prime }{\mathbf{p}}_{\perp }\\ a{a}^{\prime }b{b}^{\prime }\{n_i\}\end{array}}{C}_{\{n_i\}}{\hat{\mathbf{J}}}_{n_1{n}_4}^{ab}({\mathbf{k}}_{p-},{\mathbf{k}}_{p+})\cdot {\hat{\mathbf{J}}}_{n_2{n}_3}^{a^{\prime }{b}^{\prime }}({\mathbf{k}}_{p+}^{\prime },{\mathbf{k}}_{p-}^{\prime })$$

$$\begin{array}{rl}& {\hat{\mathbf{J}}}_{n_1{n}_4}^{ab}({\mathbf{k}}_{p-},{\mathbf{k}}_{p+})={\varphi }_{n_1,a}^{†}({\mathbf{k}}_{p-}){\hat{\mathbf{F}}}_{ab}{\varphi }_{n_4,b}({\mathbf{k}}_{p+})\\ & {\hat{\mathbf{J}}}_{n_2{n}_3}^{a^{\prime }{b}^{\prime }}({\mathbf{k}}_{p+}^{\prime },{\mathbf{k}}_{p-}^{\prime })={\varphi }_{n_2,a^{\prime }}^{†}({\mathbf{k}}_{p+}^{\prime }){\hat{\mathbf{F}}}_{a^{\prime }{b}^{\prime }}{\varphi }_{n_3,b^{\prime }}({\mathbf{k}}_{p-}^{\prime })\end{array}$$

$$\begin{gathered} \hat{H}=G_{\mathrm{sw}}+\frac{1}{2}\sum _{\mathbf{k}\ne 0}{\mathbf{X}}_{\mathbf{k}}^{†} \\ \left(\begin{array}{cc}{\mathbf{E}}_1& \mathbf{D}\\ {\mathbf{D}}^{†}& {\mathbf{E}}_{\bar{1}}\end{array}\right){\mathbf{X}}_{\mathbf{k}}. \end{gathered}$$

$${\mathbf{E}}_a=\left(\begin{array}{cc}{E}_g({\mathbf{k}}_{\perp })+c& f{e}^{i2{\Phi }_a}\\ f{e}^{-i2{\Phi }_a}& E_g({\mathbf{k}}_{\perp })+c\end{array}\right),$$

$$\mathbf{D}=\left(\begin{array}{cc}d{e}^{i({\Phi }_1-{\Phi }_{\bar{1}})}& g{e}^{i({\Phi }_{\bar{1}}+{\Phi }_1)}\\ g{e}^{-i({\Phi }_1+{\Phi }_{\bar{1}})}& d{e}^{-i({\Phi }_1-{\Phi }_{\bar{1}})}\end{array}\right),$$

$$\begin{array}{rl}& {\widetilde{ϵ}}_{b,1}({\mathbf{k}}_{\perp })=\sqrt{{\left[{\tilde{E}}_g({\mathbf{k}}_{\perp })+(\tilde{c}+\tilde{d})\right]}^2-(\tilde{f}+\tilde{g}{)}^2},\\ & {\widetilde{ϵ}}_{b,2}({\mathbf{k}}_{\perp })=\sqrt{{\left[{\tilde{E}}_g({\mathbf{k}}_{\perp })+(\tilde{c}-\tilde{d})\right]}^2-(\tilde{f}-\tilde{g}{)}^2},\end{array}$$

$$\hat{H}=G_{\mathrm{dw}}+\frac{1}{2}\sum _{\mathbf{k}\ne 0}{\mathbf{Y}}_{\mathbf{k}}^{†}\left(\begin{array}{cc}{\mathbf{M}}_{LL}& {\mathbf{M}}_{LR}\\ {\mathbf{M}}_{RL}& {\mathbf{M}}_{RR}\end{array}\right){\mathbf{Y}}_{\mathbf{k}},$$

$${\hat{H}}_{IP}^{\prime }=\int {\mathrm{d}}^2{r}_{\perp }\left(\begin{array}{cc}{\hat{\psi }}_1^{\ast }& {\hat{\psi }}_2^{\ast }\end{array}\right)\left(\begin{array}{cc}\frac{{\mathbf{p}}_{\perp }^2}{2m}+q{\hat{F}}_z^2& \Omega e^{-i{k}_{T}x}{\hat{F}}_z^2\\ \Omega e^{i{k}_{T}x}{\hat{F}}_z^2& \frac{{\mathbf{p}}_{\perp }^2}{2m}+q{\hat{F}}_z^2\end{array}\right)\left(\begin{array}{c}{\hat{\psi }}_1\\ {\hat{\psi }}_2\end{array}\right),$$

$${\hat{\mathcal{H}}}_{\mathrm{I}}=\frac{c_0}{2}{\left[|{\mathbf{\Psi }}_1(\mathbf{r}){|}^2+|{\mathbf{\Psi }}_{\bar{1}}(\mathbf{r}){|}^2\right]}^2+\frac{c_2}{2}{\left[|{\mathbf{\Psi }}_1(\mathbf{r}){|}^2-|{\mathbf{\Psi }}_{\bar{1}}(\mathbf{r}){|}^2\right]}^2,$$

$${\bf\Psi}_{a}({\bf r})={\cal A}_{\rm sw}e^{-i\frac{\vartheta}{2}}\bigg{[}e^{i\big{(}\frac{k_{T}}{2}x-\frac{\vartheta}{2}\big{)}}\varphi_{2}(z)-e^{-i\big{(}\frac{k_{T}}{2}x+\frac{\vartheta}{2}\big{)}}\varphi_{1}(z)\bigg{]},$$

$${\bf\Psi}_{a}({\bf r})={\cal A}^{\prime}_{\rm dw}\sum_{\begin{subarray}{c}j=\pm\\ l=\pm\end{subarray}}\bigg{[}u_{j\alpha}(l{\widetilde{k}}_{0})e^{i\big{[}(lk_{0}+j\frac{k_{T}}{2})x-j\frac{\vartheta}{2}+l\frac{\vartheta_{LR}}{2}\big{]}}\bigg{]}\varphi_{j}(z)$$

$${\hat{H}}_{\mathrm{A}}=A\hat{\mathbf{I}}\cdot \hat{\mathbf{S}}+\frac{g{\mu }_B}{\hslash }\hat{\mathbf{S}}\cdot \left[B_0\hat{\mathbf{z}}+{\mathbf{B}}_1(x,z,t)\right],$$

$$H_{\mathrm{A}}=\sum _{m_1}(\hslash {\omega }_L{m}_1+\frac{1}{2}\hslash {\Delta }_0)|1,m_1⟩⟨1,m_1|+\sum _{m_2}(\hslash {\omega }_L{m}_2-\frac{1}{2}\hslash {\Delta }_0)|2,m_2⟩⟨2,m_2|+\sum _{m_1,m_2}\left[\frac{1}{2}\hslash {\Omega }_{1,m_1}^{2,m_2}|2,m_2⟩⟨1,m_1|+h.c.\right]$$

$${\Omega }_{1,m_1}^{2,m_2}=\frac{2{\mu }_B}{\hslash }⟨2,m_2|{\mathbf{B}}_1\cdot \hat{\mathbf{S}}|1,m_1⟩$$

$${\Delta }_{1,m_1}^{2,m_2}={\Delta }_0+(m_1-m_2){\omega }_L.$$

$${\mathbf{B}}_1(x,z,t)=\sum _{l=1}^{\infty }{B}_l{e}^{-l\pi z/d}\left(\sin \left(\frac{l\pi x}{d}\right)\hat{\mathbf{x}}-\cos \left(\frac{l\pi x}{d}\right)\hat{\mathbf{z}}\right)\cos {\omega }_{c}t$$

$${\mathbf{B}}_1(x,z,t)\approx B_1{e}^{-\pi z/d}\left(\sin \left(\frac{\pi x}{d}\right)\hat{\mathbf{x}}-\cos \left(\frac{\pi x}{d}\right)\hat{\mathbf{z}}\right)\cos {\omega }_{c}t$$

$${\mathbf{B}}_1(x,z,t)=B_1{e}^{-\pi h/d}\left(1-\frac{\pi z^{\prime }}{d}\right)\left[\sin \left(\frac{\pi (x-vt)}{d}\right)\hat{\mathbf{x}}-\cos \left(\frac{\pi (x-vt)}{d}\right)\hat{\mathbf{z}}\right]\cos {\omega }_{c}t$$

$$E_m(x,z)=\sum _{m_2}\frac{\hslash |{\Omega }_{1,m}^{2,m_2}(x,z){|}^2}{4{\Delta }_{1,m}^{2,m_2}},$$

$$\left(\begin{array}{ccc}{E}_1(x,z)& 0& 0\\ 0& E_0(x,z)& 0\\ 0& 0& E_{\bar{1}}(x,z)\end{array}\right)=\alpha (x,z)\hat{1}+\beta (x,z){\hat{\mathbf{F}}}_z+q_M(x,z){\hat{\mathbf{F}}}_z^2,$$

$$q_M(x,z)=\frac{1}{2}\left[E_1(x,z)+E_{\bar{1}}(x,z)-2E_0(x,z)\right]$$

$$q_{DC}=\frac{(g{\mu }_B{)}^2}{\Delta W(1+2I{)}^2}\times B_0^2\approx {277Hz/G}^2\times B_0^2$$

$$q(x,z,t)=q+2{\Omega }_c(z)\cos (k_{T}x-\omega t)$$

$$\hat{H}=\frac{{\bf p}^{2}}{2m}{\hat{\bf 1}}+V_{trap}(z){\hat{\bf 1}}+\bigg{[}q+2\Omega_{c}(z)\cos(k_{T}x-\omega t)\bigg{]}\hat{F}_{z}^{2},$$

Full text

Ultracold atoms are a unique platform for exploring multi-faceted quantum magnetic behavior associated with spin. Some of the success stories in this arena include spinor BECs ueda-2013 , where magnetic interactions play an important role, as well as systems with artificial spin-orbit coupling spielman-2009 ; dalibard-2010 ; spielman-2011 ; sademelo-2011 ; pan-2014 ; ketterle-2017 ; zhai-2015 ; demarco-2015 ; zhang-2019 ; ye-2017 ; ye-2018 ; campbell-2016 , where independent-particle effects are primarily involved. Yet a comprehensive experimental framework linking these two disparate regimes of spin physics in ultracold gases has been lacking. In part, this is due to the fact that some of the richest behavior in spinor gases involves the dynamics of spin-nematic phases ketterle-1998 ; machida-1998 ; zhou-2004 ; demler-2003 ; affleck-2004 ; lett-2007 ; lett-2009 ; raman-2011 ; gerbier-2012 ; gerbier-2016 ; borgh-2014 ; symes-2017 ; kang-2019 . These phases are special because they have a vanishing total magnetization vector $\langle\hat{\bf F}\rangle=0$ and their order parameter is tensorial. For a spin-1 system, the expectation value of the spin-quadrupole tensor operator ${\hat{\bf Q}}_{ij}=\frac{1}{2}\left({\hat{\bf F}}_{i}{\hat{\bf F}}_{j}+{\hat{\bf F}}_{j}{\hat{\bf F}}_{i}\right)$ may act as an order parameter, where $i,j$ are the $\{x,y,z\}$ components of the spin-operator ${\hat{\bf F}}$ andreev-1984 . Through interactions between atoms, such tensor objects naturally generate spin entanglement and strong correlations. An important example of this is the reaction between two $|F=1,m=0\rangle$ alkali atoms through $s$-wave scattering, that is $|1,0\rangle+|1,0\rangle\leftrightarrow|1,1\rangle+|1,-1\rangle$, which conserves $m_{1}+m_{2}=0$ of atoms 1 and 2 sadler-2006 ; lucke-2011 ; gross-2011 ; bookjans-2011 ; vinit-2013 ; vinit-2018 . By contrast, the spin-orbit coupling achieved using Raman laser schemes does not readily lend itself to the study of pure spin-nematic objects, although a variety of other interacting many-body phases have been predicted galitski-2008 ; ho-2011 ; stringari-2012 ; baym-2012 ; stringari-2013 ; yamamoto-2017 .

In contrast to spin-orbit coupling, in this work we explore nematic-orbit coupling, where the linear momentum of spin-1 bosonic atoms is coupled to the spin-nematic degrees of freedom. Nematic spinor states have a zero expectation value for the spin vector $\langle{\hat{\bf F}}\rangle$ and nonzero quadrupole tensor $\langle{\hat{\bf Q}}_{ij}\rangle=\delta_{ij}-d_{i}d_{j}$, where ${\bf d}$ is the director. Easy axis or easy plane states correspond to ${\bf d}$ aligned with either the $z$ direction or lying in the $xy$-plane, respectively. Here, we propose an experimental setup to create nematic-orbit coupling between the center of mass of spin-1 bosons and the $zz$ component of the spin-quadrupolar operator ${\hat{\bf Q}}_{zz}={\hat{\bf F}}_{z}^{2}$, as shown in Fig. 1.

In the setup shown in Fig. 1, a spatio-temporally varying quadratic Zeeman shift $q({\bf r},t){\hat{\bf F}}_{z}^{2}$ is created using a combination of a static bias field and a microwave field that is produced by a monolithic microwave integrated circuit (MMIC) treutlein-2009 . After eliminating constant and linear terms in $\hat{\bf F}_{z}$ (see supplementary-material ), the effective independent particle Hamiltonian is

[TABLE]

where $\psi^{\dagger}_{a}({\bf r})$ is the creation operator of bosons at position ${\bf r}$ with spin components $a=\{\pm 1,0\}$, ${\bf p}^{2}/2m$ is the kinetic energy, $V({\bf r})=V_{trap}(z)$ is the trap potential, $q({\bf r},t)=q+2\Omega_{c}(z)\cos(k_{T}x-\omega t)$ is the resulting spatio-temporal modulation of the quadratic Zeeman shift with period $\lambda_{T}=2\pi/k_{T}$ and $\hat{\bf 1}$ is the identity matrix. The modulation amplitude $\Omega_{c}(z)=\Omega_{0}+\Omega_{1}z$ defines the strength of the nematic-orbit coupling. Since $\Omega_{c}(z)$ varies linearly with the $z$-coordinate, it couples two discrete energy levels $\epsilon_{1},\epsilon_{2}$ with different parity, which are defined by the spin-independent trapping potential $V({\bf r})$. A resonance condition for the magnetic traveling wave can be achieved when $\omega\approx\omega_{12}\equiv(\epsilon_{2}-\epsilon_{1})/\hbar$ supplementary-material . Given the discrete nature of the spectrum along $z$, we write the field operators as $\psi_{a}({\bf r})=\sum_{n}\varphi_{n}(z)\psi_{n,a}(x,y)$, where $\varphi_{n}(z)$ is the eigenfunction of trap state $n=\{1,2\}$. Within the rotating wave approximation (RWA) and zero detuning $\omega-\omega_{12}=0$, the Hamiltonian can then be rewritten in momentum space as (see supplementary-material ):

[TABLE]

Here, $\hat{\phi}^{\dagger}_{{\bf k}_{\perp}n}=\left[\phi_{n,1}({\bf k}_{\perp}),\phi_{n,0}({\bf k}_{\perp}),\phi_{n,\bar{1}}({\bf k}_{\perp})\right]$ is the spinor creation operator with subscript $\bar{1}$ as a shorthand for $-1$, ${\bf k}_{\perp}=(k_{x},k_{y})$, ${\bf H}_{D}=\varepsilon_{\bf k}{\hat{\bf 1}}+q{\hat{\bf F}}^{2}_{z},$ where $\varepsilon_{\bf k}=\hbar^{2}k_{\perp}^{2}/(2m)$ is the kinetic energy with $k_{\perp}=|{\bf k}_{\perp}|$, and ${\bf k}_{\pm}={\bf k}_{\perp}\pm(k_{T}/2){\hat{\bf x}}$ are shifted momenta. The Hermitian conjugate (H.c.) term is $\Omega\hat{\phi}^{\dagger}_{{\bf k}_{+},2}{\hat{\bf F}}^{2}_{z}\hat{\phi}_{{\bf k}_{-},1}$, where $\Omega=\int\mathrm{d}z\varphi^{*}_{1}(z)[\Omega_{1}z]\varphi_{2}(z)$ plays the role of a Rabi frequency (see supplementary-material ). The diagonalization of Eq. (2) leads to a trivial eigenvalue $E_{0}=\hbar^{2}k_{\perp}^{2}/(2m)$ corresponding to spin component $a=0$, and to non-trivial eigenvalues

[TABLE]

The lower (higher) energy branch is labeled by $\alpha$ $(\beta)$, with corresponding eigenvectors

[TABLE]

written as linear combinations of $\phi_{1,a}({\bf k}_{-})$ and $\phi_{2,a}({\bf k}_{+})$. Expressions for the coefficients $u_{\pm\alpha}({\bf k}_{\perp})$ and $u_{\pm\beta}({\bf k}_{\perp})$ are found in supplementary-material . The absolute minimum of all eigenvalues, where Bose-Einstein condensation occurs, depends on parameters $q$ and $\Omega$, and is found in the lower band $\alpha$. We locate the minima of these energy bands by extremizing with respect to $k_{x}$. We work with dimensionless variables and set $k_{T}$ as the unit of momentum and $E_{T}=\hbar^{2}k^{2}_{T}/(2m)$ as the unit of energy. The scaled parameters are $\widetilde{q}=q/E_{T}$, $\widetilde{\Omega}=\Omega/E_{T}$ and ${\widetilde{\bf k}}_{\perp}={\bf k}_{\perp}/k_{T}$.

In Fig. 2, we show the phase diagram of $\widetilde{q}$ versus $\widetilde{\Omega}$ arising from Eq. (3). The dashed-green line corresponds to the phase boundary ${\widetilde{q}}_{c}({\widetilde{\Omega}})={\widetilde{\Omega}}^{2}$ for ${\widetilde{\Omega}}<1/2$, that separates an easy-axis nematic BEC at ${\widetilde{\bf k}_{\perp}}={\bf 0}$ for spin component $a=0$, when ${\widetilde{q}}>{\widetilde{q}}_{c}({\widetilde{\Omega}})$, from a double-well easy-plane nematic BEC for spin components $a=\pm 1$, when ${\widetilde{q}}<{\widetilde{q}}_{c}({\widetilde{\Omega}})$. The dotted-red line describes the phase boundary ${\widetilde{q}}_{c}({\widetilde{\Omega}})={\widetilde{\Omega}}-1/4$ for ${\widetilde{\Omega}}>1/2$, that separates an easy-axis BEC at ${\widetilde{\bf k}_{\perp}}={\bf 0}$ for spin component $a=0$, when ${\widetilde{q}}>{\widetilde{q}}_{c}({\widetilde{\Omega}})$, from a single-well easy-plane nematic BECs for spin components $a=\pm 1$, when ${\widetilde{q}}<{\widetilde{q}}_{c}({\widetilde{\Omega}})$. The solid-blue line ${\widetilde{\Omega}}=1/2$ separates the easy-plane nematic BECs in the $\alpha$ band into two sectors: a) a double-well phase where condensation occurs at finite momenta $({\widetilde{k}}_{x},{\widetilde{k}}_{y})=(\pm{\widetilde{k}}_{0},0)$, with ${\widetilde{k}}_{0}=\sqrt{1/4-{\widetilde{\Omega}}^{2}}$, and b) a single-well phase where condensation occurs at zero momentum $\widetilde{\bf k}_{\perp}={\bf 0}$. The solid-black dot at coordinates $({\widetilde{q}},{\widetilde{\Omega}})=(1/4,1/2)$ represents a triple point.

Next, we discuss the interaction Hamiltonian ${\hat{H}}_{\rm int}={\hat{H}}_{0}+{\hat{H}}_{2}$. The first term is the spin-independent interaction ${\hat{H}}_{0}=(c_{0}/2L^{2}_{\perp}){\widetilde{H}}_{0}$, with

[TABLE]

where the subscripts $\{n_{i}\}$ denote the set of trapped states with quantum numbers $(n_{1},n_{2},n_{3},n_{4})$ that label the coefficients $C_{\{n_{i}\}}=\int\mathrm{d}z\varphi^{*}_{n_{1}}(z)\varphi^{*}_{n_{2}}(z)\varphi_{n_{3}}(z)\varphi_{n_{4}}(z)$. In Eq. (5), the momenta are ${\bf k}_{p\pm}={\bf k}_{\perp}\pm{\bf p}_{\perp}/2$ and ${\bf k}^{\prime}_{p\pm}={\bf k}^{\prime}_{\perp}\pm{\bf p}_{\perp}/2$, and the operators are

[TABLE]

In the interaction Hamiltonain, the second term is the spin-dependent interaction ${\hat{H}}_{2}=(c_{2}/2L^{2}_{\perp}){\widetilde{H}}_{2}$, with

[TABLE]

where the vector operators

[TABLE]

contain the spin-1 matrices ${\hat{\bf F}}$.

The Hamiltonians ${\hat{H}}_{\rm IP}+{\hat{H}}_{\rm int}$ preserve the magnetization $m_{z}=n_{+1}-n_{-1}$, where $n_{\pm 1}$ is the density of bosons with spin component $a=\pm 1$, that is, $m_{z}$ is a conserved quantity of the total Hamiltoninan. From now on, we consider only $m_{z}=0$, in which case a phase transition occurs at $\widetilde{q}_{c}=0$ between the easy-plane nematic state $|\zeta_{P}\rangle$ $({\widetilde{q}}<{\widetilde{q}}_{c})$ with spin-densities $n_{0}=0$, $n_{+1}=n_{-1}\neq 0$, and the easy-axis nematic state $|\zeta_{A}\rangle$ $({\widetilde{q}}>{\widetilde{q}}_{c})$ with spin-densities $n_{0}\neq 0$, $n_{+1}=n_{-1}=0$, as shown in Fig. 2, when ${\widetilde{\Omega}}=0$ ketterle-1998 ; lett-2007 ; raman-2011 ; ueda-2013 ; gerbier-2016 .

The effects of nematic-orbit coupling are also present in the collective excitations. First, we investigate the easy-axis nematic phase, where condensation occurs at $\widetilde{\bf k}_{\perp}={\bf 0}$ for spin projection $a=0$. The Bogoliubov excitation spectrum is then identical to a scalar condensate, $\varepsilon_{b}({\bf k}_{\perp})=\left[\varepsilon_{\bf k}\left(\varepsilon_{\bf k}+2c_{0}n_{c}\right)\right]^{1/2},$ where $n_{c}$ is the total particle density and $\varepsilon_{\bf k}=\hbar^{2}k_{\perp}^{2}/(2m)$ is the kinetic energy.

Next, we consider the easy-plane nematic phase in the single-well regime when $\widetilde{q}\ll\widetilde{\Omega}-1/4$ and $\widetilde{\Omega}>0.5$. We write the field operators $\phi_{n,a}$ in terms of $\chi_{a\alpha},\chi_{a\beta}$ as shown in supplementary-material . Condensation occurs at ${\widetilde{\bf k}}_{\perp}={\bf 0}$ for the $\alpha$-band only, thus we drop the $\alpha$ index from our notation. The resulting Bogoliubov Hamiltonian is

[TABLE]

The matrices for spin preserving processes are

[TABLE]

where $a=\{+1,-1\}$ is represented by $\{1,{\bar{1}}\}$, $E_{g}({\bf k}_{\perp})=E_{\alpha}({\bf k}_{\perp})-E_{\alpha}(0)$ is a measure of the excitation energy of independent particles with respect to the minimum of the $\alpha$-band, $\Phi_{a}$ is the spin-dependent phase of the condensate in the $\alpha$-band at ${\bf k}_{\perp}={\bf 0}$ and $c,f$ are proportional to the spin-preserving interaction energy $(c_{0}+c_{2})n_{c}$. The matrices for spin-flip processes are

[TABLE]

and ${\bf D}^{\dagger}$, where $d$ and $g$ are proportional to the spin-flip interaction energy $(c_{0}-c_{2})n_{c}$. Lastly, in Eq. (9), $G_{\rm sw}$ is the ground state energy and ${\bf X}^{\dagger}_{\bf k}=\begin{pmatrix}\chi^{\dagger}_{1}({\bf k}_{\perp})&\chi_{1}(-{\bf k}_{\perp})&\chi^{\dagger}_{\bar{1}}({\bf k}_{\perp})&\chi_{\bar{1}}(-{\bf k}_{\perp})\end{pmatrix}$ is a vector operator, where $\chi^{\dagger}_{a}$ represents the creation operator in the $\alpha$-band.

The positive eigenvalues in units of $E_{T}$ are

[TABLE]

where $\widetilde{E}_{g}({\bf k}_{\perp})=E_{g}({\bf k}_{\perp})/E_{T}$ is a dimensionless independent particle energy, $\widetilde{c}=(c_{0}+c_{2})n_{c}A_{\alpha}({\bf k}_{\perp})/(4E_{T}),$ $\widetilde{f}=(c_{0}+c_{2})n_{c}B_{\alpha}({\bf k}_{\perp})/(4E_{T})$ are dimensionless spin-preserving interaction energies and $\widetilde{d}=(c_{0}-c_{2})n_{c}A_{\alpha}({\bf k}_{\perp})/(4E_{T}),$ $\widetilde{g}=(c_{0}-c_{2})n_{c}B_{\alpha}({\bf k}_{\perp})/(4E_{T}),$ are dimensionless spin-flip interaction energies. Here, $A_{\alpha}({\bf k}_{\perp})=5/2+|\widetilde{\Omega}|/\left[\sqrt{{\widetilde{k}}_{x}^{2}+{\widetilde{\Omega}}^{2}}\right]$ and $B_{\alpha}({\bf k}_{\perp})=2+3|\widetilde{\Omega}|/\left[2\sqrt{{\widetilde{k}}_{x}^{2}+{\widetilde{\Omega}}^{2}}\right]$ describe the anisotropic nature of the interactions induced by the nematic-orbit coupling. When ${\widetilde{d}}={\widetilde{g}}=0$, that is, $c_{0}=c_{2}$, the matrix ${\bf D}$ of spin-flip processes vanishes and the spin-sectors $\{1,{\bar{1}}\}$ are uncoupled leading to two degenerate linear modes at low momenta. Assuming that $c_{0}>c_{2}>0$ like in ${}^{23}{\rm Na}$, we can understand a few limits from Eq. (12). In the first mode, the sum ${\widetilde{c}}+{\widetilde{d}}$ and ${\widetilde{f}}+{\widetilde{g}}$ are proportional to the spin-independent interaction parameter $c_{0}$, while in the second mode, the difference ${\widetilde{c}}-{\widetilde{d}}$ and ${\widetilde{f}}-{\widetilde{g}}$ are proportional to the spin-dependent interaction parameter $c_{2}$. Thus, the first mode is associated with density-density interactions $c_{0}$, while the second is associated with spin-spin interactions $c_{2}$. We plot the excitation spectra ${\widetilde{\epsilon}}_{b,1}({\bf k}_{\perp})$ and ${\widetilde{\epsilon}}_{b,2}({\bf k}_{\perp})$ versus $k_{x}$ in Fig. 3(a) and versus $k_{y}$ in Fig. 3(b), with $c_{0}$ and $c_{2}$ values for ${}^{23}{\rm Na}$ ho-1998 .

Lastly, we consider the easy-plane nematic phase in the double-well region, when $\widetilde{q}\ll{\widetilde{\Omega}}^{2}$ and $\widetilde{\Omega}<0.5$. Condensation occurs in two degenerate minima at $\pm k_{0}{{\bf\hat{x}}}$ of the $\alpha$-band. There are four excitation modes involving left $(L)$ and right $(R)$ wells and spin sectors $a=\{1,{\bar{1}}\}$. The Bogoliubov Hamiltonian becomes

[TABLE]

where ${\bf Y}^{\dagger}_{\bf k}=\begin{pmatrix}{\bf X}^{\dagger}_{L}({\bf k}_{\perp})&{\bf X}^{\dagger}_{R}({\bf k}_{\perp})\end{pmatrix}$ is an eight-dimensional vector with four dimensional components ${\bf X}^{\dagger}_{j}({\bf k}_{\perp})=\begin{pmatrix}\chi^{\dagger}_{j1}({\bf k}_{\perp})&\chi_{j1}(-{\bf k}_{\perp})&\chi^{\dagger}_{j\bar{1}}({\bf k}_{\perp})&\chi_{j\bar{1}}(-{\bf k}_{\perp})\end{pmatrix}$ in the $j=\{L,R\}$ sectors, and $G_{\rm dw}$ is the ground state energy. The ${\bf M}_{ij}$ matrices are given in supplementary-material and the excitation spetrum is obtained numerically, but a qualitative understanding is possible. In each well there are equal numbers of atoms with spin components $a=\{1,{\bar{1}}\}$, that is, $n_{1L}=n_{1R}$ and $n_{{\bar{1}}L}=n_{{\bar{1}}R}$. When all interactions are present and all atoms oscillate in phase, this excitation corresponds to a center-of-mass motion with linear dispersion and lowest energy at low momenta, which is also anisotropic since the effective mass is heavier along $k_{x}$. When atoms with the same spin-projection $a$ oscillate in phase in both L and R wells, but out of phase with respect to their spin-projections, then a second linear mode arises with larger (larger) velocity along $k_{x}$ $(k_{y})$ in comparison to the center-of-mass mode. When the spin-spin interactions are neglected and atoms with spin-projection $a$ oscillate out of phase in $L$ and $R$ wells they produce two degenerate linearly dispersing modes. However, when spin-spin interactions are included the degeneracy of these modes is lifted producing a linearly dispersing mode with lower (higher) energy when the relative motion of $1$ and ${\bar{1}}$ is in (out of) phase. All four modes ${\widetilde{\epsilon}}_{b,1}({\bf k}_{\perp})$, ${\widetilde{\epsilon}}_{b,2}({\bf k}_{\perp})$, ${\widetilde{\epsilon}}_{b,3}({\bf k}_{\perp})$ and ${\widetilde{\epsilon}}_{b,4}({\bf k}_{\perp})$ of the excitation spectrum are shown in Fig. 3(c) and 3(d) for ${}^{23}{\rm Na}$ parameters.

Next, we analyze manifestations of the nematic-orbit coupling in real space and focus on the easy-plane nematic phases with $n_{0}=0$ and $n_{+1}=n_{-1}\neq 0$. Far below the phase boundary ${\widetilde{q}}_{c}({\widetilde{\Omega}})$, the effective Hamiltonian is $\hat{H}_{\rm EP}={\hat{H}}^{\prime}_{IP}+{\hat{H}}_{\rm I},$ with

[TABLE]

where ${\hat{\psi}}^{*}_{n}=\left[\psi^{*}_{n,1}({\bf r}_{\perp}),\psi^{*}_{n,0}({\bf r}_{\perp}),\psi^{*}_{n,{\bar{1}}}({\bf r}_{\perp})\right]$ represents the 2D condensate wavefunction in trap states with quantum number $n$. The interaction Hamiltonian is ${\hat{H}}_{\rm I}=\int\mathrm{d}^{3}r{\hat{\cal H}}_{I}$, where

[TABLE]

with $c_{0}>c_{2}>0$ as in ${}^{23}{\rm Na}$, leading to the same local condensate densities, that is, $|{\bf\Psi}_{1}({\bf r})|^{2}=|{\bf\Psi}_{\bar{1}}({\bf r})|^{2}$.

In the single-well phase, condensation occurs in the $\alpha$-band at ${\widetilde{\bf k}}_{\perp}={\bf 0}$. However, the wavefunction $\psi_{a}({\bf r})$ in real space is a linear combination of momentum shifted $(\pm(k_{T}/2){\hat{\bf x}})$ condensates with relative phase $\vartheta$ supplementary-material , resulting in a spatial variation of the form

[TABLE]

where $\varphi_{1,2}(z)$ are the trap states along $z$ direction and its period $\lambda_{h}=2\pi/(k_{T}/2)=2\lambda_{T}$ commensurate to the period $\lambda_{T}$ of the periodic potential $q({\bf r},t)$. The phase $\vartheta=0$ supplementary-material is detemined by minimization of the free energy and ${\cal A}_{\rm sw}$ is obtained by normalizing the condensate density $n_{C}({\bf r})=\sum_{a=\{1,{\bar{1}}\}}|\psi_{a}({\bf r})|^{2}$ to the total number of condensed particles $N_{C}$ supplementary-material . Therefore, the dimensionless local condensate density ${\widetilde{n}}_{C}(\widetilde{x})$ at some fixed $\widetilde{z}_{0}$, describing a easy-plane single-period nematic density wave (SPNDW), can be obtained by squaring the norm of Eq. (16) supplementary-material . ${\widetilde{n}}_{C}({\widetilde{x}})$ for $\sigma=0.7$, ${\widetilde{\Omega}}=1$ and $\widetilde{z}=\pi/8$ is plotted in Fig. 4(a), where $\widetilde{x}=k_{T}x$, $\widetilde{z}=(2\pi/L_{z})z$ and $\sigma=N_{C}/N$ is the condensate fraction. It is uniform apart from the periodic variation at the lattice period $\lambda_{T}$.

In the double-well phase, condensation occurs in the $\alpha$-band at ${\widetilde{\bf k}}_{\perp}=\pm{\widetilde{k}}_{0}{\hat{\bf x}}$. Thus, the wavefunction $\psi_{a}({\bf r})$ in real space is a linear combination of two single-well condensates with momenta $(k_{0}\pm k_{T}/2){\hat{\bf x}}$ and phases $\vartheta$, $\vartheta_{LR}$ supplementary-material , resulting in a spatial variation of the form

[TABLE]

with two periods $\lambda_{\pm}=2\pi/|k_{0}\pm k_{T}/2|$, which are generically incommensurate with $\lambda_{T}$. Here, we denote ${\cal A}^{\prime}_{\rm dw}={\cal A}_{\rm dw}e^{-i\frac{\vartheta+\vartheta_{LR}}{2}}$, $\varphi_{-}(z)=\varphi_{1}(z)$ and $\varphi_{+}(z)=\varphi_{2}(z)$ for simplicity. The relative phase $\vartheta$, $\vartheta_{LR}$ were determined by minimizing the free energy numerically supplementary-material , resulting in $\vartheta=0$. The energy functional contains a rapid oscillation at the underlying period $\lambda_{T}$ as the system size $L_{\perp}$ is varied supplementary-material . We chose $k_{T}L_{\perp}=250$ and $\vartheta_{LR}=0$ to minimize the energy over this oscillation, with the results shown in Fig. 4(b). $\vartheta_{LR}=\pi$ achieved similar results for other $k_{T}L_{\perp}$. By squaring the wave function of Eq. (17), this leads to the dimensionless condensate density describing a double-period nematic density wave (DPNDW) along $x$ direction shown in Fig. 4(b) for $\widetilde{z}=\pi/8$ (see supplementary-material ).

In conclusion, we have proposed a mechanism for the creation of nematic-orbit coupling in spin-1 condensates and uncovered their phase diagram and excitation spectra. Our work connects orbital motion of atoms to the rich physics of spin-nematics, and opens up a new direction to explore strongly correlated spin-nematic states. Future work may include higher spin systems and coupling to other tensor components ${\hat{\bf Q}}_{ij}$. Extension to higher dimensions could allow nontrivial topology to be explored, analogous to half-quantum vortices in ordinary nematics seo-2015 , which have parallels in solid state systems  lagoudakis-2009 ; maeno-2011 .