Thermal transitions of the modulated superfluid for spin-orbit coupled correlated bosons in an optical lattice

TL;DR

This paper explores how spin-orbit coupling influences the thermal phase transitions of correlated bosons in an optical lattice, revealing various superfluid phases and their thermal behaviors.

Contribution

It introduces a comprehensive analysis of thermal transitions in a spin-orbit coupled Bose-Hubbard model, including phase boundaries and experimental implications.

Findings

Spin-orbit coupling promotes finite wavevector condensation.

Thermal broadening affects momentum distribution and magnetic textures.

Critical interactions and temperatures decrease with increasing spin-orbit coupling.

Abstract

We investigate the thermal physics of a Bose-Hubbard model with Rashba spin-orbit coupling starting from a strong coupling mean-field ground state. The essential role of the spin-orbit coupling $\left(\gamma\right)$ is to promote condensation of the bosons at a finite wavevector ${k}_{0}$. We find that the bosons display either homogeneous or phase-twisted or orbital ordered superfluid phases, depending on $\gamma$ and the inter-species interaction strength $\lambda$. We show that an increase of $\gamma$ leads to suppression of the critical interaction $U_c$ for the superfluid to Mott insulator transition in the ground state, and a reduction of the $T_c$ for superfluid to Bose-liquid transition at a fixed interaction. We capture the thermal broadening in the momentum distribution function, and the real space profiles of the thermally disordered magnetic textures, including their homogenization for $T \gtrsim T_{c}$. We provide a Landau theory based description of the ground state phase boundaries and thermal transition scales, and discuss experiments which can test our theory.

Equations

$$\displaystyle H\leavevmode\nobreak\ \leavevmode\nobreak\$$

$$H_{kin}$$

$$\displaystyle H_{U}\leavevmode\nobreak\$$

$$\text{leavevmode}\text{leavevmode}\text{leavevmode}-\sum _{i;\alpha }(\mu +\Omega {\sigma }_z)n_{i\alpha }.$$

$$H_{kin}=\sum _{\mathbf{k}}$$

$$\text{mathbbm}{h}_{\mathbf{k}}=-2t$$

$$\text{leavevmode}\text{leavevmode}\text{leavevmode}+\sin \gamma$$

$$E_{\mathbf{k}}^{\pm }$$

$$\mp (\sin \gamma \sqrt{{\sin }^2{k}_x+{\sin }^2{k}_y})]$$

$${\chi }_{\mathbf{k}}^{\pm }$$

$$Z$$

$$S^{loc}$$

$$\displaystyle+\lambda U\sum\limits_{i}n_{i1}n_{i2}-\sum\limits_{i;\alpha}(\mu+\Omega\sigma_{z})n_{i\alpha}\Big{]}$$

$$S^{hop}$$

$${\psi }_{\mathbf{k}}^{+}$$

$${\psi }_{\mathbf{k}}^{-}$$

$$S$$

$${\tilde{S}}^{hop}$$

$$H^{\mathrm{eff}}$$

$$H_U$$

$$H_{kin}^{\mathrm{eff}}$$

$${\Gamma }_i$$

$${\mathcal{M}}_{ji}$$

$${\Phi }_i=\left(\begin{array}{c}{\varphi }_{k_0}^{+}\\ {\varphi }_{k_0}^{-}\end{array}\right)\exp (\iota {\mathbf{k}}_0.{\mathbf{r}}_i)$$

$${\Phi }_i=\left(\begin{array}{c}{\varphi }_{k_0}^{+}\\ {\varphi }_{k_0}^{-}\end{array}\right)\cos ({\mathbf{k}}_0.{\mathbf{r}}_i)$$

$${\Phi }_i=\left(\begin{array}{c}{\varphi }_{k_0}^{+}\\ {\varphi }_k^{-}\end{array}\right)$$

$$S_{\mathbf{q}}=⟨\frac{1}{zV}\sum _{\mathbf{i},\mathbf{j}}Tr\left[{\Phi }_{\mathbf{i}}^{†}{\Phi }_{\mathbf{j}}\right]e^{\iota \mathbf{q}\cdot (\mathbf{i}-\mathbf{j})}⟩$$

$${\mathbf{m}}_i=⟨\frac{1}{Z}\sum _{\mu ,\nu }Tr\left[e^{-\beta H_{eff}}{b}_{i\mu }^{†}{\mathbf{\sigma }}_{\mu \nu }{b}_{i\nu }\right]⟩$$

$$n_{\mathbf{k}}=\frac{1}{N}⟨\frac{1}{ZV}\sum _{\mathbf{i},\mathbf{j},\mu }Tr\left[e^{-\beta H_{eff}}{b}_{\mathbf{i}\mu }^{†}{b}_{\mathbf{j}\mu }\right]e^{\iota \mathbf{k}\cdot (\mathbf{i}-\mathbf{j})}⟩$$

$$S^{hop}=S^{neg}+S^{pos}$$

$$S^{neg}$$

$$S^{pos}$$

$$e^{-S^{neg}}$$

