Quantum phases of spin-orbital-angular-momentum coupled bosonic gases in optical lattices

Before exploring the whole model, described by Eq. (6), we first discuss the competition between conventional nearest-neighbor hopping $t$ and orbital-angular-momentum-induced hopping $t_{soc}$. As shown in Fig. 2, the orbital-angular-momentum-induced hopping can be the order of the conventional one even for $l=1$, where the hopping amplitudes are obtained from band-structure simulations ban . By neglecting the Raman-induced spin-flip hopping $\Omega$, Eq. (6) is reduced to

	$\displaystyle H$	$\displaystyle=$	$\displaystyle-\sum_{\langle i,j\rangle,\sigma}\left(tc_{i\sigma}^{\dagger}c_{j\sigma}+it_{ij}(c_{i,\uparrow}^{\dagger}c_{j,\uparrow}-c_{i,\downarrow}^{\dagger}c_{j,\downarrow})+{\rm H.c.}\right)$
		$\displaystyle+$	$\displaystyle\underset{i,\sigma\sigma^{\prime}}{\sum}\frac{1}{2}U_{\sigma\sigma^{\prime}}n_{i\sigma}\left(n_{i\sigma^{\prime}}-\delta_{\sigma\sigma^{\prime}}\right)+V^{i}_{\rm trap}n_{i,\sigma}-\mu_{i\sigma}n_{i\sigma},$

with $m_{z}=0$.

As shown in Fig. 3, a many-body phase diagram is shown as a function of hopping amplitudes $t$ and $t_{soc}$ for interactions $U_{\uparrow\uparrow}=U_{\downarrow\downarrow}=1.01U_{\uparrow\downarrow}$, based on RBDMFT. To distinguish quantum phases, we introduce superfluid order parameters $\phi_{\sigma}$, pseudospin operators $S^{z}_{i}=\frac{1}{2}({b}^{\dagger}_{i,\uparrow}b_{i,\uparrow}-{b}^{\dagger}_{i,\downarrow}b_{i,\downarrow})$, $S^{x}_{i}=\frac{1}{2}({b}^{\dagger}_{i,\uparrow}b_{i,\downarrow}+{b}^{\dagger}_{i,\downarrow}b_{i,\uparrow})$ and $S^{y}_{i}=\frac{1}{2i}(b^{\dagger}_{i,\uparrow}b_{i,\downarrow}-{b}^{\dagger}_{i,\downarrow}b_{i,\uparrow})$, and winding number $w=\frac{1}{2\pi}\sum_{\mathcal{C}}{\rm arg}(M^{\ast}_{i}M_{j})$ with $M_{i}=S^{x}_{i}+iS^{y}_{i}$ and $\mathcal{C}$ being a closed loop around the center of the trap Zhou et al. (2020); Zhang et al. (2017). We observe five different quantum phases. When $t_{soc}\ll t$, the system demonstrates a ferromagnetic phase ($\rm MI_{I}$), identical with the system without SOAM coupling. With the growth of $t_{soc}$, a spin-vortex phase appears in the Mott-insulating regime with $w\neq 0$ ($\rm MI_{II}$), since the growth of $t_{soc}$ is equivalent to the growth of orbital angular momentum $l$. As shown in Fig. 3(b), SOAM-induced spin rotation appears around the trap center with winding number $w=1$, indicating that the spin rotates slowly with the corresponding response being mainly around the center of the trap. The physical reason is that the SOAM-induced hopping is site-dependent, and pronounced around the trap center, as indicated by Eq. (3). Further increasing $t_{soc}$, the winding number $w$ grows as well, as shown in the inset of Fig. 3(a), and finally we observe the whole system rotating in the regime $t_{soc}\gg t$ ($\rm MI_{III}$). Interestingly, this spin-vortex phase is actually a composite vortex defect, which supports a nonrotating core of antiferromagnetic spin texture, with the nearest-neighbor spins being antiparallel in the trap center, as shown in Fig. 3(b).

