Bosons with incommensurate potential and spin-orbit coupling

TL;DR

This paper explores the phase diagram of spin-half bosons in a 1D optical lattice with incommensurate potential and spin-orbit coupling, revealing novel localized phases, spin-split momentum distributions, and spectral statistics transitions.

Contribution

It provides new insights into the effects of spin-orbit coupling and Aubry-André potential on bosonic phases, including the identification of a spin-split momentum distribution as a transition signature.

Findings

Spin-split momentum distribution in localized phase

Transition of level statistics from GUE to GOE

Identification of density wave phase with NNI

Abstract

We chart out the phase diagram of ultracold `spin-half' bosons in a one-dimensional optical lattice in the presence of Aubry-Andr\'e (AA) potential and with spin-orbit (SO) and Raman couplings investigating the transition from superfluid (SF) to localized phases and the existence of density wave phase for nearest-neighbor interaction (NNI). We show that the presence of SO coupling and AA potential leads to a novel spin-split momentum distribution of the bosons in the localized phase near the boundary with the SF phase, which can act as a signature of such a transition. We also obtain the level statistics of the bosons in the superfluid phase with finite NNI and demonstrate its change from Gaussian Unitary Ensemble (GUE) to Gaussian Orthogonal Ensemble (GOE) as a function of the Raman coupling. We discuss experiments which can test our theory.

Equations

$$\hat{H}$$

$$E_k^{\pm }=-2\cos k\cos q\pm 2{\left[{\sin }^{2}k{\sin }^{2}q+{\Omega }^2/4\right]}^{1/2}$$

$${\Omega }_c$$

$$f_s=N_s^2(E[\theta ]-E[0])/(N_p{\theta }^2),$$

$$I=\sum _l(|{\psi }_{l,\uparrow }{|}^2+|{\psi }_{l,\downarrow }{|}^2{)}^2$$

$$n_{k\sigma }=\sum _{l,l^{\prime }}\exp \{ik(l-l^{\prime })\}⟨c_{l\sigma }^{†}{c}_{l^{\prime }\sigma }⟩/N_s$$

$$H_{SO}=\left(\begin{array}[]{cc}-2\cos(k+q)&\Omega\\ \Omega&-2\cos(k-q)\end{array}\right)$$

$$m_{\Omega >{\Omega }_c}^{\ast }\equiv m_{>}^{\ast }={\left(\frac{{\partial }^2{E}_k^{-}}{\partial k^2}\right)}_{k=0}=(1-{\Omega }_c/\Omega )\cos q.$$

$$H_{\theta }=-t\sum _{l,\sigma }\left({\hat{b}}_{l,\sigma }^{†}{e}^{iq{\hat{\sigma }}_z}{\hat{b}}_{l+1,\sigma }{e}^{-i\theta /N_s}+h.c.\right)+\Omega \sum _{l,\sigma }{\hat{b}}_{l,\sigma }^{†}{\hat{\sigma }}_x{\hat{b}}_{l,\sigma }+\lambda \sum _{l,\sigma }\cos (2\pi \beta l){\hat{b}}_{l,\sigma }^{†}{\hat{b}}_{l,\sigma }.$$

$$H_{\theta }=H_0+\frac{\theta }{N_s}\hat{J}-\frac{{\theta }^2}{2N_s^2}\hat{T}$$

$$f_s=-\frac{1}{2t}⟨{\psi }_0|\hat{T}|{\psi }_0⟩-\frac{1}{t}\sum _{\nu \ne 0}\frac{|⟨{\psi }_{\nu }|\hat{J}|{\psi }_0⟩{|}^2}{E^{\nu }-E^0}$$

$$H_1=\frac{\lambda }{2}\sum _k{\hat{b}}_{k\sigma }^{†}({\hat{b}}_{k+2\pi \beta \sigma }+{\hat{b}}_{k-2\pi \beta \sigma }).$$

$$\lambda =\delta \Omega (1-\tan q/\sqrt{{\sin }^2\beta {\sin }^{2}q})+2\cos q(1-\cos \beta )+|2\sin q\tan q|-|\sin (q)|\sqrt{{\sin }^2\beta +{\tan }^{2}q}.$$

$$\lambda \simeq 2\cos q(1-\cos \beta )$$

$$-({\psi }_{l+1,\uparrow }{e}^{iq}+{\psi }_{l-1,\uparrow }{e}^{-iq})+\lambda \cos (2\pi \beta l){\psi }_{l,\uparrow }+\Omega {\psi }_{l,\downarrow }+U(|{\psi }_{l,\uparrow }{|}^2+|{\psi }_{l,\downarrow }{|}^2){\psi }_{l,\uparrow }=\mu {\psi }_{l,\uparrow }$$

$$-({\psi }_{l+1,\downarrow }{e}^{-iq}+{\psi }_{l-1,\downarrow }{e}^{iq})+\lambda \cos (2\pi \beta l){\psi }_{l,\downarrow }+\Omega {\psi }_{l,\uparrow }+U(|{\psi }_{l,\uparrow }{|}^2+|{\psi }_{l,\downarrow }{|}^2){\psi }_{l,\downarrow }=\mu {\psi }_{l,\downarrow }$$

$$\psi_{l}=\mathcal{N}e^{-|l|/\xi}\left(\begin{array}[]{c}e^{ikl}\\ -e^{-ikl}\\ \end{array}\right),\quad\mathcal{N}=\sqrt{\frac{\tanh(1/\xi)}{2}}$$

$$E(k)=-\left[\frac{\cos (k-q)}{\cosh (1/\xi )}+\Omega \tanh (1/\xi )\frac{\sinh (2/\xi )}{\cosh (2/\xi )-\cos 2k}\right].$$

$$\psi^{l}=\sqrt{n_{0}^{l}}\left(\begin{array}[]{c}\cos\theta^{l}e^{i\phi_{\uparrow}^{l}}\\ \sin\theta^{l}e^{i\phi_{\downarrow}^{l}}\\ \end{array}\right)$$

$$|\psi (t)⟩$$

$$⟨\psi (t)|O|\psi (t)⟩=\sum _{m,n}{c}_m^{\ast }{c}_n{e}^{i({ϵ}_m-{ϵ}_n)t/\hslash }⟨m|O|n⟩$$

$$\mathcal{I}=\frac{N_o-N_e}{N_{tot}}$$

Full text

The study of localization phenomena in correlated systems has regained a new interest recently in the context of many-body localization (MBL) mbl1 ; mbl2 . Ultracold atoms in optical lattices, which act as emulators of strongly correlated model Hamiltonians boserev1 , can serve as test beds for such phenomena boserev2 . In this context, systems with quasiperiodic potentials, which have posed several interesting theoretical challenges over may decades quasiref1 ; quasiref2 ; quasiref3 , turn out to be particularly relevant. A model Hamiltonian describing such a quasiperiodic system is the well-known Aubry-André (AA) model aapaper , which, unlike the Andreson model, exhibits localization transition in 1D aapaper ; aapaper1 . This property of the AA model has generated an impetus to study MBL quasirev1 ; mblaa2 . Moreover, experimental realization of the AA model in bichromatic optical lattice has led to observation of localization of both lightloclight1 and ultracold matter wave inguscio ; kushref1 .

In recent past, extensive research on the Bose-Hubbard (BH) model using ultracold bosonic atoms in optical lattices paved the way for studying the effect of interactions on localization phenomenon leading to possible glassy phases inguscio1 ; demarco ; bhref1 ; Roth . In addition, intense theoretical studies has also been carried out on the BH model in the presence of Abelian and non-Abelian gauge fields; such gauge fields have been experimentally realized in atom-laser systems ab1 ; nonab1 . Such systems allow for observation of several exciting phenomena abph1 ; abph2 ; nonabph1 ; most interestingly, they enable us to study strongly interacting bosons in the presence of tunable spin-orbit (SO) coupling sorefs1 ; sorefs2 ; sorefs3 ; sorefs4 . The realization of the AA model in bichromatic lattice and the creation of SO interactions for ultracold bosons therefore provides an unique opportunity to study localization phenomenon induced by the AA potential in presence of tunable SO interactions.

In this work, we study a two-species Bose-Hubbard model coupled by Raman frequency $\Omega$, in the presence of an AA potential and a SO coupling and show that such a system leads to several novel features which appear only in the presence of both the AA potential and the SO coupling. The central results of our study are as follows. First, we chart out the phase diagram of 1D ultracold bosons in an optical lattice and demonstrate the existence of density wave (DW), superfluid (SF), and localized phases and study the transition between these phases. Second, we show that for sufficiently high $\Omega$, the bosons in the presence of both the AA potential and the SO coupling exhibits a spin-split momentum distribution in the localized phase, near the boundary with the SF phase, irrespective of the strength of their interaction. Such a splitting can therefore serve as a signature of this transition. We note that this spin splitting does not occur in the absence of either the AA potential or the SO coupling. Third, we study the level statistics of the bosons in the strongly interacting regime, where the presence of AA potential and Raman coupling $\Omega$ between the spins play a crucial role in changing the spectral statistics between different universality classes of random matrix theory (RMT). Apart from poissonian level spacing distribution in the localized regime, we find that the level statistics change continually from GUE ($\Omega=0$) to GOE as a function of $\Omega$. We identify the additional symmetry at the $\Omega=0$ point which is behind this change. Finally, we discuss experiments which can test our theory.

The Hamiltonian of the bosons in a bi-chromatic 1D lattice with AA potential and SO coupling is given by

[TABLE]

where, $\hat{b}_{l,\sigma}^{\dagger}$ and $\hat{n}_{l\sigma}=\hat{b}_{l,\sigma}^{\dagger}\hat{b}_{l,\sigma}$ are the creation and the density operator of the bosons of (pseudo)spin $\sigma$ at the lattice site $l$, $\hat{n}_{l}=\sum_{\sigma}\hat{n}_{l\sigma}$, $\bar{\sigma}=\downarrow(\uparrow)$ for $\sigma=\uparrow(\downarrow)$, $t$ is the hopping strength, $q$ is the SO coupling strength, $\Omega$ is the Raman frequency and $\lambda$ denotes the strength of the quasiperiodic potential. In the rest of the paper we consider nearest neighbor and on-site interactions with coupling strengths: $V=\mathcal{V}_{l,l+1}$ and $U=\mathcal{V}_{l,l}$ respectively; in what follows, we shall scale all energies in unit of $t$.

Non-interacting limit: We first look into the the non-interacting bosons by setting $U=V=0$ in Eq. 1. For $\lambda=0$, this reduces to a pure SO coupled bosonic system with single particle spectrum

[TABLE]

and the eigenstates are given by, $\psi_{k}^{\pm}=e^{ikx}(\cos(\theta_{k}-\pi/4\mp\pi/4),\sin(\theta_{k}-\pi/4\mp\pi/4))^{T}$, where $\cos\theta_{k}=[1/2+[4+\Omega^{2}/(\sin^{2}k\sin^{2}q)]^{-1/2}]^{1/2}$. We find that there exists a critical Raman coupling given by,

[TABLE]

below which the ground state is doubly degenerate associated with the finite momenta $k_{0}=\pm\cos^{-1}[\cos q(1+\Omega^{2}/(4\sin^{2}q))^{1/2}]$. The doubly degenerate ground states are related by $\psi_{-k}=\hat{\sigma}_{z}\hat{\mathcal{T}}\psi_{k}$, where $\hat{\mathcal{T}}=-i\hat{\sigma}_{y}\hat{\mathcal{C}}$ is the time reversal symmetry (TRS) operator and $\hat{\mathcal{C}}$ is the complex conjugation operator.

Next we turn on $\lambda$ keeping $U=V=0$. For $q=\Omega=0$, above Hamiltonian is reduced to a two component AA model which undergoes a localization transition above a critical coupling strength $\lambda_{c}=2$. For two extreme regimes $q\neq 0,\Omega=0$ (pure SO coupling) and $q=0,\Omega\neq 0$ (strong Raman coupling) the single particle Hamiltonian preserves the self duality at $\lambda_{c}=2$ and all states are localized above $\lambda_{c}$. To study localization transition in the intermediate regime with $q\neq 0,\Omega\neq 0$, we numerically diagonalize the single particle Hamiltonian to obtain the ground state and the excitation spectrum. Since the duality does not hold in this regime a mobility edge appears and energy dependent localization occurs for eigenstates aaso . We focus on the localization transition of the ground state in presence of SO interaction and the variation of the critical disorder strength $\lambda_{c}$ on Raman coupling. We locate the change from the SF to the localized phase in two ways. First, we measure the superfluid fraction(SFF) by applying a phase twist Fisher $\theta\ll\pi$ at the boundary. In the presence of such a twist $t_{l,l+1}\rightarrow te^{i\theta/N_{s}}$ (Eq. 1). The SFF can then be computed as Roth ,

[TABLE]

where $E[\theta]$ is the ground state energy in presence of twist, $N_{s}$($N_{p}$) is the number of sites(particles).

The second measure of localization is the inverse participation ratio (IPR) of the ground state wavefunction defined as

[TABLE]

where $|\psi_{l,\sigma}|^{2}$ is the boson density of spin $\sigma$ at site $l$. As expected, we find that SFF decreases and the IPR increases with increasing $\lambda$ around the transition.

The phase diagram obtained from these computations is shown in Fig. 1(a) in $\lambda-\Omega$ plane for $q=0.3\pi$. We note that $\lambda_{c}$ decreases from its self-dual value $\sim 2$ for $q\neq 0$ and shows a dip at $\Omega_{c}$, which demarcates the delocalized phase in two regime. Below $\Omega_{c}$ the degeneracy of the ground state is lifted by the quasi-periodic potential; however the ground state has a net momentum and $s_{z}$ polarization. For $\Omega>\Omega_{c}$, the ground state wavefunction is spin-polarized along $\hat{x}$ and has vanishing net momentum. The behavior of $\lambda_{c}$ with $\Omega$ can be understood from the enhancement of effective mass of bosons in the lower branch $m^{\ast}=\partial^{2}E_{k}^{-}/\partial k^{2}|_{k=k_{0}}$ due to the combined effect of SO and Raman couplings. This in turn reduces the effective hopping strength $t_{\rm eff}=t/m^{\ast}$ of underlying AA model for which the critical strength for localization transition can be estimated as $\lambda_{c}\sim 2t_{\rm eff}\sim 2/m^{\ast}$ supp1 . We note that the idea of $m^{\ast}$ also quantitatively explains the the variation of SFF with $\Omega$ for $\lambda=0$ and that SFF decreases and the IPR increases with increasing $\lambda$ around the transition as expected supp1 .

To elucidate the role of the AA potential and the SO coupling in the transition, we compute the spin-resolved momentum distribution defined as

[TABLE]

where $k=2\pi m/N_{s}$, with $m\in[-N_{s}/2,(N_{s}-1)/2]$. For $\Omega>\Omega_{c}$, both $n_{k\uparrow}$ and $n_{k\downarrow}$ is peaked at $k=0$ in the delocalized phase as seen for standard superfluids. In contrast, in the localized regime near the transition, $n_{k\sigma}$ becomes spin dependent and is peaked at $k=k_{\sigma}^{\rm max}\neq 0$ preserving the symmetry $n_{\uparrow}(k)=n_{\downarrow}(-k)$ (see Fig. 1(b)). The splitting of these peaks are given by $\Delta_{k^{\rm max}_{\sigma}}=k_{\uparrow}^{\rm max}-k_{\downarrow}^{\rm max}\sim q$ leading to the conclusion that the split in $n_{k\sigma}$ arises from a finite SO coupling. As shown in Fig. 1(c), (d), $\Delta_{k^{\rm max}_{\sigma}}$ vanishes for either $q=0$ or $\lambda=0$; this shows the necessity of both the AA potential and the SO coupling for the peak splitting. We find that this splitting can be qualitatively understood from a variational wavefunction calculation and is associated with spin dephasing showing a fluctuation of relative phase between two spin components of the ground state wavefunction supp1 .

Hardcore limit: To explore the effect of interaction on this phenomenon we now set $U\to\infty$, keeping $V=0$. This limit facilities computation by imposing the constraint $n_{l}\leq 1$ at each site and allows us to perform exact diagonalization within a restricted Hilbert space of three states per site. We restrict our calculation to half-filled HC bosons, $\sum_{l}n_{l}=N_{s}/2$ so that we are always in the SF phase for $\lambda=V=0$. In addition to SFF we also compute the boson condensate fraction (BCF) since BCF and SFF are quite different for strongly interacting bosons and are important for characterizing the localized phases. We construct the one-body density matrix from the ground state $|\psi_{0}\rangle$: $\rho(l,\sigma;l^{\prime},\sigma^{\prime})=\langle\psi_{0}|\hat{b}_{l^{\prime},\sigma^{\prime}}^{\dagger}\hat{b}_{l,\sigma}|\psi_{0}\rangle$legg1 ; the largest eigenvalue $N_{c}$ of which gives the BCF $f_{c}=N_{c}/{\rm Tr}(\hat{\rho})$.

A plot of $f_{s}$ and $f_{c}$ in the $\Omega-\lambda$ plane for a fixed $q$ is shown in Fig. 2(a,b). These plots clearly indicate a regime for $\lambda>\lambda_{c}$ where $f_{s}$ vanishes but $f_{c}$ remains finite indicating localized phase of the bosons. Although in finite system there is no transition, the behavior of $\lambda_{c}$ obtained from SFF is similar to that for non-interacting bosons; however $\Omega_{c}$ shifts to a lower value. Near this boundary, particularly for $\Omega\leq\Omega_{c}$, there is clear indication of Bose-glass (BG) phase with $f_{s}=0$ and $f_{c}\neq 0$.

Finally, we compute the $n_{k\sigma}$ of the hardcore bosons. As shown in Figs. 2(c),(d), the splitting of the spin momentum peak occurs in the localized phase and survives in the hardcore limit. We have checked that $n_{k\sigma}$ are peaked at $k=0$ in the delocalized regime and at $k_{\sigma}^{\rm max}$ in the localized regime near the transition. Thus we find that the shift in $n_{k\sigma}$ due to presence of a finite $q$ survives in the presence of strong on-site interaction. Similar conclusions can be drawn for weakly interacting bosons for which $U/t\ll 1$supp1 .

Phase diagram at finite $V$: Next we turn on a finite $V$ for the hardcore bosons and obtain the phase diagram by computing SFF and BCF as a function of $V/t$ and $\lambda/t$ for a fixed $q$ and $\Omega$ (see Fig. 3(a,b)). For small $V$ we find that an increase of $\lambda$ leads to a depletion of superfluid density keeping the condensate fraction finite indicating a finite-size crossover from a SF to a localized phase. Similarly for a fixed small $\lambda$, an increase in $V$ leads to an analogous depletion of superfluid density; this indicates the onset of the DW phase with broken translational symmetry. For $\lambda=\Omega=0$, Eq.1 reduces to the well studied XXZ model which exhibits the SF to DW transition at $V/t=1$ XXZ . Similarly for the SF-DW transition at small $\lambda$, a first order transition is expected since the DW state breaks translational invariance while the SF states breaks $U(1)$ gauge symmetry.

The phase diagram of the bosons in the large $V/t$ regime as a function of $\lambda$ can not be completely understood from Fig. 3(a) and (b) since $\rho_{s}=0$ for all $\lambda$ in this regime. To have an understanding of the nature of the boson phase with increasing $\lambda$, we study the structure factor $S(k)=4\sum_{l,l^{\prime}=1}^{N_{s}}e^{ik(l-l^{\prime})}\langle\hat{n}_{l}\hat{n}_{l^{\prime}}\rangle/N_{s}^{2}$. We first note that in the limit $\lambda/t<<1$ the ground state forms a DW leading to $S(\pi)\simeq 1$ and $S(k)\simeq 0$ for $k\neq\pi$ SPati . This DW state is expected to melt with increasing $\lambda$ leading to a vanishing of peak of $S(k)$ at $k=\pi$. A plot of $S(\pi)$ in the $\lambda-V$ plane, shown in Fig. 3(c), indicates the melting with increasing $\lambda$. The dynamical signature of such melting may be obtained by studying boson dynamics following quench of $\lambda$ across its melting value supp1 .

GUE-GOE spectral statistics change: Finally, we show that the present model with $V\neq 0$ hosts a change of spectral statistics from GUE-GOE in the superfluid phase at finite $\lambda$. To this end, we first note that for $\Omega=0$, $[\hat{H},\hat{S}_{z}]=0$ and the boson ground state lies in the $S_{z}^{\rm total}=N_{s}/2$ sector. However, for states within this sector, one does not have TRS since $[\hat{H},\hat{\mathcal{T}}]\neq 0$ for a fixed $S_{z}$ sector. Thus for $\Omega=0$ with a fixed $S_{z}^{\rm max}\neq 0$ sector, one has $[\hat{H},\hat{S}_{z}\hat{\mathcal{T}}]\neq 0$. In contrast for $\Omega\neq 0$, it is easy to see using Eq. 1, $[\hat{H},\hat{S}_{z}\hat{\mathcal{T}}]=0$. The latter symmetry is a consequence of invariance of $\hat{H}$ under TRS followed by a $\pi$ rotation in spin space about the $z$ axis. The presence of this additional symmetry leads to GOE to GUE crossover as $\Omega$ is turned on and increased Haake_p ; evec_dist .

To show this, we first calculate the level spacing ratio quasirev1 ; Bogomolny $r_{\nu}={\rm min}(\delta_{\nu+1},\delta_{\nu})/{\rm max}(\delta_{\nu+1},\delta_{\nu})$, where $\delta_{\nu}=E_{\nu+1}-E_{\nu}$, $E_{\nu}$ being the $\nu^{\rm th}$ energy eigenvalue. We compute the quantity $\langle r\rangle=\sum_{\nu}r_{\nu}/{\mathcal{N}}$, where ${\mathcal{N}}$ is the total number of levels. For $\Omega=0$, working with the energy levels in the maximal $S_{z}$ sector, we find that $\langle r\rangle$ shows a crossover from its GUE value of $\approx 0.58$ to that for Poisson statistics $\approx 0.38$ with increasing $\lambda$ (see Fig. 4(a)). In contrast, for large $\Omega=2$, a similar analysis shows that $\langle r\rangle$ crosses over from its GOE value of $\approx 0.527$ to Poisson with increasing $\lambda$ (see Fig. 4(b)).

In finite-sized systems with no strict symmetry breaking, the level statistics can not be captured for small but finite $\Omega$ values. We therefore concentrate on the variation of the Shannon and structural entropy for studying the crossover between GUE-GOE statistics. The eigenvector corresponding to the $\nu^{\rm th}$ eigenmode can be written as $|\Phi_{\nu}\rangle=\sum_{\chi}c_{\nu}^{\chi}|\chi\rangle$ where $|\chi\rangle$ are the basis states. The corresponding Shannon entropy is given by $S_{\nu}^{\rm Shn}=-\sum_{\chi}|c_{\nu}^{\chi}|^{2}\ln|c_{\nu}^{\chi}|^{2}$. It is well known that $S^{\rm Shn}=\sum_{\nu}S_{\nu}^{\rm Shn}$ has the value $S_{GOE}^{\rm Shn}=\Psi(N/2+1)-\Psi(3/2)$ for GOE and $S_{GUE}^{Shn}=\Psi(N+1)-\Psi(2)$ for GUE Haake_b ; Izrailev_rmt . Here $\Psi$ is the Digamma function and $N$ is the system dimension. The structural entropy for the $\nu$th eigenmode is defined as follows $S_{\nu}^{Str}=S_{\nu}^{Shn}-\ln\xi_{\nu}$ where $\xi_{\nu}$ is the IPR corresponding to the $\nu$th eigenmode. It is known that $S_{\nu}^{Str}\approx 0.37(0.27)$ for GOE(GUE) Varga ; Izrailev_rmt . In Fig. 4(c), we have plotted the distribution of $S_{\nu}^{Str}$ showing that the peak of the distribution shifts from it’s GUE value ($\approx 0.27$) to its GOE value ($\approx 0.37$) as $\Omega$ is changed from $0.6$ to $2.0$. In Fig. 4(d), we plot the variation of $S_{\nu}^{\rm shn}$ showing a smooth crossover from its value for GUE to that for GOE with increasing $\Omega$.

Discussion: Apart from a rich phase diagram, our analysis shows that in a system of 1D ultracold bosons in an optical lattice the interplay between SO interaction, Raman coupling and AA potential leads to novel effects, particularly the splitting in the spin resolved momentum distribution. Above $\Omega_{c}$, such a split happens only when both $\lambda$, $q\neq 0$ and may serve as an indicator of the localization transition. The experimental verification of this splitting would involve preparing a system of bosons with SO coupling nonab1 in the presence of a 1D bichromatic lattice to model AA potential inguscio ; finally the spin-resolved momentum distribution of these bosons can be measured by usual Stern-Gerlach technique exp3 . Our prediction is that such an experiment would observe a spin-split momentum distribution near the localization transition which increases with increasing $\lambda$ or $q$. We note that typically experiments are done with finite lattice sites $N_{s}\sim 12$ greiner1 ; thus our numerical results are expected to be of direct relevance for experimental systems. We have also shown that the spectral statistics of the present model follows Poissonian distribution for large $\lambda$ indicating localization and hosts a GUE-GOE crossover as a function of $\Omega$ in the delocalized regime. Finally we have identified the existence of localized glassy and DW phases as a result of the interaction and the AA potential.

Acknowledgement: BM thanks A. Dutta and S. Mukherjee for discussion.

Appendix A Non-interacting limit

The single particle Hamiltonian of a spin-orbit(SO) coupled bosonic system in an optical lattice (as given by Eq. 1 of the main text with $\lambda=\mathcal{V}=0$) can be written as a $2\times 2$ matrix in the momentum representation as

[TABLE]

The energy dispersion of $H_{\rm SO}$ is given by Eq. 2 of the main text. From this dispersion, which provides the expression of the lower branch of the spectrum $E_{k}^{-}$, we find that there exists a critical value $\Omega_{c}$ above which the ground state is doubly degenerate and the energy minima shifts to finite momenta $k_{0}=\pm\cos^{-1}[\cos q\sqrt{1+\Omega^{2}/(4\sin^{2}q)}]$ (see Fig. 5).

The effective mass (or the band mass) of the bosons is thus given by $m^{*}=\partial^{2}_{k}E_{k}^{-}|_{k=k_{0}}$. For $\Omega>\Omega_{c}$, the expression of $m^{\ast}$ can be written as,

[TABLE]

We note that the effective masses $m^{*}_{>}(m^{*}_{<})$ in the two regimes $\Omega>\Omega_{c}(<\Omega_{c})$ both vanish at $\Omega=\Omega_{c}$. Furthermore, in the absence of disorder the superfluid fraction (SFF) turns out to be the boson effective mass. Thus $m^{\ast}$ captures the behavior of the SFF obtained numerically as a function of $\Omega$ (see Fig. 6(a)); this situation is similar to that obtained in the continuum limit stringari .

In presence of the AA potential we numerically diagonalize the single particle Hamiltonian to obtain the ground state and the full excitation spectrum. The localization transition of the ground state is characterized by the vanishing of superfluid fraction (SFF) which can be obtained using Eq. 4 of the main text. Alternatively, one can also adopt the perturbative approach to calculate the SFF. To this end, we note that in the presence of a small phase twist $\theta$ across the boundary, the original Hamiltonian becomes,

[TABLE]

An expansion of the $H_{\theta}$ (Eq. 9) to ${\rm O}(\theta^{2})$ yields,

[TABLE]

where, $H_{0}$ is the unperturbed Hamiltonian, $\hat{T}=-t\sum_{l,\sigma}(\hat{b}_{l+1,\sigma}^{\dagger}e^{iq\hat{\sigma}_{z}}\hat{b}_{l+1,\sigma}+h.c.)$ is the kinetic energy operator and $\hat{J}=it\sum_{l,\sigma}(\hat{b}_{l+1,\sigma}^{\dagger}e^{iq\hat{\sigma}_{z}}\hat{b}_{l,\sigma}-h.c.)$ is the current operator. So, to $O(\theta^{2})$ the superfluid fraction is given by,

[TABLE]

where [math] and $\nu$ stands for the lowest and $\nu$th eigenmode respectively. In Fig. 6(b) we have plotted $f_{S}$ using the above prescription; we note that the SFF vanishes at $\lambda_{c}<2$ at which the IPR starts rising indicating the localization transition.

The localization transition can also be qualitatively understood from the vanishing of the energy gap $\Delta E$ at the critical disorder strength $\lambda_{c}$. The energy gap ($\Delta E$) between the ground state and the 1st excited state is expected to vanish at the localization transition point. Using this fact, one can obtain a qualitative understanding of the phase diagram for both small and large $\Omega$. We note that for small $\delta\Omega=\Omega-\Omega_{c}>0$ at $\lambda=0$, the ground state is at $k=0$ and has an energy $E(k=0)=-(\Omega+2\cos q)$. Now let us turn on $\lambda$ which leads to a perturbation term, which can be written in momentum space as

[TABLE]

Such a perturbation term leads to a hybridization of the ground state at $k=0$ with the one at $k=\beta$ which has energy $E(k=\beta)=-2\cos\beta\cos q-|2\sin\beta\sin q|+{\rm O}(\Omega^{2})$. Thus the simplest qualitative estimate of the transition line for small $\delta\Omega$ occurs when $\lambda\simeq E(k=\beta)-E(k=0)$ leading to

[TABLE]

We note that this reproduces the linear behavior of the phase boundary for small $\delta\Omega$. A similar analysis can also be carried out at $\Omega\gg 1$. Here the ground state is again at $k=0$ for $q<\pi/2$. An exactly similar analysis as the one charted out above shows that for this case $E[k=\beta]-E[k=0]=2\cos q(1-\cos\beta)+{\rm O}(1/\Omega)$ which leads to

[TABLE]

Thus the phase boundary becomes a horizontal line in the $\lambda-\Omega$ plane, as also seen in exact numerics. In Fig. 6(c) we have shown $\Delta E$ as a function of $\lambda$ for two different system sizes corroborating these qualitative features and justifying the assumption of vanishing of $\Delta E$ at the transition point.

Appendix B Localization of weakly interacting bosons

In the weakly interacting limit, i.e., for $U/t<<1$ and $V=0$, we replace the quantum field operator $\hat{b}_{l,\sigma}$ by the classical field operator $\psi_{l,\sigma}$ assuming the existence of a 1D quasi-condensate shlyapnikov1d . By minimizing the energy functional calculated thereby, we obtain the discrete non-linear Schrödinger (DNLS) equation for the condensate wave function $\psi_{l,\sigma}$ given by,

[TABLE]

where $\mu$ is the chemical potential. We then obtain the ground state wavefunction $\psi_{l\sigma}$ numerically and use it to compute all relevant quantities such as IPR and $f_{s}$. In what follows we have shown the results of such numerical study which are shown in Fig. 7.

In Fig. 7(a) we plot the ground state IPR as a function of $\lambda$ for different interaction strength $UN_{p}/t$. We see that on increasing $\lambda$ beyond the localization transition, the growth of IPR decreases. This is due to the fact that the ground state wavefunction becomes multi-site localized due to weak repulsive interaction (see Fig. 7(c)). We further calculate the superfluid fraction which vanishes in the localized phase as depicted in Fig. 7(b).

To gain a better understanding of the localization transition, we further study the spin resolved momentum distribution of bosons in the regime $\Omega<\Omega_{c}$. In contrast to the non-interacting case, the superfluid with finite $U$ the chooses one of the two symmetry broken states with spins polarized along the z-axis stringari1 . As a result the momentum distribution corresponding to the spin polarization of the ground state becomes highly peaked at the nonvanishing momentum of the ground state as depicted in Fig. 8(a). With increasing disorder strength, other momentum modes get gradually occupied and spin-momentum distributions are peaked at equal and opposite momentum with a net spin polarization indicating symmetry breaking (see Fig. 8). Finally in the localized phase, the momentum distributions become symmetric and peaked around finite momentum with $n_{k,\uparrow}=n_{-k,\downarrow}$. To verify this we plot the order parameter $m=\sum_{k}(n_{k,\uparrow}-n_{-k,\downarrow})^{2}$ and the total magnetization $M=\sum_{k}(n_{k,\uparrow}-n_{k,\downarrow})$ which decreases with increasing $\lambda$ and finally vanishes in the localized phase (see Fig. 7(d)).

Next we investigate the momentum distribution in the regime $\Omega>\Omega_{c}$; similar to the non-interacting case, we see that in the delocalized regime the momentum distributions for both up and down spins are peaked at zero momentum, whereas, in the localized phase they are peaked at finite momentum and other momentum modes get gradually occupied (see Fig. 9(a,b)).

Appendix C Peak splitting and spin dephasing near the localization

transition

This effect of spin-splitted momentum distribution of the localized wavefunction in the regime $\Omega>\Omega_{c}$ arises due to the interplay between the SO interaction and Raman coupling. Here we provide a simple variational calculation to understand this effect. First we consider the variational wavefunction given by

[TABLE]

where $l$ is the site index and $\xi$ represents localization length which is assumed to subsume the effect of AA potential and the interaction. The spinor part is chosen in such a way that up(down) spin momentum distribution is peaked at $\pm k$ and for $k=0$ it reduces to the usual form of the ground state for $\Omega>\Omega_{c}$. For the aforementioned wavefunction, the parameter $k$ is treated as the variational parameter and we investigate its dependence on $\xi$ and $\Omega$. Considering the single particle Hamiltonian of a spin-orbit(SO) coupled bosonic system in an optical lattice (Eq. 1 of the main text with $\lambda=\mathcal{V}=0$), the energy can be written as,

[TABLE]

From the structure of $E(k)$ we note that the $E(k)\to-(\cos(k-q)+\Omega\delta_{k0})+{\rm O}(1/\xi)$ in the delocalized limit where $\xi\to\infty$. This implies that this functional reproduces the correct $k=0$ ground state for $\Omega>\Omega_{c}$ in the absence of the AA potential. Thus in this case the momentum distribution of both the spin-up and spin-down components are peaked at $k=0$. In the strongly localized phase, were $\xi\ll 1$, the second term dominates and in the limit of single site localization $k$ looses its meaning. However, in between these two limits for finite $\xi$, the ground state minima shifts to finite $k$ provided $q\neq 0$. This is seen by minimizing $E(k)$ to obtain $k_{min}$ and by plotting its variation as a function of $\xi^{-1}$ as shown in Fig. 10. As seen from Fig. 10, it is evident that $k_{min}$ decreases with increasing $\xi$ and finally it vanishes in the delocalized regime i.e. $\xi^{-1}\rightarrow 0$. We further notice for a fixed $\xi$ the spin splitting (characterized by $k_{min}$) decreases with decreasing strength of SO interaction ($q$) and eventually vanish for $q=0$. This simple variational calculation elucidates how the combined effect of localization and SO interaction gives rise to the spin splitted momentum distribution in the regime $\Omega>\Omega_{c}$.

Spin-phase diffusion: The effect of spin-splitting in momentum distribution near the localization transition is also accompanied with the phase fluctuation of the wavefunction. In general, the wavefunction can be written as,

[TABLE]

where $\phi_{l}=\phi_{\uparrow}^{l}-\phi_{\downarrow}^{l}$ is the phase angle of the spinor at site $l$. For $\Omega>\Omega_{c}$, we find that $\cos\theta=\sin\theta\approx 1/\sqrt{2}$ and $\phi\approx\pi$ in the delocalized phase, whereas, near localization transition due to increasing phase fluctuations, the phase angle fluctuates significantly from $\pi$ at different sites. We quantify the phase fluctuation by calculating $|\langle e^{i\phi}\rangle|$, where the average is taken over all the lattice sites. In Fig. 11 we have shown the behavior of $|\langle e^{i\phi}\rangle|$ as a function of the disorder potential strength $\lambda$ which shows that near the localization transition it decreases from $1$ with increasing strength of the disorder $\lambda$. In case of hard core bosons we consider the eigenvector corresponding to the largest eigenvalue of the density matrix defined in the main text and calculate the similar quantity.

Appendix D Non-equilibrium dynamics in the strongly interacting regime

To elucidate the localization transition of the HCB, we now look at into the non-equilibrium dynamics of the bosons. We start from the density wave state at $\lambda=0$ denoted by $|\psi(0)\rangle$. Next we quench $\lambda$ to a finite value $\lambda_{f}$ so that the system Hamiltonian after the quench is given by $H[\lambda_{f}]$. Let us denote the eigenfunctions and eigenvalues of $H[\lambda_{f}]$ as $|m\rangle$ and $\epsilon_{m}$ respectively. The time evolved wavefunction $|\psi(t)\rangle$ at any instant of time $t$ after the quench can be obtained by solving the Schrodinger equation $i\hbar|\psi(t)\rangle=H[\lambda_{f}]|\psi(t)\rangle$ and is given by

[TABLE]

The expectation value of any operator $O(t)$ can be obtained from $|\psi(t)\rangle$ as

[TABLE]

Using Eq. 19, we calculate the time evolution of the imbalance factor which is defined as,

[TABLE]

where, $N_{o[e]}=\langle\psi(t)|\sum_{r\in{\rm odd[even]sites}}{\hat{b}}_{i}^{\dagger}{\hat{b}}_{i}|\psi(t)\rangle$ and $N_{tot}=N_{o}+N_{e}$. Note that at $t=0$, we have density wave state with ${\mathcal{I}}=1$ and it approaches zero for a delocalized state. In Fig. 12(a) we have shown the time evolution of $\mathcal{I}(t)$ for the up spin species (the same feature can be observed for the down spin species as well) for different $\lambda$ values. We note that for small $\lambda$ which corresponds to the delocalized regime, $\mathcal{I}$ vanishes to zero with time showing the ergodic dynamics in that regime, whereas for larger $\lambda$ value which corresponds to the localized regime, $\mathcal{I}$ doesn’t vanish and saturates to some positive value which indicates the non-ergodic regime and the density wave ordering is retained in the course of time evolution. In Fig. 12(b) the final density distribution after the time evolution has been shown for different values of $\lambda$.

We repeated the same numerical experiment starting from a different initial state where the atoms are loaded on one half of the lattice and study the imbalance factor $\mathcal{I^{\prime}}=(N_{l}-N_{r})/N_{tot}$ as a function of time $N_{l}$ and $N_{r}$ being the total number density of bosons at the left and the right halves of the lattice respectively. In Fig. 12(c,d) we have plotted the time evolution of the imbalance factor and the final density distribution of the up spin species for different values of the disorder strength.