$$=\prod _{k\sigma n}\{\int \mathcal{D}\left[{\varphi }_{k\sigma n}^{\ast },{\varphi }_{k\sigma n}\right]e^{{\varphi }_{k\sigma n}^{\ast }{\tilde{E}}_{k\sigma }^{-1}{\varphi }_{k\sigma n}}$$

$$\times e^{-\left({\psi }_{k\sigma n}^{\ast }{\varphi }_{k\sigma n}+{\varphi }_{k\sigma n}^{\ast }{\psi }_{k\sigma n}\right)}\}$$

$$\stackrel{\varphi \to \sqrt{-\tilde{E}}\varphi }{=}$$

$$\times e^{-\sqrt{-{\tilde{E}}_{k\sigma }}\left({\psi }_{k\sigma n}^{\ast }{\varphi }_{k\sigma n}+{\varphi }_{k\sigma n}^{\ast }{\psi }_{k\sigma n}\right)}\}$$

$$\displaystyle\delta E^{(2)}=-\frac{1}{2U\tilde{\mu}}\sum\limits_{i}\biggl{[}f\left(\tilde{\mu},\lambda\right)\left(|\Gamma_{i1}|^{2}+|\Gamma_{i2}|^{2}\right)$$

$$\displaystyle+\sqrt{\biggl{(}g\left(\tilde{\mu},\lambda\right)\left(|\Gamma_{i1}|^{2}-|\Gamma_{i2}|^{2}\right)\biggr{)}^{2}+\left(\frac{2\lambda\tilde{\mu}|\Gamma_{i1}\Gamma_{i2}|}{\lambda-\tilde{\mu}}\right)^{2}}\,\,\biggr{]}$$

$$+\sum _{i\sigma }|{\varphi }_{i\sigma }{|}^2$$

$$f\left(\tilde{\mu },\lambda \right)$$

$$g\left(\tilde{\mu },\lambda \right)$$

$$|{\Gamma }_{\mathbf{i}\alpha }{|}^2=\sum _{\mathbf{j}\sigma ;\mathbf{l}\delta ;\mathbf{k},\mathbf{q}}$$

Full text

spin-orbit coupling, Bose-Hubbard model

I Introduction

The physics of strong correlation in ultracold atom systems has been a subject of intense theoretical and experimental research in the recent past rev1 ; greiner1 ; expt1 ; expt2 ; bht1 ; bht2 . The initial studies in this field concentrated on single boson species. This choice is motivated by the experimental ease of realizing the superfluid (SF) and Mott insulating (MI) states of these bosons. Indeed, the first experimental study of SF-MI quantum phase transition used ${}^{87}{\rm Rb}$ bosons in their $F=1$ state greiner1 . More recently, there have been concrete proposals to realize artificial Abelian gauge fields for such bosons gauge1 ; gauge2 . The phase diagram of strongly correlated bosons in the presence of such gauge fields have also been investigated abeth1 ; abeth2 and reveal a rich structure.

Several recent cold atomic experiments tune Raman processes to create artificial spin-orbit couplings in multicomponent Bose systems spielman2011 ; ketterle2016 ; ketterle2017 . Most of these experimental procedures produce an equal mixture of Rashbha and Dresselhaus coupling, which leads to an effective Abelian gauge field for the bosons. However, there have been concrete proposals to experimentally realize purely Rashba type spin-orbit coupling spielman-review . This is equivalent to a non-Abelian gauge-field for two component bosons.

The ground state phase diagram of such systems have been theoretically studied iskin ; nandini-prl ; hofstetter ; kush1 ; saptarshi1 . These studies employed several theoretical techniques such as mean field theoriesiskin , simulated annealing of effective quantum spin modelsnandini-prl , real space bosonic dynamical mean field theory (BDMFT) hofstetter , and strong coupling expansion kush1 ; saptarshi1 . They have unearthed a rich ground state phase diagram for these systems. Some of the unconventional phases found include those with long range magnetic order in the Mott ground state nandini-prl and the possibility of a boson condensate at finite momentum kush1 ; saptarshi1 . Such studies have also been supplemented by their weak-coupling counterparts in the continuum where there is no Mott transition. The weakly interacting condensates have been studied using the Bogoliubov-Hartree-Fock approximation baym2014 .

In spite of several studies on the ground state, only limited theoretical work exists on the thermal phases of spin-orbit coupled systems. For Abelian systems with equal mixture of Rashba and Dresselhaus coupling, Ref. hickey, derives an effective $t-J$ model for the bosons and studies the thermal phases of this effective model. The study reveals a stripe superfluid order at low temperature and a two step melting upon increasing temperature, leading first to a striped normal phase of the bosons and then to a homogeneous state. Similar studies were carried out for two component fermions in optical lattices ref51 . However, to the best of our knowledge, the thermal phases of Bose-Einstein condensates (BECs) in the presence of Rashba spin-orbit coupling have not been studied before. This is particularly pertinent since an equal mixture of Rashba and Dresselhaus terms breaks the four-fold rotation symmetry of the lattice, while the Rashba spin-orbit term keeps it intact. This leads to the possibility of superfluid phases with lower symmetry than that of the lattice.