To understand the underlying physics in the Mott regime, we treat the hopping as perturbations and derive an effective exchange model at half filling. With defining the projection operators $\mathcal{P}$ and $\mathcal{Q}=1-\mathcal{P}$, we can project the system into the Hilbert space consisting of both singly occupied sites and the states being at least one site with double occupation, and obtain an effective exchange model $H_{eff}=-\mathcal{P}H_{t}\mathcal{Q}(\frac{1}{\mathcal{Q}H_{U}\mathcal{Q}-E})\mathcal{Q}H_{t}\mathcal{P}$ Pinheiro et al. (2013); Pinheiro (2016); Mila and Schmidt (2011). The effective exchange model is finally given by:

	$\displaystyle H_{eff}$	$\displaystyle=$	$\displaystyle\sum_{\langle i,j\rangle}J_{z}S_{i}^{z}S_{j}^{z}+J\left(S_{i}^{x}S_{j}^{x}+S_{i}^{y}S_{j}^{y}\right)$		(20)
		$\displaystyle+$	$\displaystyle D\left(S_{i}\times S_{j}\right)_{z}.$

Here, $J_{z}=-4(\frac{2}{U}-\frac{1}{U_{\uparrow\downarrow}})(t^{2}+t^{2}_{ij})$, $J=-\frac{4(t^{2}-t^{2}_{ij})}{U_{\uparrow\downarrow}}$, and $D=-\frac{8tt_{ij}}{U_{\uparrow\downarrow}}$. The details of derivation are given in the Appendix C.

In the absence of SOAM coupling, this effective model is reduced to the conventional XXZ model, where the system prefers ferromagnetic and anti-ferromagnetic orders Duan et al. (2003); Kuklov and Svistunov (2003); Altman et al. (2003); Isacsson et al. (2005). In the presence of SOAM coupling, a Dzyaloshinskii-Moriya term Fert et al. (2013); Luo et al. (2014) appears in the $z$ direction. This term competes with the normal Heisenberg exchange interactions, resulting spin-vortex defects in the Mott-insulating regime. Interestingly, the Heisenberg exchange term $J$ also depends on the SOAM-induced hopping $t_{soc}$, and dominates in the regime $t_{soc}\gg t$, resulting an antiferromagnetic texture. This texture is consistent with our numerical results, as shown in Fig. 3(b). We remark that the SOAM-induced Dzyaloshinskii-Moriya term $D$ preserves rotational symmetry with spin texture rotating along the azimuthal direction, in contrast to the spin-linear-momentum coupling by breaking lattice-translational symmetry Cole et al. (2012); He et al. (2015); Stanescu et al. (2008); Graß et al. (2011); Mandal et al. (2012).

With the increase of hopping amplitudes, atoms delocalize and the superfluid phase appears. We characterize the superfluid phase with superfluid order parameters $\phi_{\sigma}$. In the superfluid region, we observe two quantum many-body phases, with one being a phase with phase rotating ($\rm SF_{II}$), and the other with conventional phase ($\rm SF_{I}$).

Now we turn to study the whole system, described by Eq. (6), and focus on the stability of spin-vortex texture in the strongly interacting regime. Generally, SOAM coupling preserves rotational symmetry and favors spin-vortex defects Chen et al. (2020b, a), whereas the one-dimensional Raman-induced spin-flip hopping prefers the stripe phase and breaks translational symmetry Wang et al. (2010); Zhai (2015). It is expected that more exotic many-body phases appear, due to the competition between SOAM coupling and Raman-induced spin-flip hopping. Here, we choose two hyperfine states of a ⁸⁷Rb gas as examples, where all the Hubbard parameters are obtained from band-structure simulations ban . To emphasize the influence of SOAM coupling, we consider the orbital angular momenta $l=1$ and $l=2$. Note here that $t\approx t_{soc}$ for orbital angular momentum $l=1$ in the deep lattice, as shown in Fig. 2.

Rich phases are found in Fig. 4, including Mott-insulating phases with ferromagnetic ($\rm MI_{I}$), vortex ($\rm MI_{II}$), and canted-antiferromagnetic ($\rm MI_{IV}$) orders Cooley et al. (2020), and superfluid phases with vortex ($\rm SF_{II}$) and stripe ($\rm SF_{III}$) patterns. In the limit $\Omega\ll t_{soc}$, the many-body phases develop spin-vortex textures ($\rm MI_{II}$ and $\rm SF_{II}$), whereas the system prefers density-wave orders ($\rm MI_{IV}$ and $\rm SF_{III}$) in the limit $\Omega\gg t_{soc}$. This conclusion is consistent with our general discussion above, as a result of the interplay of SOAM coupling and Raman-induced spin-flip hopping. Note here that the region of the spin-vortex phase is enlarged for larger orbital angular momentum, as shown in Fig. 4(b), indicating large opportunity for observing this spin texture for larger orbital angular momentum.

To characterize these different phases, we choose winding number, real-space spin texture, spin-structure factor $S_{\vec{q}}=\left|\vec{S}_{i}e^{i\vec{q}\cdot\vec{r}_{i}}\right|$  Cole et al. (2012), and local phase of superfluid order parameter, as shown in Fig. 5. Here, we choose the orbital angular momentum $l=2$, and different lattice depths $V=14\,E_{R}$ with hopping $t\approx 0.011$ [Fig. 5(a)], and $V=11.5\,E_{R}$ with $t\approx 0.022$ [Fig. 5(b)]. For small $\Omega$, a spin-vortex phase ($\rm MI_{II}$) develops with winding number $w=4$, as shown in Fig. 5(a). Increasing $\Omega$, the spin texture changes to ferromagnetic ($\rm MI_{I}$) and canted-antiferromagnetic ($\rm MI_{IV}$) textures with vanishing winding number, which are characterized both by magnetic spin-structure factor $S_{\vec{q}}$ and real-space spin texture, as shown in the inset of Fig. 5(a). We remark here that the $\rm MI_{IV}$ phase possesses spin-density-wave in $S_{z}$, ferromagnetic order in $S_{x}$, and antiferromagnetic order in $S_{y}$. The local phase of superfluid order parameter is shown in Fig. 5(b). For small $\Omega$, a nonzero winding number of the local phase develops in the vortex superfluid ($\rm SF_{II}$). With $\Omega$ larger, we find the local phase demonstrates a stripe order instead ($\rm SF_{III}$).

To understand the physical phenomena in the Mott regime, we derive an effective exchange model of the system at half filling, described by Eq. (6),

	$\displaystyle H_{eff}$	$\displaystyle=$	$\displaystyle\underset{\langle i_{x},j_{x}\rangle}{\sum}J^{\prime}_{z}S_{i_{x}}^{z}S_{j_{x}}^{z}+J^{\prime}_{x}S_{i_{x}}^{x}S_{j_{x}}^{x}+J^{\prime}_{y}S_{i_{x}}^{y}S_{j_{x}}^{y}$		(21)
		$\displaystyle+$	$\displaystyle\underset{\langle i_{y},j_{y}\rangle}{\sum}J_{z}S_{i_{y}}^{z}S_{j_{y}}^{z}+J\left(S_{i_{y}}^{x}S_{j_{y}}^{x}+S_{i_{y}}^{y}S_{j_{y}}^{y}\right)$
		$\displaystyle+$	$\displaystyle\sum_{\langle i,j\rangle}D\left(S_{i}\times S_{j}\right)_{z},$

where $J_{z}^{\prime}=-4\left(\frac{2\left(t^{2}+t_{ij}^{2}\right)}{U}+\frac{\Omega^{2}}{U_{\uparrow\downarrow}}-\frac{\left(t^{2}+t_{ij}^{2}\right)}{U_{\uparrow\downarrow}}-\frac{2\Omega^{2}}{U}\right)$, $J_{x}^{\prime}=-4\left(\frac{\left(t^{2}-t_{ij}^{2}\right)}{U_{\uparrow\downarrow}}+\frac{\Omega^{2}}{U_{\uparrow\downarrow}}\right)$, and $J_{y}^{\prime}=-4\left(\frac{\left(t^{2}-t_{ij}^{2}\right)}{U_{\uparrow\downarrow}}-\frac{\Omega^{2}}{U_{\uparrow\downarrow}}\right)$. We observe that the Raman-induced spin-flip hopping does not influence Dzyaloshinskii-Moriya interactions, but induce an anisotropy for Heisenberg exchange interactions in the $x$ and $y$ directions. When $\Omega$ is large enough, the Raman-induced hopping can induce $J^{\prime}_{y,z}$ to be positive and $J^{\prime}_{x}$ negative. It indicates that a spin-density-wave and canted-antiferromagnetic order develops for large $\Omega$, consistent with our numerical results, as shown in Fig. 5(a). In the intermediate regime of $\Omega$, the $J^{\prime}_{x}$ term dominates and a ferromagnetic order appears. For small $\Omega$, the effective model reduces to Eq. (20), and the spin vortex pattern dominates.