In this work, we study the thermal phases of a two-orbital Bose-Hubbard model in the presence of a Rashba spin-orbit coupling. Our study thus involves bosons in the presence of an effective non-Abelian gauge field. In what follows, we use an auxiliary field decomposition of the kinetic energy followed by a ‘classical’ approximation to the auxiliary field. We then carry out a Monte-Carlo study of the resulting model, sampling the auxiliary field configurations. The method has been used in the past for the single species Bose-Hubbard model joshi-thermal . It retains the key low energy thermal fluctuations and yields accurate thermal transition scales.

We start by deriving an effective Hamiltonian whose mean field ground state coincides, in the main, with earlier results nandini-prl . Our results on this problem are the following: (i) We find that the ground state is either a Mott insulator, or a superfluid with condensation either at a single wavevector $({\bf k}_{0})$ or two wavevectors $(\pm{\bf k}_{0})$. The $\pm{\bf k}_{0}$ condensate constitutes a orbital density wave, while the finite ${\bf k}_{0}$ condensate is a phase twisted superfluid saptarshi1 . (ii) The superfluid has associated ‘magnetic’ textures - related to the spatially varying orbital occupancy. (iii) Increasing temperature leads to the simultaneous loss of superfluidity and order in the magnetic textures. We establish the $T_{c}$ scale for varying Hubbard interaction, interspecies coupling and spin-orbit interaction using our Monte Carlo scheme joshi-thermal . (iv) The momentum distribution function, $n_{\bm{k}}$, evolves from its ‘low symmetry’ character at low temperature to four-fold symmetry as $T\rightarrow T_{c}$, providing a detectable thermal signature of Rashba coupling. Finally, (v) we construct an effective Landau theory which provides some analytic understanding of the thermal scales, and discuss experiments which can test our theory.

The plan of the rest of this work is as follows. In Sec. II, we introduce the Bose-Hubbard Hamiltonian in the presence of Rashba spin-orbit coupling and describe the method used for our calculation. This is followed by Sec. III, where we study the ground state phase diagram. We study the finite temperature effect on different phases in Sec. IV. Finally, we discuss our main results, chart out experiments which can test our theory, and conclude in Sec. V. Some details of our calculation and the construction of the effective Landau theory are presented in the Appendices.

II Model and Method

In this section, we shall present the model we use and also discuss the details of the method used for computation.

II.1 Model

We begin by defining a Rashba spin-orbit coupled two-orbital Bose-Hubbard Hamiltonian on a square lattice in 2D:

[TABLE]

Here $\mathcal{R}(i,j)=-t\exp[\iota\bm{A}.(\bm{i}-\bm{j})]/2$ is the real space hopping matrix, $\bm{A}=(\gamma\sigma_{y},-\gamma\sigma_{x},0)$ is the synthetic gauge field. $U$ is the on-site repulsion, $\lambda$ denotes the ratio between inter-orbital and intra-orbital on-site repulsion, and $\Omega$ is the Zeeman field which arises due to the coupling of the Raman laser to the bosonic atom spielman2011 . This term depends on the strength of the atom-laser coupling and can be tuned to the extent that the spin-orbit physics does not get completely masked. In this work, following Refs. nandini-prl, , we shall later set $\Omega$ to zero in order to have a clean demonstration of the effects of spin-orbit coupling. In what follows, we also neglect another additional on-site term $H\sim\delta\sigma_{y}/2$ which depends on the detuning parameter $\delta$ of the Raman laser and can be made small by sufficient reduction of the detuning. For the rest of this work, we set the lattice spacing $a_{0}=1$.

The kinetic part $H_{kin}$ can be mode separated and can be written as

[TABLE]

Here $+(-)$ denotes the upper(lower) bands in Fig.1. The band structure respects $\pi/2$ rotational symmetry of the square lattice. Since the local interaction terms do not break this symmetry, this degeneracy should remain intact even in the many-body spectrum. For Rashba type spin-orbit coupling the band minima always lie on the diagonals of the two-dimensional (2D) Brillouin zone (BZ). The locations are at $(\pm k_{0},\pm k_{0})$ where $k_{0}$ is determined by the strength of the SO coupling: $k_{0}=\tan^{-1}[\tan(\gamma)/\sqrt{2}]$. The noninteracting density of states (DOS) has been shown in Fig.1. As the spin-orbit coupling strength $\gamma$ is varied from [math] to $1/2$, the DOS develops additional van Hove singularities at finite energies, while the singular peak at $\omega=0$ turns into a dip with a linear rise.

II.2 Effective Hamiltonian

In order to simulate the finite temperature physics of this model we introduce auxiliary fields and implement an approximation that maintains a positive definite stiffness for these fields. The usual mean-field decomposition bht1 of the kinetic term does not meet this requirement.

We start by writing the imaginary time coherent state path integral using the Hamiltonian above bht2

[TABLE]