In a realistic system, one can tune the balance of the two-spin components, which actually acts as an effective magnetic field $m_{z}$. Here, we can control the chemical potential difference of the two components to study the effect of the magnetic field. In Fig. 6, we fix the depth of optical lattice $V=9\,E_{R}$ with hopping $t\approx 0.045$, and study the many-body phase diagram as a function of the effective magnetic field $m_{z}$ and Raman-induced spin-flip hopping $\Omega$. When the effective magnetic field is large and negative, the spin-$\downarrow$ component supports a vortex structure, indicated by the local phase of the superfluid order parameter, as shown in inset of Fig. 6, or vice versa. For large enough Raman coupling, the two components are mixed, and the system supports a stripe pattern in the $x$ direction. We remake here that the phase diagram is similar to the one achieved by the SOAM experiments in continuous space Zhang et al. (2019), where the difference is vortex-antivortex pair phase for larger Raman coupling, instead of stripe order Chen et al. (2020a, b), since we essentially include a Raman lattice in the $x$ direction.

In our scheme, the Raman transition is $\Lambda$-type configuration with $\Omega_{1}(r)=\Omega_{1}{\rm cos}(kx)$ along the $x$ direction and $\Omega_{2}(r)=\Omega_{2}e^{-2r^{2}/\rho^{2}}e^{-i2l\phi}$ along the $z$ direction. In the regime $\left|\Delta\right|\gg\Omega_{1,2}$, the single-photon transition between the ground and excited states is suppressed. We can adiabatically remove the excited state, and the system is effectively regarded as two-ground-state mixtures coupled by two-photon Raman processes. Including the two-photon Raman processes, we can obtain an effective spin-$1/2$ Hamiltonian

$\displaystyle H=\left(\begin{matrix}\frac{p^{2}}{2m}+\frac{\delta}{2}&\Omega^{\prime}\left(r\right)\cos kx\cdot e^{-i2l\phi}\\ \Omega^{\prime}\left(r\right)\cos kx\cdot e^{i2l\phi}&\frac{p^{2}}{2m}-\frac{\delta}{2}\\ \end{matrix}\right),$

(S1)

where $\delta$ is a two-photon detuning, and $\Omega^{\prime}(r)=\Omega_{1}\Omega_{2}e^{-2r^{2}/\rho^{2}}/\Delta$ denotes the Raman Rabi frequency. After introducing the unity transformation to the single-particle wave function

$\displaystyle U$

$\displaystyle=$

$\displaystyle\left(\begin{matrix}e^{-il\phi}&0\\ 0&e^{il\phi}\\ \end{matrix}\right),$

(S2)

and Pauli matrix $\sigma$, we finally obtain

$\displaystyle H$

$\displaystyle=$

$\displaystyle-\frac{\hbar^{2}}{2mr}\frac{\partial}{\partial_{r}}\left(r\frac{\partial}{\partial_{r}}\right)+\frac{\delta}{2}\sigma_{z}+\frac{\left(L_{z}-l\hbar\sigma_{z}\right)}{2mr^{2}}^{2}+\Omega^{\prime}\left(r\right)\cos kx\cdot\sigma_{x},$

where $L_{z}=-i\hbar\partial_{z}$ is the orbital angular momentum along the $z$ axis. Normally, the plan-wave laser is a LG beam with the intensity being suppressed near the trap center, which can influence experimental observations. Here, we instead consider a Gaussian-type Raman beam with orbital angular momentum, where such a Gaussian beam can be obtained by a quarter-wave plate Chen et al. (2020c); Rafayelyan and Brasselet (2017). The waist of the plane-wave laser is set to $\rho=20$.

The angular momentum $t_{ij}$ is given by

	$\displaystyle\mathcal{H}_{L_{z}}$	$\displaystyle=$	$\displaystyle\int dx\Psi^{\dagger}(x)\frac{l\hbar}{mr^{2}}L_{z}\Psi(x)$		(S4)
		$\displaystyle=$	$\displaystyle\frac{l\hbar}{m}\int dx\Psi^{\dagger}(x)\frac{1}{r^{2}}L_{z}\Psi(x).$

For a sufficient deep lattice, the field operator $\Psi(x)$ can be expanded in the lowest-band Wannier basis $\omega\left(\boldsymbol{x}-\boldsymbol{R}_{i}\right)$. Eq. (S4) can be rewritten as

$\displaystyle\mathcal{H}_{L_{z}}$

$\displaystyle=$

$\displaystyle\frac{l\hbar}{m}\sum_{\langle i,j\rangle}c^{\dagger}_{i}c_{j}\int{dx^{3}}\omega^{\ast}\left(\boldsymbol{x}-\boldsymbol{R}_{i}\right)\frac{1}{r^{2}}L_{z}\omega\left(\boldsymbol{x}-\boldsymbol{R}_{j}\right).$

(S5)

For nearest neighbors $i$ and $j$, we have the relation that

$\displaystyle\frac{1}{r_{1}^{2}}K_{i,j}<\int{dx^{3}}\omega^{\ast}\left(\boldsymbol{x}-\boldsymbol{R}_{i}\right)\frac{1}{r^{2}}L_{z}\omega\left(\boldsymbol{x}-\boldsymbol{R}_{j}\right)<\frac{1}{r_{2}^{2}}K_{i,j},$

where $r_{1}={\rm max}(r_{i},r_{j})$, and $r_{2}={\rm min}(r_{i},r_{j})$, with $r_{i}$ ($r_{j}$) being the distance between site $i$ ($j$) and the trap center. For simplicity, we take the distance between the midpoint of sites $i$, $j$, and the trap center as the approximation of $r$, and denote it by $r^{\prime}$. Eq. (S4) reads

$\displaystyle\mathcal{H}_{L_{z}}$

$\displaystyle\approx$

$\displaystyle\frac{l\hbar}{mr^{\prime 2}}\sum_{\langle i,j\rangle}c^{\dagger}_{i}c_{j}K_{i,j},$

(S7)

with $K_{i,j}=\int{dx^{3}}\omega^{\ast}\left(\boldsymbol{x}-\boldsymbol{R}_{i}\right)L_{z}\omega\left(\boldsymbol{x}-\boldsymbol{R}_{j}\right)=-i\hbar\int{dx^{3}}\omega^{\ast}\left(\boldsymbol{x}-\boldsymbol{R}_{i}\right)(x\partial_{y}-y\partial_{x})\omega\left(\boldsymbol{x}-\boldsymbol{R}_{j}\right)$. Generally, the Wannier function can be factorized into the $x$- and $y$-dependent parts for the deep lattice, and we finally obtain

	$\displaystyle K_{i,j}$	$\displaystyle=$	$\displaystyle-i\hbar\left(\int dx\omega^{\ast}(x-x_{i})x\omega(x-x_{j})\int dy\omega^{\ast}(y-y_{i})\partial_{y}\omega(y-y_{j})\right.$		(S8)
		$\displaystyle-$	$\displaystyle\left.\int dx\omega^{\ast}(x-x_{i})\partial_{x}\omega(x-x_{j})\int dy\omega^{\ast}(y-y_{i})y\omega(y-y_{j})\right).$

For the discrete lattice system, $K_{ij}$ can be written in a product form

$\displaystyle K_{i,j}$

$\displaystyle=$

$\displaystyle-i\hbar\left(\frac{x_{i}y_{j}-x_{j}y_{i}}{d}\alpha\right),$

(S9)