Next, we wish to implement a Hubbard-Stratonovich decomposition of the hopping part of the action. To this end, we segregate the negative energy part of the bands ($\tilde{E}_{\bm{k}}^{\pm}$), and introduce an auxiliary field decomposition of the negative-band action using two fields $\{\phi^{+}_{i,n}\}$, $\{\phi^{-}_{i,n}\}$ for each lattice point and Matsubara frequency, $(i,n)$. The effects of positive energy part of the bands can be built back perturbatively, and should not affect the low-energy physics significantlyjoshi-thermal . The resulting action is given by

[TABLE]

Next, we note that an effective Hamiltonian can be derived from Eq. 4 if we retain only the zero Matsubara frequency mode of the auxiliary fields $\{\phi^{+}_{i,0}\}$, $\{\phi^{-}_{i,0}\}$. For the single orbital problem this approximation reproduces the mean-field sheshadri ground state exactly, and captures thermal scales which agree well with full quantum Monte-Carlo joshi-thermal . The effects of the finite-frequency modes can be built back perturbatively as quantum corrections over the static background. This has been accomplished for the single orbital problem joshi-spectral and such corrections are known to leave the qualitative nature of the thermal phase and phase transitions unchanged. For bosons coupled via spin-orbit coupling, this turns out to be more cumbersome and we defer computation of such corrections to a future work.

The effective Hamiltonian obtained by retaining only $\{\phi^{+}_{i,0}\}$, $\{\phi^{-}_{i,0}\}$ fields is given by

[TABLE]

where $\Phi_{i}\equiv\left(\begin{smallmatrix}\phi^{+}_{i}\\ \phi^{-}_{i}\end{smallmatrix}\right)$ is a local spinor composed of zero mode of the auxiliary fields $\{\phi^{+}_{0}\}\equiv\{\phi^{+}\}$ and $\{\phi^{-}_{0}\}\equiv\{\phi^{-}\}$. $\Psi_{i}\equiv\left(\begin{smallmatrix}b^{1}_{i}\\ b^{2}_{i}\end{smallmatrix}\right)$ is a local spinor involving the bosons in the two orbitals. $\mathcal{M}_{ji}$ are 2x2 matrices which couple the chiral auxiliary fields with the orbital bosonic fields, with coefficients picked up in the band truncation process. The information of the spin-orbit coupling enters the effective Hamiltonian through these coefficient matrices. Here $H_{U}$ is the local interaction part as in the original Hamiltonian 1b and $\Omega$ has been set to zero in the subsequent calculations. The details of the procedure leading to $H_{\rm eff}$ can be found in the AppendixA.

II.3 Methods

The effective Hamiltonian obtained in the last section, can be treated using several approximation schemes. In this work, we are going to use two such schemes. The first of these, used to obtain zero temperature phases of the system, involves treating $\{\Phi_{i}\}$ as variational parameters and subsequent minimization of the energy obtained from the effective Hamiltonian. In this scheme, the energy for a configuration of $\Phi$s is obtained by diagonalizing the boson Hamiltonian $H_{eff}[\Phi_{i}]$. This yields the optimal ground state configuration of $\Phi_{i}$ fields. In this work, we restrict ourselves to four families of such variational wavefunctions given by

Single mode:

[TABLE] 2. 2.

Two mode:

[TABLE] 3. 3.

Four mode:

[TABLE] 4. 4.

Vortex:

[TABLE]

where ${\bf r}_{i}=(x_{i},y_{i})$ are the coordinates of site $i$. A sketch of these variational profiles of $\Phi_{i}$ and the corresponding magnetic texture of the bosons is given in Fig. 4. We note that the local Hilbert space for the bosons needs to be restricted for the problem to be numerically tractable. This is done by choosing a cutoff, $N_{i}$, in number of boson occupation per site. In what follows, we have ensured that the cut-off is chosen such that including more states beyond it does not have any effect on the energy of the system, up to a desired accuracy. The variational calculation gives us the mean field ground state of our effective model5.

Having obtained the ground state configuration of the bosons, the second method we use yields information about its thermal behavior. To this end, we use a classical Monte-Carlo scheme by starting from the ground state configuration and successively increasing the temperature. The free energy for a configuration of $\{\Phi_{i}\}$s is again obtained by diagonalizing the boson Hamiltonian $H_{eff}$ for every attempted update of the auxiliary fields. The equilibrium $\{\Phi_{i}\}$ configurations are generated by implementing a Metropolis based update scheme. In this scheme, at any given site $i$, we have two complex scalar auxiliary fields, $\phi^{+}_{i}$ and $\phi^{-}_{i}$. For each of the fields, the amplitude fluctuations are considered to be within twice their ground state amplitude. In contrast, arbitrary phase fluctuations of these fields are allowed. The local hybridization $\Gamma_{i}$ depends on the $\Phi_{i}$ configurations on all sites, as defined in equation II.2. For a given $\{\Phi_{i}\}$ configuration the bosonic Hamiltonian is written in Fock space after truncating the local Hilbert space within $N_{i}$ particle states, as in the variational calculation. The resulting matrix is then diagonalized exactly to obtain the free energy for the configuration.

II.4 Indicators

To detect the presence of spatial order we compute the structure factor:

[TABLE]

where $V$ is the volume of the system, $z$ is the coordination number and $\Phi_{i}$s are the auxiliary fields introduced in sec.II.1.

The local magnetic texture of the two-orbital bosons is defined by the vector,

[TABLE]

where $Z$ is the partition function and the angular brackets denote thermal averaging.

The momentum distribution of the bosons given by:

[TABLE]

where $N$ is the total no. of bosons, $Z$ is the partition function, $V$ is system volume, and the angular brackets denote thermal averaging.

III Variational Ground State

In this section, we shall use the variational scheme outlined earlier to obtain the mean-field ground state phase diagram of the bosons. In what follows, we have numerically implemented this scheme on a $16\,\times\,16$ unit cell with $4\leq N_{i}\leq 10$ hybridization states per site. The chosen value of $N_{i}$ depending on the value of the on-site interaction $U$. For every parameter point $N_{i}$ have been fixed at its optimal value, so that increasing it does not affect the results. Unless otherwise mentioned, the filling should be considered as fixed to one boson per site.

Due to the symmetry in the problem, we can restrict $\gamma$ to the interval $[0,0.5]$. Moreover, we notice that in the atomic limit, where the problem becomes independent of $\gamma$, the level schemes differ qualitatively if one tunes $\lambda$ across unity, as shown in Appendix B (see Fig. 11). This allows us to segregate the two parameter regimes - $\bm{\lambda<1}$ and $\bm{\lambda>1}$. We present our results for a characteristic value of $\lambda$ in each of these intervals ($\lambda=0.5$ and $1.5$ respectively), and expect qualitatively similar trends for other values of $\lambda$ in the respective intervals. At each parameter point we first classify the ground state phases using expectation values of linear bosonic operators like $\langle b^{\dagger}_{i\alpha}\rangle$. This allows us to demarcate the ground state superfluid (SF) - Mott insulator (MI) phase boundary (FIG. 3). The order parameter vanishes in the MI phase, as a result, the kinetic part of the Hamiltonian has no contribution in the energy and we recover the atomic limit. In the SF phase a non-vanishing amplitude of $\langle b^{\dagger}_{i\alpha}\rangle$ survives throughout the system, while in the MI phase it vanishes on all sites. We further classify the superfluid phases by using expectation values of bosonic bilinears as defined in Eq. 9. This yields a classification of the superfluid phases into the following subcategories:

• •

Homogeneous - where $\langle b^{\dagger}_{i\alpha}\rangle$ and the bilinears remain constant throughout the system.

• •

Phase-twisted - where the amplitude of $\langle b^{\dagger}_{i\alpha}\rangle$ as well as the bilinears remain constant throughout the system, but the phase of $\langle b^{\dagger}_{i\alpha}\rangle$ varies from site to site.

• •

Z-FM - in which $\langle b^{\dagger}_{i1}\rangle$ retains a homogeneous nonzero value, but $\langle b^{\dagger}_{i2}\rangle$ vanishes throughout the system; $m_{z}$ remains pinned to 1, while $m_{x}$ and $m_{y}$ vanish.

• •

Stripe - in which both the amplitude as well as the phase of $\langle b^{\dagger}_{i\alpha}\rangle$ vary from site to site, and the bilinears show stripe like patterns across the system.

The effect of increasing $\gamma$ at fixed $U$ and $\lambda$ can be understood as follows. The effective bandwidth of the system varies with $\gamma$ as $W\left(\gamma\right)=4t\sqrt{2\left(1+\cos^{2}\gamma\right)}$. Thus one requires progressively larger bare hopping $t/U$ to compensate for the $\cos^{2}\gamma$ term in order to stabilize the superfluid phase. Thus we expect $t_{c}$ to increase with $\gamma$ for fixed $U$ and $\lambda$. This expectation is verified in our numerics as can be seen from both panels of Fig. 3. Within the superfluid phase, the phase diagram can be broadly classified into three separate regimes. In the first of these, where $\bm{\lambda<1}$ (FIG.3(a)), single mode variational profile minimizes $H_{\rm eff}$. For any finite value of $\gamma$ this leads to a phase-twisted superfluid with uniform density in both the orbitals throughout the system, while for $\gamma=0$ it reduces to the homogeneous superfluid phase. The fact that any finite $\gamma$ would necessarily lead to a phase twisted superfluid can be understood in terms of an effective Landau functional, which has been discussed in Appendix B.

In the second regime where $\bm{\lambda>1}$ (FIG. 3(b)), for low values of $\gamma$ we get condensation in only one of the orbitals, leading to a $z$-polarized ferromagnetic texture as shown in . In contrast, for larger values of $\gamma$, the two mode variational state wins over others in the superfluid phase, leading to a stripe-like orbital order with modulating density in each orbital. The pitch of the orbital density wave depends of $\gamma$, and for $\gamma>0.4$ it leads to a Z-AFM order. The complete phase diagram in the superfluid phase as a function of $\gamma$ and $\lambda$ is shown in Fig. 4. The superfluid-Mott phase boundary is governed by the vanishing of the second order coefficient of the Landau functional obtained by tracing out the bosons in the strong coupling limit. We discuss this procedure in detail and chart out the analytic intuition obtained from it in AppendixB.