with $\alpha=\int dx\omega^{\ast}\left(x-d\right)\partial_{x}\omega\left(x\right)$ and $d$ being the lattice constant, as shown in Fig. S1. We finally obtain the orbital angular momentum in the Wannier basis

	$\displaystyle\mathcal{H}_{L_{z}}$	$\displaystyle=$	$\displaystyle\sum_{\langle i,j\rangle}-i\frac{l\hbar^{2}}{mr^{\prime 2}}\frac{x_{i}y_{j}-x_{j}y_{i}}{d}\left(c^{\dagger}_{i}c_{j}-c_{i}c^{\dagger}_{j}\right)\int dx\omega^{\ast}\left(x-d\right)\partial_{x}\omega\left(x\right)$		(S10)
		$\displaystyle=$	$\displaystyle\sum_{\langle i,j\rangle}i\frac{x_{i}y_{j}-x_{j}y_{i}}{r^{\prime 2}}\left(-\frac{l\hbar^{2}}{dm}\int dx\omega^{\ast}\left(x-d\right)\partial_{x}\omega\left(x\right)\right)\left(c^{\dagger}_{i}c_{j}-c_{i}c^{\dagger}_{j}\right).$

The system can be described by an effective exchange model in the deep Mott-insulating regime. To derive the effective model, we first divide the Hilbert space according into site occupations for filling $n=1$. We define the operators $\mathcal{P}$ and $\mathcal{Q}$ , which denote the projection into the Mott state subspace $\mathcal{H}_{P}$ and the perpendicular subspace $\mathcal{H}_{Q}$. For a Hamiltonian $H$, the Schrödinger equation reads

$\displaystyle H\left|\psi\right>=E\left|\psi\right>.$

(S11)

With the unity operator $1=\mathcal{P}+\mathcal{Q}$, we obtain

$\displaystyle H(\mathcal{P}+\mathcal{Q})\left|\psi\right>=E(\mathcal{P}+\mathcal{Q})\left|\psi\right>.$

(S12)

Multiplying by $\mathcal{P}$ and $\mathcal{Q}$ the left side of Eq. (S12) results in

	$\displaystyle(\mathcal{P}H\mathcal{P}+\mathcal{P}H\mathcal{Q})|\psi\rangle$	$\displaystyle=$	$\displaystyle E\mathcal{P}|\psi\rangle,$		(S13)
	$\displaystyle(\mathcal{Q}H\mathcal{P}+\mathcal{Q}H\mathcal{Q})|\psi\rangle$	$\displaystyle=$	$\displaystyle E\mathcal{Q}|\psi\rangle.$		(S14)

Eq. (S14) can be rewritten with the projection operator relation $\mathcal{Q}^{2}=\mathcal{Q}$,

$\displaystyle\mathcal{Q}\left|\psi\right>=\frac{1}{E-\mathcal{Q}H\mathcal{Q}}\mathcal{Q}H\mathcal{P}|\psi\rangle.$

(S15)

Inserting (S15) into (S13), we obtain an equation for $\mathcal{P}|\psi\rangle$,

$\displaystyle(\mathcal{P}H\mathcal{Q}\frac{1}{E-\mathcal{Q}H\mathcal{Q}}\mathcal{Q}H\mathcal{P})\mathcal{P}|\psi\rangle=E\mathcal{P}|\psi\rangle.$

(S16)

For Hamiltonian $H$, we can divide it into two parts $H=H_{t}+H_{U}$, where $H_{t}$ and $H_{U}$ are the hopping and interaction parts of $H$, respectively. Eq. (S16) can be rewritten as

$\displaystyle(\mathcal{P}H_{t}\mathcal{Q}\frac{1}{E-\mathcal{Q}H_{U}\mathcal{Q}-\mathcal{Q}H_{t}\mathcal{Q}}\mathcal{Q}H_{t}\mathcal{P})\mathcal{P}|\psi\rangle=E\mathcal{P}|\psi\rangle,$

(S17)

with $\mathcal{P}H_{U}\mathcal{P}$, $\mathcal{P}H_{U}\mathcal{Q}$ and $\mathcal{P}H_{t}\mathcal{P}$ being zero. Because $E\sim\frac{t^{2}}{U}\sim 0$ in the Mott phase, we take $\frac{1}{E-\mathcal{Q}H_{U}\mathcal{Q}-\mathcal{Q}H_{t}\mathcal{Q}}\approx\frac{1}{-\mathcal{Q}H_{U}\mathcal{Q}-\mathcal{Q}H_{t}\mathcal{Q}}$. Using the expansion $\frac{1}{A-B}=\frac{1}{A}\Sigma^{\infty}_{n=0}(B\frac{1}{A})^{n}$, with $A=-\mathcal{Q}H_{U}\mathcal{Q}$ and $B=\mathcal{Q}H_{t}\mathcal{Q}$, we finally obtain

$\displaystyle\mathcal{H}_{eff}$

$\displaystyle=$

$\displaystyle\mathcal{P}H_{t}\mathcal{Q}\frac{1}{-\mathcal{Q}H_{U}\mathcal{Q}}\cdot\Sigma^{\infty}_{n=0}(\mathcal{Q}H_{t}\mathcal{Q}\frac{1}{\mathcal{Q}H_{U}\mathcal{Q}})^{n}\mathcal{Q}H_{t}\mathcal{P}.$

(S18)

Normally, we only need to include nearest-neighbor terms in the effective Hamiltonian with $n=0$, i.e.

$\displaystyle\mathcal{H}_{eff}$

$\displaystyle=$

$\displaystyle\mathcal{P}H_{t}\mathcal{Q}\frac{1}{-\mathcal{Q}H_{U}\mathcal{Q}}\mathcal{Q}H_{t}\mathcal{P}.$

(S19)

In the tight-binding regime, the subspace $\mathcal{H}_{P}$ of the states with half filling for a two-site problem is

$\displaystyle\mathcal{H}_{P}:\left\{\left|\uparrow,\uparrow\right>,\left|\uparrow,\downarrow\right>,\left|\downarrow,\downarrow\right>,\left|\downarrow,\downarrow\right>\right\},$

(S20)

where $\left|\sigma,\sigma^{\prime}\right>$ denotes the spin state $\sigma$ in the left site and $\sigma^{\prime}$ in the right one. The subspace $\mathcal{H}_{P}$ for doubly occupied sites is

$\displaystyle\mathcal{H}_{Q}:\left\{\left|\uparrow\uparrow\right>,\left|\uparrow\downarrow\right>,\left|\downarrow\downarrow\right>\right\}.$

(S21)

In a matrix form, $\left(\mathcal{Q}H_{U}\mathcal{Q}\right)^{-1}$ is given by

$\displaystyle\left(\mathcal{Q}H_{U}\mathcal{Q}\right)^{-1}$

$\displaystyle=$

$\displaystyle\left(\begin{matrix}\frac{1}{U}&0&0\\ 0&\frac{1}{U}&0\\ 0&0&\frac{1}{U_{\uparrow\downarrow}}\\ \end{matrix}\right).$

(S22)

According to Eq. (S19), we obtain the final effective Hamiltonian

	$\displaystyle\mathcal{H}_{eff}$	$\displaystyle=$	$\displaystyle-\left(\frac{4\left(t^{2}+t_{ij}^{2}\right)}{U}+\frac{2\Omega^{2}}{U_{\uparrow\downarrow}}\right)\left(n_{i,\uparrow}n_{j,\uparrow}+n_{i,\downarrow}n_{j,\downarrow}\right)-\left(\frac{2\left(t^{2}+t_{ij}^{2}\right)}{U_{\uparrow\downarrow}}+\frac{4\Omega^{2}}{U}\right)\left(n_{i,\uparrow}n_{j,\downarrow}+n_{i,\downarrow}n_{j,\uparrow}\right)$
		$\displaystyle-$	$\displaystyle\frac{2\left(-t+it_{ij}\right)^{2}}{U_{\uparrow\downarrow}}c_{i,\downarrow}^{\dagger}c_{i,\uparrow}c_{j,\uparrow}^{\dagger}c_{j,\downarrow}-\frac{2\left(-t-it_{ij}\right)^{2}}{U_{\uparrow\downarrow}}c_{i,\uparrow}^{\dagger}c_{i,\downarrow}c_{j,\downarrow}^{\dagger}c_{j,\uparrow}$
		$\displaystyle-$	$\displaystyle\frac{2\Omega^{2}}{U_{\uparrow\downarrow}}\left(c_{i,\uparrow}^{\dagger}c_{i,\downarrow}c_{j,\uparrow}^{\dagger}c_{j,\downarrow}+c_{i,\downarrow}^{\dagger}c_{i,\uparrow}c_{j,\downarrow}^{\dagger}c_{j,\uparrow}\right)-(-1)^{i_{x}}\left(\frac{2\Omega\left(-t-it_{ij}\right)}{U}+\frac{2\Omega\left(-t-it_{ij}\right)}{U_{\uparrow\downarrow}}\right)\left(c_{j,\downarrow}^{\dagger}c_{j,\uparrow}+c_{i,\uparrow}^{\dagger}c_{i,\downarrow}\right)$
		$\displaystyle-$	$\displaystyle(-1)^{i_{x}}\left(\frac{2\Omega\left(-t+it_{ij}\right)}{U_{\uparrow\downarrow}}+\frac{2\Omega\left(-t+it_{ij}\right)}{U}\right)\left(c_{j,\uparrow}^{\dagger}c_{j,\downarrow}+c_{i,\downarrow}^{\dagger}c_{i,\uparrow}\right).$