We note here that in our calculations we find that the four mode and vortex configurations do not feature in the ground state, although at certain parameter points their energies come very close to the ground state energy. This is in contrast to the phase diagram obtained in previous works nandini-prl ; hofstetter using other techniques. This might be an artifact of band truncation in our implementation of the mean-field approximation, although it is not entirely clear whether other mean-field approaches can actually capture those phases iskin . Nevertheless, at larger values of $\lambda\left(\gtrsim 1.5\right)$, our ground state phase diagram matches qualitatively with that in Ref. nandini-prl, . In this region, we wish to highlight our finite temperature results, since the merit of our technique is in capturing the thermal scales nonperturbatively, which could not have been possible, to this extent, using other techniques.

Next, we study the magnetic structure of the ground state. The magnetic texture, shown in Fig. 5 arises from the relative boson density modulation between the two orbitals over different lattice sites. We find that in the ground state, for $\lambda<1$, $m_{zi}=0$ which indicates that there is no local population imbalance between the two orbitals throughout the lattice as shown in Fig. 5(a). The planar components, which encapsulate the relative phase between the two orbitals, are also same on all sites. In contrast, for $\lambda>1$, the ground state, for $\gamma=0$, has $|m_{zi}|=1$ which means that the bosons condense in only one of the orbitals and the density in the other orbital remains zero on all sites. Increasing $\gamma$ leads to a diagonal stripe-like order with $|m_{zi}|<1$ indicating population imbalance between the two orbitals. This imbalance varies in space leading to the stripe-like order as shown in Fig. 5(b).

At $T=0$ and in the superfluid phase, $n_{\bm{k}}$ is sharply peaked as shown in Fig. 6. The peak height represents the condensate fraction, which depends on the strength of interaction $U$ and the spin-orbit coupling $\gamma$. The condensate gets depleted with increasing $U$ (keeping $\gamma$ and $\lambda$ fixed) leading to diminished peak height. For $\lambda<1$, the position of the momentum distribution peak shifts from $k=0$ to $\left(k_{0},k_{0}\right)$ where $k_{0}$ is given by the band minima. This is shown in the top panel of Fig. 6. Note that the position of this minima is controlled by the spin-orbit coupling. For $\lambda>1$ the single peak at $\gamma=0$ splits into two peaks at $\left(\pm k_{0},\pm k_{0}\right)$ with equal heights as shown in the bottom panel of Fig. 6. This indicates that the ground state is a superposition of Bose condensates at two distinct wavevectors. The peak heights diminish with increasing $\gamma$, keeping U fixed. This can be attributed to the fact that the band stiffness about the minimum decreases as the spin-orbit strength is increased. We note that such a superposition state may be unstable in the presence of a trap potential and we shall not address this issue further here.

IV Finite Temperature Results

In this section we chart out the finite temperature phases starting from the variational mean-field ground states obtained in the previous section. We use the classical Monte Carlo scheme described in Sec. II.3 and run the simulation on a 16$\times$16 lattice with two fluctuating fields, $\phi^{\pm i}$ and $\phi^{-}$ at each site $i$. Both the amplitude and the phase interval of the fields are discretized in hundred subintervals. The amplitude interval is restricted to twice the saddle point value while full phase fluctuation has been allowed. The real space $\{\Phi_{i}\}$ configurations are obtained by sampling over four thousand MC sweeps for each temperature. In each these sweeps, all the sites of the system are updated once. A total of $N_{0}=100$ configurations are saved at every temperature, which are subsequently used to calculate thermal averages of observables.

The finite temperature phase diagram is shown below in Fig.7. The low temperature state is the variational ground state which we have discussed at length in Sec. III. As we heat up the system it gets thermally disordered and finally makes transition to a normal state. The normal state is a Bose liquid with no long range order, but short range spatial correlations. The critical temperature $T_{c}$ varies non-monotonically with $U$. As $U$ is lowered stating from $U_{c}$, $T_{c}$ grows linearly up to quite low values of U ($\sim 2-6$ depending on $\lambda$ and $\gamma$) after which it falls suddenly. For $\gamma=0$ the fall is sharp and is easily discernible in Fig. 7, while for finite $\gamma$, it is quite gradual. The low $U$ part of the phase diagram is numerically inaccessible due to large number fluctuations in the condensate, for which one needs to retain enormously high number of local hybridization states. For this reason we could access results only up to $U/t=2$ . With increasing $\gamma$ the $T_{c}$ scales get suppressed at all values of $\lambda$ and U. This can again be attributed to suppression of effective bandwidth by the spin orbit coupling as discussed in Sec. III.