Introducing the pseudospin operator as follows

	$\displaystyle S^{x}_{i}$	$\displaystyle=$	$\displaystyle\frac{1}{2}\left(c^{\dagger}_{i,\uparrow}c_{i,\downarrow}+c^{\dagger}_{i,\downarrow}c_{i,\uparrow}\right)$		(S24)
	$\displaystyle S^{y}_{i}$	$\displaystyle=$	$\displaystyle\frac{1}{2i}\left(c^{\dagger}_{i,\uparrow}c_{i,\downarrow}-c^{\dagger}_{i,\downarrow}c_{i,\uparrow}\right)$		(S25)
	$\displaystyle S^{z}_{i}$	$\displaystyle=$	$\displaystyle\frac{1}{2}\left(n_{i,\uparrow}-n_{i,\downarrow}\right),$		(S26)

Eq. (LABEL:effec1) can be rewritten as

	$\displaystyle\mathcal{H}_{eff}$	$\displaystyle=$	$\displaystyle-4\left(\frac{2\left(t^{2}+t_{ij}^{2}\right)}{U}+\frac{\Omega^{2}}{U_{\uparrow\downarrow}}-\frac{\left(t^{2}+t_{ij}^{2}\right)}{U_{\uparrow\downarrow}}-\frac{2\Omega^{2}}{U}\right)S_{i}^{z}S_{j}^{z}-4\left(\frac{\left(t^{2}-t_{ij}^{2}\right)}{U_{\uparrow\downarrow}}+\frac{\Omega^{2}}{U_{\uparrow\downarrow}}\right)S_{i}^{x}S_{j}^{x}$		(S27)
		$\displaystyle-$	$\displaystyle 4\left(\frac{\left(t^{2}-t_{ij}^{2}\right)}{U_{\uparrow\downarrow}}-\frac{\Omega^{2}}{U_{\uparrow\downarrow}}\right)S_{i}^{y}S_{j}^{y}-\frac{8tt_{ij}}{U_{\uparrow\downarrow}}\left(S_{i}\times S_{j}\right)_{z}$
		$\displaystyle+$	$\displaystyle(-1)^{i_{x}}\Omega t\left(\frac{4}{U}+\frac{4}{U_{\uparrow\downarrow}}\right)\left(S_{i}^{x}+S_{j}^{x}\right)-(-1)^{i_{x}}\Omega t_{ij}\left(\frac{4}{U}+\frac{4}{U_{\uparrow\downarrow}}\right)\left(S_{i}^{y}-S_{j}^{y}\right).$

The result can be easily extended to the case with vanishing $\Omega=0$, and it reads

	$\displaystyle\mathcal{H}_{eff}$	$\displaystyle=$	$\displaystyle-4\left(\frac{2\left(t^{2}+t_{ij}^{2}\right)}{U}-\frac{\left(t^{2}+t_{ij}^{2}\right)}{U_{\uparrow\downarrow}}\right)S_{i}^{z}S_{j}^{z}-\frac{4\left(t^{2}-t_{ij}^{2}\right)}{U_{\uparrow\downarrow}}S_{i}^{x}S_{j}^{x}$
		$\displaystyle-$	$\displaystyle\frac{4\left(t^{2}-t_{ij}^{2}\right)}{U_{\uparrow\downarrow}}S_{i}^{y}S_{j}^{y}-\frac{8tt_{ij}}{U_{\uparrow\downarrow}}\left(S_{i}\times S_{j}\right)_{z}.$