Next, we address the effect of finite temperature on the momentum distribution functions. The results are shown in Fig. 8. The peaks in the ground state momentum distribution show significant thermal broadening with increasing $T$. This is best appreciated by looking at the $\gamma=0$ behavior (top panel in Fig. 8). The condensate fraction remains almost constant up to T = 0.1$T_{c}$, after which particles start getting excited out of the condensate. For $T\simeq T_{c}$ there is significant broadening of the peak even though the superfluid order still survives. Beyond $T_{c}$ phase fluctuations destroy the coherence giving uniform Bose liquid. For finite $\gamma$ one can notice thermal weights developing in the symmetry related $k$-points when the system is close to $T_{c}$, for both the $\lambda$ values. These weights signify the presence of low energy states at certain $k-$ points, which is reminiscent of the band structure symmetry. At temperatures close to $T_{c}$ thermal fluctuations excite particles out of the condensate to these low energy states, without destroying the overall phase coherence in the system. As the system is heated up further the populations in these symmetry related $k-$points tend to homogenize at the cost of destroying superfluidity.

Next, we consider the behavior of the magnetic texture as a function of temperature. As the system is heated from the ground state the magnetic textures start getting disordered. The thermal behavior of the magnetic texture is shown in Fig. 9. We observe that for a temperature $T<T_{c}$ the planar moments become more disordered as compared to $m_{z}$ (shown in color). This can be attributed to the fact that the planar moments capture the gapless phase fluctuations of the superfluid, whereas $m_{z}$ captures their population difference. Finally, for $T>T_{c}$, we find that the planar moments become completely disordered while the $z$ component homogenizes.

We track the peak in the structure factor $S_{\bf q}$ with temperature to locate the onset of long range order as shown in Fig. 10. We find that as the system is heated from its ground state, the auxiliary fields start fluctuating about their saddle point; consequently, the distribution of the $\{\Phi_{i}\}$s broaden. At each site the two variables (per species), i.e.  the amplitude and the phase of the auxiliary field $\phi_{i\sigma}$ get disordered with temperature. It is the fluctuations of the phase degree of freedom which ultimately kill superfluidity in the system. The transition temperature $T_{c}\left(\lambda,\gamma,U\right)$ can be inferred from the ”knee” of the $S_{\bm{q}}$ peak vs temperature curve. Thus this measurements allow us to locate $T_{c}$ which may be relevant in realistic experiments.

V Discussion

In this work we have studied the thermal phases and phase transitions for bosons with Rashba spin-orbit coupling. Our starting point has been a strong coupling mean-field phase of these bosons in the SF phase near the SF-MI critical point. We find that the result of our mean-field study lead to homogeneous, phase-twisted, and orbital density-wave ordered SF phases depending on the strength of spin-orbit coupling. The phase diagram that we find agrees qualitatively with earlier studies using more sophisticated methods nandini-prl . Using these ground states as the starting point, we then perform a finite temperature Monte Carlo study of the thermal properties of the bosons. The thermal phase diagram for the bosons shows reduction of the critical temperature $T_{c}$ with increasing strength of the spin orbit coupling $\gamma$ at a fixed value of the Hubbard interaction $U$. This can be interpreted as spin-orbit coupling introducing an effective frustration in the system leading to reduction of order parameter stiffness and hence $T_{c}$. We also obtain the thermal broadening in the momentum distribution and the presence of satellite peaks at the band minima which reflects the four-fold symmetry of the Rashba term. We note that such four-fold symmetric momentum distribution would be absent in earlier studies which studies an effective Abelian theory involving an equal mixture of Dresselhaus and Rashba spin-orbit terms. We find that the orbital density waves survive to temperatures close to $T_{c}$. Finally, we also study the magnetic textures of these bosons via computation of the magnetization ${\bf m}_{i}$. In particular, we provide a clear description of the thermal evolution of these textures and their subsequent homogenization for $T>T_{c}$.

The present study neglects the quantum fluctuations of the auxiliary fields completely. This leads to an overestimation of $U_{c}$ on one hand, but more importantly, leads to loss of any dynamics in the Mott phase at zero temperature. A scheme for building back the finite frequency quantum modes already exists, and has been used to capture quantum dynamics in the single orbital problemjoshi-thermal . Using that method, in this problem one hopes to recover the vortex-like magnetic textures close to the Mott phasenandini-prl . We leave this issue as a subject of future study.

The simplest experimental verification of our work would be measurement of the momentum distribution of the bosons in the SF phase at finite temperature. We provide a detailed thermal broadening of the momentum distribution function which could be verified by standard experiments. In addition, we also predict that $n_{\bf k}$ would reflect the four-fold symmetry of the Rashba coupling term at finite temperature. This property involves peak positions of the momentum distribution which is easily measured in standard experiments.

Conclusion: We have studied strongly correlated two-component bosons on a square 2D lattice in the presence of Rashba spin-orbit coupling. We focus on the finite temperature problem and use a recently developed auxiliary field based Monte Carlo tool, that retains all the classical thermal fluctuations in this correlated system, to address the thermal state. We establish, to the best of our knowledge for the first time, the superfluid critical temperature $T_{c}$ for varying intra- and inter-species repulsion and spin orbit coupling. We study the momentum distribution and ‘magnetic textures’ as the temperature is increased through $T_{c}$ and highlight the loss of coherence and spatial order. We have predicted experimentally verifiable signatures of the Rashba coupling in the finite temperature superfluid.

We acknowledge use of the HPC clusters at HRI.

Appendix A Derivation of effective action

The full partition function is defined in Eq.3. Keeping $S^{loc}$ intact we wish to decompose the $S^{hop}$ by a Hubbard-Stratonovich (HS) transformation. In order to implement it we need to segregate the negative part of the bands, so that the bosonic Gaussian integral remains well defined. This leads to

[TABLE]

In this work, we neglect the $S^{pos}$ part and implement a HS transformation on the $S^{neg}$.

[TABLE]

where $\{\phi^{+}_{n}\}$ and $\{\phi^{-}_{n}\}$ are the auxiliary fields which couple with the respective chiral bosonic modes. This procedure therefore leads to Eq. 4 of the main text.

Appendix B Landau functional close to $U_{c}$

We derive an effective spin model for the bosons in the SF phase near the SF-MI transition. To this end, note that at large $U/t$, close to the Mott phase, the original boson fields can be integrated out to give an effective description of the bosons in terms of the auxiliary fields. It leads to a Landau energy functional, with coefficients depending on the parameters of the theory. This procedure is similar in spirit to well-known derivation of such effective spin models in the Mott phases of the bosonsdemler1 ; issacson1 ; however, here we obtain such a model for their SF phase.

For the single orbital problem one can derive the free energy functional by performing a cumulant expansion of the SPA functional joshi-thermal . In the two-orbital problem the ground state in the atomic limit is degenerate as shown in Fig. 11). Thus one needs to use degenerate perturbation theory about the atomic limit. The Landau energy functional after second order correction in $\{\Gamma_{i\alpha}\}$ is given by:

[TABLE]

with $\tilde{\mu}\equiv\frac{\mu}{U}$, and

[TABLE]

Notice that the square root term lifts the degeneracy of the ground state. We now express the hybridization fields $\{\Gamma_{i\alpha}\}$ in terms of the auxiliary fields $\{\phi_{i\sigma}\}$ using Eq. II.2.

[TABLE]

If we choose the $\{\phi_{i\sigma}\}$ from the single mode variational family and use the fact that the amplitude for the $\{\phi_{+}\}$ field vanishes in the ground state, then the energy functional can be written as:

[TABLE]

where $V$ is the volume of the system. The condensation wavevector in the ground state is given by the $\bm{k}_{0}$ for which $\alpha^{(2)}$ becomes maximally negative. In the expression of $\alpha^{(2)}$ the factor in brackets remains positive definite for the region of parameter space in which the single mode solution dominates. Hence, the maximally negative value of $\alpha^{(2)}$ occurs at the minima of the lower band, which are given by $\left(\pm k_{0},\pm k_{0}\right)$, with $k_{0}=\tan^{-1}\left[\tan\gamma/\sqrt{2}\right]$. From this, we can also conclude that the presence of an arbitrarily small $\gamma$ would lead to a phase-twisted superfluid. At the optimal $\bm{k}_{0}$, the SF-Mott phase boundary is determined by the zeros of $\alpha^{(2)}$. At $\lambda=0.5$, for which the single mode variational state dominates, we have matched the phase boundary obtained through numerical minimization, with that obtained from the effective Landau theory. We find excellent agreement between the two, as is evident in Fig. 3. A similar match was also found for $\lambda=1.5$ where we have stripe and z-FM like order in the ground state.

Notice that at this level we have truncated the Landau expansion to second order. The energy functional obtained above is quadratic in $\{|\phi_{i\sigma}|\}$, and hence the amplitudes would vanish at the minimum. So, unless we compute the $\delta E^{(4)}$ correction, this scheme cannot be used to optimize over the amplitudes. However, once the optimal amplitudes are fixed from the variational calculation, this functional may be used to anneal the phase of the auxiliary fields, assuming that the amplitude variation with temperature is small close to $U_{c}$. This would allow us to compare the $T_{c}\left(U\right)$ curves of the bosonic theory with the effective spin model. The expectation is that they would coincide at strong coupling, as in Ref. joshi-thermal, , allowing us to describe the physics in terms of the low energy degrees of freedom. For a crude estimate, one can ignore the terms inside the square root to derive a more explicit looking functional in terms of the phase degrees of freedom.

[TABLE]

The couplings $\mathcal{A}$ and $\mathcal{B}$ depend on the band structure, and rapidly decay to zero with increasing distance. This allows us to approximate the lattice sum by just the sum over nearest neighbors (or the next-nearest neighbors, in case the nearest neighbor coupling vanishes). Hence, under all these assumptions, one can extract an effective exchange scale which would allow us to calculate an effective ordering temperature ($T^{*}$) for each point in our parameter space. A comparison of $T^{*}$ with the $T_{c}$ obtained from the monte-carlo has been shown in Fig. 12. The approximation gets better at lower $\gamma$ (where neglecting the terms within the square root in Eq. B can be easily justified) as expected. The match seems reasonably good, given the drastic nature of approximations made for extracting a $T^{*}$ out of the effective Landau functional.