Magnetic phase transitions of insulating spin-orbit coupled Bose atoms in one-dimensional optical lattices

TL;DR

This paper investigates magnetic phase transitions in insulating spin-orbit coupled Bose atoms in one-dimensional optical lattices, revealing how interactions and external fields influence magnetic order and quantum phases.

Contribution

It introduces a detailed analysis of magnetic phases using matrix product states, highlighting the effects of spin-orbit coupling, interactions, and transverse fields on quantum magnetism.

Findings

Gapped ferromagnetic phase with long-range order when no transverse field.

Gapless Luttinger liquid phase with algebraic correlations under certain interactions.

Emergence of long-range ordered phases with transverse field and Dzyaloshinskii-Moriya interaction.

Abstract

We consider the insulating spin-orbit coupled Bose atoms confined within one-dimensional optical lattices and explore their ground-state magnetic phase transitions. Under strong interactions, the charge degrees of atoms are frozen and the system can be described by an anisotropic XXZ Heisenberg chain with Dzyaloshinskii-Moriya interaction and transverse field. We apply the matrix product state method to obtain low-energy states and analyze the lowest energy gaps and the ground-state magnetization and correlations. We find when the transverse field is absent, the ground state is a gapped ferromagnetic phase with long-range correlation in the z direction if the interspin s-wave interacting strength is stronger than that of the intraspin one, otherwise it is a gapless Luttinger liquid (LL) phase with algebraic decaying correlation. When the transverse field is turned on, the gapless LL phase is broken, there emerges a long-range correlated phase with ferromagnetic, antiferromagnetic or spiral order, which depends on the DM interaction strength. We believe our study provides a complete understanding of the interplay between SOC and quantum magnetism of spinor atoms in optical lattices.

Equations

$${\hat{H}}_t=$$

$${\hat{H}}_U=\frac{U}{2}\sum _{j,\sigma }{\hat{n}}_{j\sigma }\left({\hat{n}}_{j\sigma }-1\right)+U_{\uparrow \downarrow }\sum _j{\hat{n}}_{j\uparrow }{\hat{n}}_{j\downarrow },$$

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

$$|\tilde{\psi }⟩$$

$$|\psi ⟩$$

$$\sum _{a_l}^{{\chi }_l}{U}_{a_{l-1}{\sigma }_l,a_l}{S}_{a_l}{V}_{a_l,{\sigma }_{l+1}\cdots {\sigma }_L}^{†}={\Psi }_{a_{l-1}{\sigma }_l,{\sigma }_{l+1}\cdots {\sigma }_L},$$

$$A_{a_{l-1},a_l}^{{\sigma }_l[l]}=U_{a_{l-1}{\sigma }_l,a_l},$$

$${\Psi }_{a_l{\sigma }_{l+1},{\sigma }_{l+2}\cdots {\sigma }_L}={\left(S{V}^{†}\right)}_{a_l,{\sigma }_{l+1}\cdots {\sigma }_L},$$

$$\sum _{{\sigma }_l}{A}^{{\sigma }_l[l]†}{A}^{{\sigma }_l[l]}=I,\left(l<L\right)$$

$$\sum _{{\sigma }_L}{A}^{{\sigma }_L[l]†}{A}^{{\sigma }_L[L]}=⟨\psi |\psi ⟩.$$

$$\sum _{a_{l-1}}^{{\chi }_{l-1}}{U}_{{\sigma }_1\cdots {\sigma }_{l-1},a_{l-1}}{S}_{a_{l-1},a_{l-1}}{V}_{a_{l-1},{\sigma }_l{a}_l}^{†}={\Psi }_{{\sigma }_1\cdots {\sigma }_{l-1},{\sigma }_l{a}_l},$$

$$B_{a_{l-1},a_l}^{{\sigma }_l[l]}=V_{a_{l-1},{\sigma }_l{a}_l}^{†},$$

$${\Psi }_{{\sigma }_1\cdots {\sigma }_{l-2},{\sigma }_{l-1}{a}_{l-1}}={\left(US\right)}_{{\sigma }_1\cdots {\sigma }_{l-1},a_{l-1}},$$

$$|\psi ⟩$$

$$\sum _{{\sigma }_l}{B}^{{\sigma }_l[l]}{B}^{{\sigma }_l[l]†}=I,\left(l>1\right)$$

$$\sum _{{\sigma }_1}{B}^{{\sigma }_1[1]}{B}^{{\sigma }_1[1]†}=⟨\psi |\psi ⟩,$$

$${\tilde{M}}_{a_{l-1}{\sigma }_l,a_l}=M_{a_{l-1},a_l}^{{\sigma }_l[l]},$$

$$\sum _{a_l^{\prime }}{U}_{a_{l-1}{\sigma }_l,a_l^{\prime }}^{[\mathrm{L}]}{S}_{a_l^{\prime }}^{[\mathrm{L}]}{V}_{a_l^{\prime },a_l}^{[\mathrm{L}]†}={\tilde{M}}_{a_{l-1}{\sigma }_l,a_l},$$

$$A_{a_{l-1},a_l^{\prime }}^{{\sigma }_l[l]}=U_{a_{l-1}{\sigma }_l,a_l^{\prime }}^{[\mathrm{L}]},$$

$$M_{a_l^{\prime },a_{l+1}}^{{\sigma }_{l+1}[l+1]}=\sum _{a_l}{\left(S^{[\mathrm{L}]}{V}^{[\mathrm{L}]†}\right)}_{a_l^{\prime },a_l}{M}_{a_l,a_{l+1}}^{{\sigma }_{l+1}[l+1]},$$

$${\tilde{M}}_{a_{l-1},{\sigma }_l{a}_l}=M_{a_{l-1},a_l}^{{\sigma }_l[l]},$$

$$\sum _{a_{l-1}^{\prime }}{U}_{a_{l-1},a_{l-1}^{\prime }}^{[\mathrm{R}]}{S}_{a_{l-1}^{\prime }}^{[\mathrm{R}]}{V}_{a_{l-1}^{\prime },{\sigma }_l{a}_l}^{[\mathrm{R}]†}={\tilde{M}}_{a_{l-1},{\sigma }_l{a}_l},$$

$$B_{a_{l-1}^{\prime },a_l}^{{\sigma }_l[l]}=V_{a_{l-1}^{\prime },{\sigma }_l{a}_l}^{[\mathrm{R}]†},$$

$$M_{a_{l-2},a_{l-1}^{\prime }}^{{\sigma }_{l-1}[l-1]}=M_{a_{l-2},a_{l-1}}^{{\sigma }_{l-1}[l-1]}{\left(U^{[\mathrm{R}]}{S}^{[\mathrm{R}]}\right)}_{a_{l-1},a_{l-1}^{\prime }},$$

$$\hat{O}$$

$${\hat{W}}^{[1]}$$

$${\hat{W}}^{[1<l<L]}$$

$$\epsilon \left[|\psi ⟩\right]=⟨\psi |\hat{H}|\psi ⟩-E⟨\psi |\psi ⟩,$$

$$⟨\psi |\hat{H}|\psi ⟩=\sum _{{\sigma }_k,{\sigma }_k^{\prime }}\sum _{a_{k-1},a_k}\sum _{a_{k-1}^{\prime },a_k^{\prime }}\sum _{b_{k-1},b_k}\left(L_{a_{k-1},a_{k-1}^{\prime }}^{b_{k-1}[k-1]}{M}_{a_{k-1},a_k}^{{\sigma }_k[k]\ast }{W}_{b_{k-1},b_k}^{{\sigma }_k,{\sigma }_k^{\prime }[k]}{M}_{a_{k-1}^{\prime },a_k^{\prime }}^{{\sigma }_k^{\prime }[k]}{R}_{a_k,a_k^{\prime }}^{b_k[k+1]}\right),$$

$$L_{a_l,a_l^{\prime }}^{b_l[l]}=\sum _{{\sigma }_l,{\sigma }_l^{\prime }}\sum _{a_{l-1},a_{l-1}^{\prime }}\sum _{b_{l-1}}{L}_{a_{l-1},a_{l-1}^{\prime }}^{b_{l-1}[l-1]}{A}_{a_{l-1},a_l}^{{\sigma }_l[l]\ast }{W}_{b_{l-1},b_l}^{{\sigma }_l,{\sigma }_l^{\prime }[l]}{A}_{a_{l-1}^{\prime },a_l^{\prime }}^{{\sigma }_l^{\prime }[l]},$$

$$R_{a_{l-1},a_{l-1}^{\prime }}^{b_{l-1}[l]}=\sum _{{\sigma }_l,{\sigma }_l^{\prime }}\sum _{a_l,a_l^{\prime }}\sum _{b_l}{B}_{a_{l-1},a_l}^{{\sigma }_l[l]\ast }{W}_{b_{l-1},b_l}^{{\sigma }_l,{\sigma }_l^{\prime }[l]}{B}_{a_{l-1}^{\prime },a_l^{\prime }}^{{\sigma }_l^{\prime }[l]}{R}_{a_l,a_l^{\prime }}^{b_l[l+1]},$$

$$\sum _{{\sigma }_k^{\prime }}\sum _{a_{k-1}^{\prime },a_k^{\prime }}\sum _{b_{k-1},b_k}{L}_{a_{k-1},a_{k-1}^{\prime }}^{b_{k-1}[k-1]}{W}_{b_{k-1},b_k}^{{\sigma }_k,{\sigma }_k^{\prime }[k]}{R}_{a_k,a_k^{\prime }}^{b_k[k+1]}{M}_{a_{k-1}^{\prime },a_k^{\prime }}^{{\sigma }_k^{\prime }[k]}-E{M}_{a_{k-1},a_k}^{{\sigma }_k[k]}=0.$$

$${\hat{H}}_{a_{k-1}{\sigma }_k{a}_k,a_{k-1}^{\prime }{\sigma }_k^{\prime }{a}_k^{\prime }}^{[k]}=\sum _{b_{k-1},b_k}{L}_{a_{k-1},a_{k-1}^{\prime }}^{b_{k-1}[k-1]}{W}_{b_{k-1},b_k}^{{\sigma }_k,{\sigma }_k^{\prime }[k]}{R}_{a_k,a_k^{\prime }}^{b_k[k+1]}.$$

$$⟨{\psi }_n|{\varphi }_m⟩=0,$$

$$\left({\hat{P}}^{[k]†}{\hat{H}}^{[k]}{\hat{P}}^{[k]}\right)|v^{[k]}⟩-E|v^{[k]}⟩=0,$$

$$F_{a_{k-1},a_k}^{{\sigma }_k[k](m)}=\sum _{a_{k-1}^{\prime },a_k^{\prime }}{\mathcal{L}}_{a_{k-1},a_{k-1}^{\prime }}^{[k-1](m)}{A}_{a_{k-1}^{\prime },a_k^{\prime }}^{{\sigma }_k[k](m)}{\mathcal{R}}_{a_k,a_k^{\prime }}^{[k+1](m)},$$

$${\mathcal{L}}_{a_l,a_l^{\prime }}^{[l](m)}=\sum _{a_{l-1},{\sigma }_l,a_{l-1}^{\prime }}{M}_{a_{l-1},a_l}^{{\sigma }_l[l]\ast }{\mathcal{L}}_{a_{l-1},a_{l-1}^{\prime }}^{[l-1](m)}{A}_{a_{l-1}^{\prime },a_l^{\prime }}^{{\sigma }_l[l](m)},$$

Full text

I introduction

Ultracold spinor atoms in optical lattice provide an excellent platform for simulating magnetic phase transitions in quantum Heisenberg models. In the Mott regime, the fluctuation of charge degree of spinor atoms is suppressed and the low-energy physics can be captured by an effective spin superexchange model Kuklov and Svistunov (2003); Duan et al. (2003); Lee (2004), in which the spin-spin coupling can be tuned via spin-dependent s-wave scattering and lattice depth. In recent years, there is a lot of interest in creating synthetic gauge fields and spin-orbital coupling (SOC) in ultracold system Galitski and Spielman (2013); Zhai (2015); Zhang and Jo (2018); Zhang and Liu (2018). In optical lattice systems, SOC is engineered via dressing optical lattice with periodic Raman potentials Liu et al. (2013); Wu et al. (2016); Song et al. (2018), in which the SOC induces nearest-neighbor spin-flip; and ladder-like system subject to gauge fields (the ladder legs labeled by either the internal spin states or the real-space lattices, which can be mapped to effective spin Hügel and Paredes (2014)) by Raman-assisted tunneling Celi et al. (2014); Stuhl et al. (2015); Mancini et al. (2015); Atala et al. (2014) or optical clock transition Wall et al. (2016); Livi et al. (2016); Kolkowitz et al. (2017), in which the SOC induces on-site spin-flip.

In addition to the realistic atom-atom interaction, the SOC plays a key role in magnetic phase transitions. It may adjust the anisotropic couplings and then lead to the so-called Dzyaloshinskii-Moriya (DM) Dzyaloshinsky (1958); Moriya (1960) exchange interaction, which induces exotic magnetic phases. In two-dimensional (2D) systems, besides normal ferromagnetic and antiferromagnetic phases, the coexistence of DM interaction and anisotropy give spiral phases, vortex crystal structure and novel Skymion Cai et al. (2012); Cole et al. (2012); Radić et al. (2012); Gong et al. (2015); Wang et al. (2016). Classical phase diagrams have been obtained by Monte Carlo simulations Cole et al. (2012); Gong et al. (2015), steepest descent minimization method Radić et al. (2012) and variational mean-field approach Wang et al. (2016). In one-dimensional systems, the effective model reduces into an anisotropic Heisenberg chain with one-component DM interaction. The novel structures in 2D systems is simplified to gapless spiral phase on the plane perpendicular to the DM vector, characterized by algebraic decaying correlation function, signaling a gapless Luttinger liquid (LL) phase Zhao et al. (2014); Piraud et al. (2014); Xu et al. (2014). The various phases and quantum phase diagrams in 1D systems are widely studied Peotta et al. (2014); Zhao et al. (2014, 2014); Piraud et al. (2014); Xu et al. (2014) using density matrix renormalization group method.

Most previous works concentrate on nearest-neighbor spin-flip induced by SOC, where the effective models do not contain any external fields, in particular, spin-flipping fields. An antiferromagnetic XXZ model with DM interaction and transverse field, which can be mapped into an isotropic antiferromagnetic XXX spin chain in spiral field, can be realized in a 1D fermionic ladderlike optical lattice penetrated by synthetic magnetic field Sun et al. (2013). The spiral field provides a new flexible freedom to tune the phase diagram. However, there are only two phases in this model: a gapless LL phase and a ferromagnetic phase, due to the absence of anisotropy. One can alternately load spinor (two-component) bosons into such lattice, where the spin-dependent s-wave scattering can cause anisotropic spin-spin coupling. Furthermore, since both fermionic Mancini et al. (2015); Wall et al. (2016); Livi et al. (2016); Kolkowitz et al. (2017) and bosonic ladder subject to gauge field Stuhl et al. (2015) have been realized in experiments, there is no reason to leave the bosonic case unexplored. It is interesting to study the interplay of anisotropy, DM interaction and spin-flipping external field.

In this article, we explore the effect of SOC on the Mott insulator of two-level Bose atoms in 1D optical lattices. Utilizing the SOC-dressed Hubbard model realized in a recent experiment Stuhl et al. (2015), we derive an effective spin model in the strong repulsive interaction limit: an XXZ model with DM interaction subject to external transverse field. The strength of the anisotropy, DM interaction and transverse field of the spin model can be tuned by the parameters of the optical lattice and the SOC strength. By employing a variational matrix product state (MPS) search Verstraete et al. (2008); Schollwöck (2011); Wall and Carr (2012), we obtain the low-energy states. Then, by calculating the lowest energy gap, and the spin-spin correlation function for the ground state (GS), we identify four typical phases characterized by long-range correlation: (I) $z$-FM phase with ferromagnetic correlation along the $\hat{z}$ direction; (II) $x$-PARA phase with only ferromagnetic correlation along the external field direction; (III) $y$-AFM phase with antiferromagnetic correlation along the $\hat{y}$ direction; and (IV) $xy$-SP phase with spiral correlation on the $\hat{x}\hat{y}$ plane. In particular, we find although the transverse field does not open the gap of the spiral phase on the $\hat{x}\hat{y}$ plane, it induces long-range correlation, distinct from the gapless LL phase found in the SOC induced nearest-neighbor spin-flip hopping model and the fermionic ladder penetrated by synthetic magnetic field. Finally, via analyzing the order parameters, we give rich GS phase diagrams in the full parameter range.

The article structure is as follows. In this section, we introduce the related background and our motivation. In Sec. II, we derive the effective Hamiltonian for our physical system. In Sec. III, we review the MPS method for obtaining the low-energy states. In Sec. IV, we calculate the lowest energy gaps, two-site spin correlation function and the order parameters for determining the different phases and phase boundaries. In the last section, we summarize and discuss our results.

II Strongly interacting spin-orbit coupled Bose atoms in 1D optical lattices

Recent experiments have realized spin-orbit coupled bosons in 1D optical lattice by coupling the three internal states (pseudospins) with Raman-assisted transition Stuhl et al. (2015). It is a virtual three-leg ladder pierced by magnetic flux $\phi$, with the extra dimension being the internal states of the bosons. It can be reduced to a two-leg ladder if the second order Zeeman shift is large enough that the upper internal level can be removed Brion et al. (2007). Moreover, such a bosonic ladder may also be realized by coupling metastable states of bosonic atoms by optical clock transition, which has been utilized to realize fermionic ladder in gauge field Wall et al. (2016); Livi et al. (2016); Kolkowitz et al. (2017). Choosing the Landau gauge, for which the phase is accumulated by intra-leg hopping and resulting a net flux of $\phi$ through each plaquette, the single-particle Hamiltonian of this system can be written as

[TABLE]

where $\hat{a}^{\dagger}_{j,\sigma}(\sigma=\uparrow,\downarrow)$ creates a particle of internal state (pseudospin) $\sigma$ at site $j$. The first term is the spin-conserved nearest-neighbor hopping with strength $te^{\pm i\frac{\phi}{2}}$, where $t$ can be tuned by the depth of the optical lattice, and the magnetic flux $\phi$ is related to the SOC momentum $k_{\mathrm{SOC}}=\phi k_{L}/\pi$ ($k_{L}$ is the lattice momentum), which is the momentum transfer of the Raman lasers and can be tuned by choosing different wavelength and/or changing the relative angle of them. The second term describes the on-site spin-flip with strength $\frac{\Omega}{2}$, with $\Omega$ being the Rabi frequency of the Raman lasers. Multiple particles in this optical lattice are described by the Hamiltonian $\hat{H}=\hat{H}_{t}+\hat{H}_{U}$, with the interaction term

[TABLE]

where $\hat{n}_{j\sigma}=\hat{a}^{\dagger}_{j\sigma}\hat{a}_{j\sigma}$ is the particle number operator. The on-site interspin and intraspin interaction strengthes are denoted as $U_{\uparrow\downarrow}$ and $U_{\sigma\sigma}$ respectively, which can be tuned by Feshbach resonance. In the following, we set repulsive interaction as $U_{\uparrow\uparrow}=U_{\downarrow\downarrow}=U>0$ and $U_{\uparrow\downarrow}=\lambda U>0$.

We are interested in magnetic properties in the deep Mott insulator regime at half filling with strong interaction $U,U_{\uparrow\downarrow}\gg t,\Omega$. We treat the tunneling Hamiltonian $\hat{H}_{t}$ as perturbations to the on-site interaction $\hat{H}_{U}$. The GS of $\hat{H}_{U}$ is a Mott insulator with exactly one particle per site. It is many-fold degenerate since the spin on every site is arbitrary. $\hat{H}_{t}$ couples the manifold ground states of $\hat{H}_{U}$ via virtual process. Up to second order, it gives an effective description of the low-energy physics by a spin-1/2 model. Defining the spin operator $\mathbf{S}_{j}=\frac{1}{2}\hat{a}^{\dagger}_{j\alpha}\bm{\sigma}_{\alpha\beta}\hat{a}_{j\beta}$ with $\bm{\sigma}$ being the Pauli matrices, the effective Hamiltonian reads

[TABLE]

where $J=-\frac{4t^{2}}{\lambda U}$ and $J_{z}=2\lambda-1$. It is an anisotropic $XXZ$ Heisenberg model with DM interaction in transverse field. The Heisenberg coupling and DM interaction proportional to $J$ are induced by flux dependent nearest-neighboring hopping; The transverse field proportional to $\Omega$ is induced by on-site pseuudospin flip. We note that the fermionic counterpart is of the same form as Eq. (3), but with $J=\frac{4t^{2}}{U_{\uparrow\downarrow}}$ and $J_{z}=1$, where the magnetic transitions is a gapless LL to ferromagnetic phase transition Sun et al. (2013).

In the following, we study the GS properties of Eq. (3). We restrict our discussion in the regime $\phi\in[0,\pi]$, since it is the realizable range in the experiments, and Eq. (3) satisfies $\hat{\Gamma}\hat{H}_{\mathrm{eff}}\left(\phi,\lambda,\Omega\right)\hat{\Gamma}^{\dagger}=\hat{H}_{\mathrm{eff}}\left(-\phi,\lambda,\Omega\right)$ with $\hat{\Gamma}=\prod_{j}2\hat{S}_{j}^{x}$, in the theoretical perspective. Besides, we scale the energy by $\frac{4t^{2}}{U}$ and define $\Omega^{\prime}=\Omega/\frac{4t^{2}}{U}$.

III Matrix product state method for searching the low-energy states

The general Hamiltonian Eq. (3) for arbitrary $\phi,\lambda$ and $\Omega$ are not exactly solvable. In this section, we apply the MPS algorithm to determine the GS and low-lying excited states under open boundary condition. It makes use of the Schmidt decomposition (SD) and treats the states and operators in the matrix product form. By discarding the related states of small-weighted singular values, the state space is reduced in block, and states are approximated in some optimal way. Through a variational search, the minimum-energy state in the reduced space is found. In the following, we will introduce how to represent the states and operators in the matrix product form, and display the variational algorithm for determining the low-energy states.

III.1 Schmidt decomposition

To express the states of the system as matrix product form, one makes use of the SD. Any pure state on a composite system $\mathcal{H}_{A}\otimes\mathcal{H}_{B}$ is read as $\left|{\psi}\right\rangle=\sum_{i_{A},i_{B}}^{N_{A},N_{B}}M_{i_{A},i_{B}}\left|{i_{A}}\right\rangle\left|{i_{B}}\right\rangle$, where ${\left|{i_{A}}\right\rangle}$ and ${\left|{i_{B}}\right\rangle}$ are the bases of subsystem $A$ and $B$ with dimension $N_{A}$ and $N_{B}$ respectively. The SD on $\left|{\psi}\right\rangle$ means it can be decomposed as $\left|{\psi}\right\rangle=\sum_{\alpha=1}^{\chi}S_{\alpha}\left|{\alpha_{A}}\right\rangle\left|{\alpha_{B}}\right\rangle$, where $\{\left|{\alpha_{A}}\right\rangle\}$ and $\{\left|{\alpha_{B}}\right\rangle\}$ are the eigenstates of the reduced density matrices $\hat{\rho}_{A}$ and $\hat{\rho}_{B}$ of the subsystem $A$ and $B$ respectively, and $S_{\alpha}^{2}$ are their shared eigenvalules, satisfying $\sum_{\alpha=1}^{\chi}S^{2}_{\alpha}=1$, with $\chi=\min(N_{A},N_{B})$. The set of $\{S_{\alpha}\}$ are referred as Schmidt coefficients, and the number of nonzero Schmidt coefficients $\chi_{s}$, which meets $1\leq\chi_{s}\leq\chi$, is referred as Schmidt rank. To relate the SD with the coefficient matrix $M$, one can expand $\{\left|{\alpha_{A}}\right\rangle\}$ and $\{\left|{\alpha_{B}}\right\rangle\}$ in the original bases as: $\left|{\alpha_{A}}\right\rangle=\sum_{i_{A}}U_{i_{A},\alpha}\left|{i_{A}}\right\rangle$ and $\left|{\alpha_{B}}\right\rangle=\sum_{i_{B}}V^{\dagger}_{\alpha,i_{B}}\left|{i_{B}}\right\rangle$, with matrices $U$ and $V$ satisfying $U^{\dagger}U=1$ and $V^{\dagger}V=1$. Defining a diagonal matrix $S$ with entries $S_{\alpha}$, $M$ can be expanded as $M_{i_{A},i_{B}}=\sum_{\alpha}U_{i_{A},\alpha}S_{\alpha}V^{\dagger}_{\alpha,i_{B}}$, which is just a singular value decomposition (SVD) of the $M$ matrix.

The form of the SD and the hint of the SVD on coefficient matrix $M$ provide an optimal way to approximate the state vector $\left|{\psi}\right\rangle$ with smaller spanned dimension. That is, if the Schmidt coefficients $S_{\alpha}$ is arranged in a descending order: $S_{1}\geq S_{2}\geq\cdots$, and the Schmidt rank is truncated to some smaller $\tilde{\chi}<\chi_{s}$ by discarding the states with small-weighted singular value, the state

[TABLE]

is the closest rank-$\tilde{\chi}$ approximation to $\left|{\psi}\right\rangle$ in the sense that the Frobenius norm between the coefficient matrices of these two states is minimized. This property is the key ingredient for the feasibility of MPS algorithm, as can be seen in the following.

III.2 Matrix product states

Now, we show how to represent the quantum states of our system in matrix product form via making use of the SD or SVD. Any pure state in 1D can be written as $\left|{\psi}\right\rangle=\sum^{d}_{\sigma_{1}\cdots\sigma_{L}}c_{\sigma_{1}\cdots\sigma_{L}}\left|{\sigma_{1}\cdots\sigma_{L}}\right\rangle$, where $\{\left|{\sigma_{i}}\right\rangle\}$ called the physical indices, are the local basis with dimension $d$, $L$ is the number of lattice sites and $c_{\sigma_{1}\cdots\sigma_{L}}$ is the complex amplitude. This pure state can be represented as an MPS

[TABLE]

where $\vec{\sigma}$ is shorted for $\{\sigma_{i}\}$, following the recursive routine:

[TABLE]

where the initial $\Psi$ is reshaped from the coefficient vector $\Psi_{a_{0}\sigma_{1},\sigma_{2}\cdots\sigma_{L}}=c_{\sigma_{1}\cdots\sigma_{L}}$ (here $a_{0}=1$ is an auxiliary index), the first equality is the SVD on $\Psi$, the second equality is just the replacement of matrix $U$ by the tensor $A^{[l]}$, and the last equality is the reshaping of $SV^{\dagger}$ into a new matrix $\Psi$. Each $A^{[l]}$ consists of $d$ matrices of bond dimension $\chi_{l-1}\times\chi_{l}$, which is determined from the SVD: $\chi_{l}=\min\left(d^{l},d^{L-l}\right)$. The property $U^{\dagger}U=I$ makes the $A^{[l]}$ satisfy the normalization relationship

[TABLE]

The MPS with all matrices satisfying this normalization condition is called left-canonical MPS. Note that the decomposition of the complex amplitudes is not unique. If the recursive procedure is started from the right side:

[TABLE]

where the initial $\Psi_{\sigma_{1}\cdots\sigma_{L-1},\sigma_{L}a_{L}}=c_{\sigma_{1}\cdots\sigma_{L}}$ (here $a_{L}=1$ is the auxiliary index), the first equality is the SVD on matrix $\Psi$, the second equality is replacing $V^{\dagger}$ by tensor $B^{[l]}$ and the last equality is the reshaping of $US$ into a new $\Psi$, the MPS reads as

[TABLE]

The tensor $B^{[l]}$ consists of $d$ matrices of bond dimension $\chi_{l-1}\times\chi_{l}$ and satisfy the normalization condition

[TABLE]

making use of the fact $V^{\dagger}V=I$ on each SVD decomposition. The MPS with all matrices satisfying the above normalization condition is called right-canonical MPS.

In fact, the degree of nonuniqueness is much higher: there is a gauge degree of freedom in writing the MPS. That is, if one inserts an invertible matrix $X$ with dimension $\chi_{l}\times\chi_{l}$ and its inverse $X^{-1}$ into two adjacent MPS matrices $M^{\sigma_{l}[l]}$ and $M^{\sigma_{l+1}[l+1]}$ and makes the transformation $M^{\sigma_{l}[l]}X\to M^{\sigma_{l}[l]},X^{-1}M^{\sigma_{l+1}[l+1]}\to M^{\sigma_{l+1}[l+1]}$, the MPS is invariant. Instructively, one can specify a general MPS $\left|{\psi}\right\rangle=\sum_{\vec{\sigma}}M^{\sigma_{1}[1]}\cdots M^{\sigma_{L}[L]}\left|{\vec{\sigma}}\right\rangle$, by choosing a site $k$, which is called the orthogonal center, that all the matrices left and right to it are left- and right-normalized respectively. This particular kind of MPS is called mixed-canonical MPS. The left and right-normalization condition can be imposed by the way quite similar to the one constructing MPS from the coefficient vector. The left-normalization condition is imposed by the recursive routine:

[TABLE]

starting from $M^{[1]}$ and stoping before reaching $M^{[k]}$. The first equality is the reshaping of tensor $M^{[l]}$ into a matrix $\tilde{M}$, the second equality is the SVD of $\tilde{M}$, the third equality is the replacement of $U^{[L]}$ by the tensor $A^{[l]}$, and the last equality is the absorption of the matrices $S^{[\mathrm{L}]}$ and $V^{[\mathrm{L}]\dagger}$ into the next $M^{[l+1]}$. Then, the tensors $A^{[l]}$ satisfy the normalization condition Eq. (III.2). The bond dimension of each tensor $A^{[l]}$ is determined from the SVD: $\chi^{\prime}_{l}=\min(d\chi_{l-1},\chi_{l})$, where $\chi_{l}$ is the bond dimension of the original tensor $M^{[l]}$. The right-normalization condition is enforced by the recursion routine:

[TABLE]

starting from $M^{[L]}$ and stoping before reaching $M^{[k]}$. Then, the tensors $B^{[l]}$ satisfy the normalization condition Eq. (III.2). The bond dimension of each tensor $B^{[l]}$ is determined from the SVD: $\chi^{\prime}_{l}=\min(d\chi_{l+1},\chi_{l})$. At last, multiply the residual $U$ and $S$ matrices resulting from the two recursion to $M^{[k]}$: $\tilde{M}^{\sigma_{k}[k]}=(S^{[\mathrm{L}]}V^{[\mathrm{L}]\dagger})M^{\sigma_{k}[k]}(U^{[\mathrm{R}]}S^{[\mathrm{R}]})$, and the norm square of the state is read as $\left\langle{{\psi}}\mathrel{\left|{\vphantom{{\psi}{\psi}}}\right.\kern-1.2pt}{{\psi}}\right\rangle=\sum_{\sigma_{k}}\mathrm{Tr}\left(\tilde{M}^{\sigma_{k}[k]\dagger}\tilde{M}^{\sigma_{k}[k]}\right)$. Definitely, if $k$ is set as $1$ or $L$, the recursion Eq. (III.2) or Eq. (III.2) gives the left- or right-canonical MPS.

The MPSs obtained in theses ways are exact, but not productive for computation, for the reason that the dimension of the matrices grows up exponentially, as can be seen from the recursion constructing the canonical MPS. One way to make the MPSs practicable is to bound the bond dimension to some maximum $\tilde{\chi}$ following Eq. (4). That is, in the process building a canonical MPS from a state vector or a general MPS, once the bond dimension grows above $\tilde{\chi}$, truncate it to $\tilde{\chi}$ following Eq. (4). As a result, the elements in the MPS are decimated in block effectively. This approximation is valid for GS in 1D without losing noticeable accuracy, due to two facts: the singular value spectra decay exponentially; the bipartite entanglement of the GS obeys an area law in the gapped phase and increases as subsystem size only logarithmically near the critical point. The MPS algorithm which introduced in the following is in fact based on such a decimation procedure.

III.3 Matrix product operators

The natural generation of writing states as matrix product form to operators is the matrix product operator (MPO). A general operator $\hat{O}$ expressed in the local bases is $\hat{O}=\sum_{\vec{\sigma},\vec{\sigma}^{\prime}}O_{\vec{\sigma},\vec{\sigma}^{\prime}}\left|{\vec{\sigma}}\right\rangle\left\langle{\vec{\sigma}^{\prime}}\right|$. Its matrix product form is defined as

[TABLE]

where each $\hat{W}^{[l]}$ can be considered as a $D_{l-1}\times D_{l}$ operator-valued matrix with elements $\hat{W}^{[l]}_{b_{l-1},b_{l}}=\sum_{\sigma_{l},\sigma^{\prime}_{l}}W^{\sigma_{l},\sigma^{\prime}_{l}[l]}_{b_{l-1},b_{l}}\left|{\sigma_{l}}\right\rangle\left\langle{\sigma^{\prime}_{l}}\right|$. The expression of the MPO is actually a sum of matrices products, which is of the same form as a general Hamiltonian. This makes it quite intuitive to express a 1D Hamiltonian as an MPO. In fact, through defining some finite state automaton rules, all 1D Hamiltonian with finite-range interaction can be written as exact MPO form Wall and Carr (2012). In considering our model Eq. (3), we can write down its MPO representation directly:

[TABLE]

III.4 Variational ground state search

We now show how to obtain the GS using an MPS as a variational ansatze. To find the optimal ground MPS $\left|{\psi}\right\rangle$ with maximum bond dimension $\chi$, we have to minimize the functional

[TABLE]

where $E$ is the Lagrangian multiplier, and $\hat{H}$ is in the MPO form. This optimization problem is hard to solve at the first glace for that the variables appear as products. Fortunately, we can get the optimal solution via an iterative algorithm: minimize the energy $\varepsilon$ with respect to the tensor $M^{[k]}$ at site $k$ with all other MPS tensors fixed, and obtain the better state lower in energy; move to the next $M^{[k+1]}$ and find the state again lower in energy; repeat sweeping through all sites until the energy is converged, and finally the minimum energy and the corresponding GS are obtained.

To minimize the energy functional $\varepsilon$ with respect to a particular $M^{[k]}$, we have to calculate $\varepsilon$ explicitly. Suppose the ansatze MPS $\left|{\psi}\right\rangle$ is of mixed-canonical form with the orthogonal center at a chosen $k$, the overlap can be directly read as $\left\langle{{\psi}}\mathrel{\left|{\vphantom{{\psi}{\psi}}}\right.\kern-1.2pt}{{\psi}}\right\rangle=\sum_{\sigma_{k}a_{k-1}a_{k}}|M^{\sigma_{k}[k]}_{a_{k-1},a_{k}}|^{2}$. The expectation value of the MPO $\hat{H}$ in $\left|{\psi}\right\rangle$ is written as

[TABLE]

where the tensors $L$ and $R$ are the partial overlap of the Hamiltonian and the state, constructed following the recursive procedure

[TABLE]

with the initial $L_{a_{0},a^{\prime}_{0}}^{b_{0}[0]}=R_{a_{L},a^{\prime}_{L}}^{b_{L}[L+1]}=1$. Now take the extremum of Eq. (28) with respect to $M^{\sigma_{k}[k]*}_{a_{k-1},a_{k}}$, we obtain

[TABLE]

It is an eigenvalue problem $\hat{H}^{[k]}\left|{v^{[k]}}\right\rangle-E\left|{v^{[k]}}\right\rangle=0$ if we view $M^{[k]}$ as a vector $\left|{v^{[k]}}\right\rangle$ with entries $v^{[k]}_{a_{k-1}\sigma_{k}a_{k}}=M^{\sigma_{k}[k]}_{a_{k-1},a_{k}}$, and introduce the effective Hamiltonian by the reshaping

[TABLE]

Thus the optimal solution $M^{[k]}$ at present can be obtained by solving the effective Hamiltonian $\hat{H}^{[k]}$ for the GS $\left|{v^{[k]}_{0}}\right\rangle$ with energy $E_{0}$ and reshaping $\left|{v^{[k]}_{0}}\right\rangle$ back to $M^{[k]}$, with $E_{0}$ being the current energy.

In summary, the iterative variational GS search algorithm is as follows:

$\left(i\right)$ Input. Input $\hat{H}$ in the MPO form, a guessed MPS $\left|{\psi}\right\rangle$ with maximum bond dimension $\chi$, and a tolerance $\varsigma$ for energy convergence.

$\left(ii\right)$ Initialization. Transform $\left|{\psi}\right\rangle$ to the right-canonical form according to Eq. (III.2). Initialize the [math]th tensor $L_{a_{0},a^{\prime}_{0}}^{b_{0}[0]}=1$. Construct all the right overlaps $R$ by Eq. (III.4).

$\left(iii\right)$ Right sweep. Construct the effective Hamiltonian according to Eq. (32), solve it for the minimum energy $E_{0}$ and state vector $\left|{v^{[k]}}\right\rangle$. Update $M^{[k]}$ by reshaping $M^{\sigma_{k}[k]}_{a_{k-1},a_{k}}=v^{[k]}_{a_{k-1}\sigma_{k}a_{k}}$. Left-normalize $M^{[k]}$ and move the orthogonal center to the right site $k+1$ by Eq. (III.2). Update the $k$th overlap $L^{[k]}$ recursively following Eq. (III.4). Continue sweeping to the right until the boundary is reached.

$\left(iv\right)$ Left sweep. Construct the effective Hamiltonian and solve it for the minimum energy $E_{0}$ and state vector $\left|{v^{[k]}}\right\rangle$. Update $M^{[k]}$ by reshaping $M^{\sigma_{k}[k]}_{a_{k-1},a_{k}}=v^{[k]}_{a_{k-1}\sigma_{k}a_{k}}$. Right-normalize $M^{[k]}$ and move the orthogonal center to the left site $k-1$ by Eq. (III.2). Update the $k$th overlap $R^{[k]}$ recursively following Eq. (III.4). Continue sweeping to the left until the boundary is reached.

$\left(v\right)$ Repeat steps $\left(iii\right)$ and $\left(iv\right)$ until the convergence is achieved $\langle\hat{H}^{2}-E_{0}^{2}\rangle<\varsigma$.

$\left(vi\right)$ Output. Output the minimum energy $E_{0}$ and the MPS $\left|{\psi}\right\rangle$ which is of right-canonical form now.

Note that the state obtained in this way is not necessarily the GS, for it may get stuck in some local minimum state. Two ways help improving such a dilemma. The first is to prepare the initial state in the desired subspace with good quantum number. The approached state must be the energy-minimized state in that subspace. This is not the case for our model, for there is no explicit conserved quantity in our system. The second way is to generalize the single site to a contiguous block during the local search and modify the algorithm accordingly, at the cost of consuming longer computational time and more computational resources.

III.5 Variational excited states search

With the GS achieved, we now show how to obtain the subsequent excited states incrementally. As the way searching the GS, the $n$th excited state ($n$ExS) is found by minimizing the energy functional $\varepsilon_{n}\left[\left|{\psi_{n}}\right\rangle\right]=\left\langle{\psi_{n}}\right|\hat{H}\left|{\psi_{n}}\right\rangle-E_{n}\left\langle{{\psi_{n}}}\mathrel{\left|{\vphantom{{\psi_{n}}{\psi_{n}}}}\right.\kern-1.2pt}{{\psi_{n}}}\right\rangle$, but under $n$ orthogonality constraints

[TABLE]

where $\{\left|{\phi_{m}}\right\rangle\}$ with $m=0,1,\cdots n-1$ are the $n$ lower-lying eigenstates. It makes the searching program constrained in the space orthogonal to the one formed by $\{\left|{\phi_{m}}\right\rangle\}$. The minimization of $\varepsilon_{n}$ with respect to local MPS tensor $M^{[k]}$ under such constrainets is equivalent to solving the eigenvalue problem

[TABLE]

where $\hat{H}^{[k]}$, defined by Eq. (32), is the effective Hamiltonian for the variational local tensor $M^{[k]}$, and $\hat{P}^{[k]}$ is the project operator into the orthogonal space of the lower-lying space $\{\left|{\phi_{m}}\right\rangle\}$.

To find the projector $\hat{P}^{[k]}$ for every local tensor, we calculate the overlaps between the lower-lying states and the variational state explicitly: $\left\langle{{\psi_{n}}}\mathrel{\left|{\vphantom{{\psi_{n}}{\phi_{m}}}}\right.\kern-1.2pt}{{\phi_{m}}}\right\rangle=\sum_{a_{k-1},\sigma_{k},a_{k}}M_{a_{k-1},a_{k}}^{\sigma_{k}[k]*}F_{a_{k-1},a_{k}}^{\sigma_{k}[k](m)}$, with

[TABLE]

where $A^{[k](m)}$ is the $k$th MPS tensor of the $m$th lower-lying state, and the tensors $\mathcal{L}^{(m)}$ and $\mathcal{R}^{(m)}$ are the partial overlap between $\left|{\phi_{m}}\right\rangle$ and $\left|{\psi_{n}}\right\rangle$, constructed following the recursive procedure

[TABLE]

with the initial $\mathcal{L}^{[0](m)}_{a_{0},a^{\prime}_{0}}=1$ and $\mathcal{R}^{[L+1](m)}_{a_{L},a^{\prime}_{L}}=1$. By viewing the tensor $M^{[k]}\left(F^{[k](m)}\right)$ as vector $\left|{v^{[k]}}\right\rangle\left(\left|{u^{[k](m)}}\right\rangle\right)$ with elements $v^{[k]}_{a_{k-1}\sigma_{k}a_{k}}=M^{\sigma_{k}[k]}_{a_{k-1},a_{k}}\left(u^{[k](m)}_{a_{k-1}\sigma_{k}a_{k}}=F_{a_{k-1},a_{k}}^{\sigma_{k}[k](m)}\right)$, the projector $\hat{P}^{[k]}$ for the local tensor space is directly read as

[TABLE]

where $\left(\mathcal{N}^{-1}\right)_{mm^{\prime}}=\mathrm{Tr}\left(F^{[k](m)\dagger}F^{[k](m^{\prime})}\right)$.

The flow of the variational $n$ExS search is the same as that for the GS, but with each step modified accordingly:

$\left(i\right)$ Input. Input $\hat{H}$ in the MPO form, the $n$ obtained lower-lying states $\{\left|{\phi_{m}}\right\rangle\}$ in the MPS form, a guessed MPS $\left|{\psi}\right\rangle$ for the $n$ExS with maximum bond dimension $\chi$, and a tolerance $\varsigma$ for energy convergence.

$\left(ii\right)$ Initialization. Transform $\left|{\psi}\right\rangle$ to the right-canonical form according to Eq. (III.2). Initialize the [math]th tensors $L_{a_{0},a^{\prime}_{0}}^{b_{0}[0]}=1$ and $\mathcal{L}_{a_{0},a^{\prime}_{0}}^{[0](m)}=1$. Construct all the right overlaps $R$ and $\mathcal{R}$ by Eq. (III.4) and Eq. (III.5) respectively.

$\left(iii\right)$ Right sweep. Construct the effective Hamiltonian by Eq. (32) and the projector by Eq. (35) and Eq. (37). Solve the projected effective Hamiltonian for the minimum energy $E_{n}$ and state vector $\left|{v^{[k]}}\right\rangle$. Update $M^{[k]}$ by reshaping $M^{\sigma_{k}[k]}_{a_{k-1},a_{k}}=v^{[k]}_{a_{k-1}\sigma_{k}a_{k}}$. Left-normalize $M^{[k]}$ and move the orthogonal center to the right site $k+1$ by Eq. (III.2). Update the $k$th overlaps $L^{[k]}$ and $\mathcal{L}^{k}$ recursively following Eq. (III.4) and Eq. (III.5) respectively. Continue sweeping to the right until the boundary is reached.

$\left(iv\right)$ Left sweep. Construct the projected effective Hamiltonian and solve it for the minimum energy $E_{n}$ and state vector $\left|{v^{[k]}}\right\rangle$. Update $M^{[k]}$ by reshaping $M^{\sigma_{k}[k]}_{a_{k-1},a_{k}}=v^{[k]}_{a_{k-1}\sigma_{k}a_{k}}$. Right-normalize $M^{[k]}$ and move the orthogonal center to the left site $k-1$ by Eq. (III.2). Update the $k$th overlaps $R^{[k]}$ and $\mathcal{R}^{[k]}$ recursively following Eq. (III.4) and Eq. (III.5) respectively. Continue sweeping to the left until the boundary is reached.

$\left(v\right)$ Repeat steps $\left(iii\right)$ and $\left(iv\right)$ until the convergence is achieved $\langle\hat{H}^{2}-E_{n}^{2}\rangle<\varsigma$.

$\left(vi\right)$ Output. Output the energy $E_{n}$ and the MPS $\left|{\psi}\right\rangle$ which is of right-canonical form now.

The variational ES search also suffers form the local minimum dilemma. Besides, the area-law of entanglement does not apply to the bulk excited states, thus the bond dimension $\chi$ has to be increased to ensure the discarded states are of small-weighted singular values. This limits the algorithm to be applicable only for low-energy states.

IV Magnetic phase transitions

In order to distinguish between different phases, we calculate the energy gaps $\Delta_{1,2}=E_{1,2}-E_{0}$ between the GS and the first two excited states, and expectation value of observables in the GS: the two-site correlation function $\langle\hat{S}_{j}^{\alpha}\hat{S}_{l}^{\alpha}\rangle$ and the associated structure factor $Q_{\alpha}(k)=\frac{1}{L}\sum_{jl}e^{ik(j-l)}\langle\hat{S}_{j}^{\alpha}\hat{S}_{l}^{\alpha}\rangle$; the order parameters, magnetization $M_{\alpha}=\frac{1}{L}|\sum_{j}\langle\hat{S}_{j}^{\alpha}\rangle|$, staggered magnetization $N_{\alpha}=\frac{1}{L}|\sum_{j}(-1)^{j}\langle\hat{S}_{j}^{\alpha}\rangle|$, and the spiral order $C_{\alpha}=\frac{1}{L}|\sum_{j}\langle[\bm{\vec{S}}_{j}\times\bm{\vec{S}}_{j+1}]^{\alpha}\rangle|$, where $\alpha=x,y,z$.

The calculation of expectation value of observables in the framework of MPS is quite feasible. For a general observable $\hat{O}=\sum_{\vec{\sigma},\vec{\sigma}^{\prime}}O_{\sigma_{1},\sigma^{\prime}_{1}}^{[1]}\cdots O_{\sigma_{L},\sigma^{\prime}_{L}}^{[L]}\left|{\vec{\sigma}}\right\rangle\left\langle{\vec{\sigma}^{\prime}}\right|$, the expectation value in the state expressed as right-orthogonal MPS is directly written as

[TABLE]

where the dominate is the norm square of the state. When reduced to single-site and two-site cases, the calculation is greatly simplified by making use of the normalization conditions. In calculating the single-site observable $\hat{S}_{j}^{\alpha}=\sum_{\sigma_{j},\sigma^{\prime}_{j}}S^{\alpha[j]}_{\sigma_{j},\sigma^{\prime}_{j}}\left|{\sigma_{j}}\right\rangle\left\langle{\sigma^{\prime}_{j}}\right|$, it is convinent to transform the right-canonical MPS into one of mixed-canonical form with the orthogonal center at $j$. Then the expectation value is just read as

[TABLE]

In calculating the expectation value of the two-site observable $\hat{S}_{j}^{\alpha}\hat{S}_{l}^{\beta}$ (here we set $j<l$ without lose of generality), we move the orthogonal center $k$ to any site in the range $j\leq k\leq l$. The evaluation of the expectation value is thus reduced to contracting the tensors in this range. This can be done in a recursive manner: we first construct a tensor at site $j$, $G^{[j]}_{a_{j},a^{\prime}_{j}}=\sum_{\sigma_{j},\sigma^{\prime}_{j},a_{j-1}}M^{\sigma_{j}[j]*}_{a_{j-1},a_{j}}S^{\alpha[j]}_{\sigma_{j},\sigma^{\prime}_{j}}M^{\sigma^{\prime}_{j}[j]}_{a_{j-1},a^{\prime}_{j}}$, then generate the next $G$ by the recursion

[TABLE]

until the $l$th site is reached, where $O^{[m]}$ is the identity matrix when $m<l$ and $O^{[l]}=S^{\beta[l]}$, and finally the expectation value is calculated as $\langle\hat{S}^{\alpha}_{j}\hat{S}^{\beta}_{l}\rangle=\frac{\mathrm{Tr}\left(G^{[l]}\right)}{\mathrm{Tr}\left(M^{\sigma_{k}\dagger[k]}M^{\sigma_{k}[k]}\right)}$.

In our MPS calculation of the GS and the expectation value of observables, we set the system size as $L=195$, unless otherwise specified, and the tolerance for energy convergence as $\varsigma=10^{-7}$. We randomly initialize an MPS ansatz with maximum bond dimension for several $\chi$ and searching for the GS and energy following the algorithm introduced in Sec. III.4. We then calculate the expectation value of observables in the GS, and find the results already converges as $\chi$ approaches $16$ for $L=195$, in the sense that increasing $\chi$ does not change the results of GS energy and the expectation value of observables. In calculating the energy gaps between the GS and the first two excited states, we fix the system size $L=45$, $75$, $105$, $135$, $165$, $195$ and $295$ respectively. We find the maximum bond dimension has to reach $24$ for results convergence.

According to the behavior of the energy gap and correlation function, we divide the GS into four typical phases with long range correlation: (I) the $z$-FM phase with ferromagnetic correlation along the $\hat{z}$ direction; (II) the $x$-PARA phase only with ferromagnetic correlation along the external filed direction; (III) the $y$-AFM phase with antiferromagnetic correlation along the $\hat{y}$ direction; (IV) the $xy$-SP phase with spiral correlation on the $\hat{x}\hat{y}$ plane.

The imposing of SOC introduces transverse field which polarizes the spin along the $\hat{x}$ axis, and DM interaction which causes spin to rotate on the $\hat{x}\hat{y}$ plane. In the following, we discuss the interplay between the anisotropic Heisenberg coupling, DM interaction and transverse field.

IV.1 Interplay between anisotropy and DM interaction

We first study the interplay between the anisotropy and DM interaction by setting $\Omega^{\prime}=0$, i.e. without Ramman-assisted tunneling. Intuitively, the effect of the SOC should be absent and the GS properties not depend on the flux $\phi$. To see this, we gauge away the DM interaction by performing a unitary transformation $\hat{H}^{\prime}_{\mathrm{eff}}=\hat{U}\hat{H}_{\mathrm{eff}}\hat{U}^{\dagger}$ with $\hat{U}=\prod_{j}e^{i\phi j\hat{S}_{j}^{z}}$, resulting in

[TABLE]

It is a standard ferromagnetic XXZ model, whose GS has been solved exactly by the Bethe ansatz Takahashi (2005). As a result, when $\lambda>1$, the GS is a gapped $z$-FM phase with long-range ferromagnetic correlation in the $\hat{z}$ direction. When $\lambda<1$, it is in a gapless phase with algebraic decaying correlation functions, signaling a gapless LL phase. Regarding the original Hamiltonian, the magnetic flux $\phi$ makes the phase finer sorted. When $\lambda>1$, the ground state is in the $z$-FM phase regardless any $\phi$. When $\lambda<1$, the gapless LL phase transforms to a gapless ferromagnet on the $\hat{x}\hat{y}$ plane for $\phi=0$, a gapless antiferromagnet on the $\hat{x}\hat{y}$ plane for $\phi=\pi$, and a gapless spiral phase for $0<\phi<\pi$. To illustrate the nature of the gapless spiral phase, we plot in Fig. 1 the finite-size scaling of the energy gaps $\Delta_{1,2}$, the correlation functions and structure factors for the representative parameter point $(\phi,\lambda)=(0.5\pi,0.75)$. The gapless spiral phase is characterized by zero energy gap shown in Fig. 1(a) and oscillated correlation functions along the $\hat{x}$ and $\hat{y}$ direction with algebraic decaying envelop shown in Fig. 1(b). What’s more, the periods of the oscillation of $\langle\hat{S}_{j}^{x}\hat{S}_{l}^{x}\rangle$ and $\langle\hat{S}_{j}^{y}\hat{S}_{l}^{y}\rangle$ are the same, as shown in Fig. 1(c), where $Q_{x}(k)$ and $Q_{y}(k)$ show the same peak at $k=\phi$. The gapless spiral phase found here is the same as that of the SOC-induced nearest-neighbor spin-flip model and the fermioic ladder, both of which belong to the LL phase.

IV.2 Interplay between anisotropy and external field

To discuss the interplay between anisotropic Heisenberg coupling and the external field, we set $\phi=0$ or $\phi=\pi$, i.e. zero or maximum magnetic flux, at which the DM interaction is zero. In these two limits, the Hamiltonian is reduced to a(n) (anti)ferromagnetic XXZ model in transverse field:

[TABLE]

where ferromagnetic coupling corresponds to $\phi=0$ and antiferromagnetic coupling corresponds to $\phi=\pi$. The phase diagram of these two models are presented in Fig. 2, in which besides the $z$-FM (I) and gapless LL phase (red solid lines), additional $x$-PARA phase (II) and $y$-AFM phase (III) appear, which will be explained in the following.

We first consider the ferromagnetic case, i.e. $\phi=0$. When $\lambda>1$, the coupling in the $\hat{z}$ direction dominates over that on the $\hat{x}\hat{y}$ plane. The GS is mainly determined by the competition between the coupling in the $\hat{z}$ direction and the transverse field. There is an Ising phase transition as $\Omega^{\prime}$ increases: when $\Omega^{\prime}$ is smaller than a threshold $\Omega^{\prime}_{c}$, the GS exhibits long range ferromagnetic correlation along the $\hat{z}$ direction, so is in the $z$-FM phase. When $\Omega^{\prime}$ increases beyond $\Omega^{\prime}_{c}$, the $z$-ferromagnetic correlation vanishes completely and only long range $x$-ferromagnetic correlation remains, thus the GS is in the so-called $x$-PARA phase. Both of the $z$-FM and $x$-PARA phase are illustrated in Fig. 3 (a1)-(a3) and (b1)-(b3) respectively, in which we plot the energy gaps, the correlation function $\langle S^{\alpha}_{j}S^{\alpha}_{l}\rangle$ and the related structure factors $Q_{\alpha}(k)$ for $(\lambda,\Omega^{\prime})=(1.25,0.2)$ and $(1.25,0.3)$ respectively. In the $z$-FM phase, the GS is twofold degenerate shown in Fig. 3(a1). Both $\langle S^{z}_{j}S^{z}_{l}\rangle$ and $\langle S^{x}_{j}S^{x}_{l}\rangle$ are finite for any range and their structure factors have a peak at $k=0$, see Fig. 3(a2) and (a3). In the $x$-PARA phase, there is a finite gap between the GS and the excited states, shown in Fig. 3(b1). The correlation functions $\langle S^{y}_{j}S^{y}_{l}\rangle$ and $\langle S^{z}_{j}S^{z}_{l}\rangle$ decay to zero rapidly and only $\langle S^{x}_{j}S^{x}_{l}\rangle$ is finite in the long range shown in Fig. 3(b2), which is signatured by the smoothness of $Q_{y}(k)$ and $Q_{z}(k)$ and the sharp peak of $Q_{x}(k)$ at $k=0$, see Fig. 3(b3). In fact, the $z$-FM-$$x-PARA phase transition is not unique for $\phi=0$, but also applies for arbitrary value of $\phi$ with $\lambda>1$. In all the following discussions, we will not refer to the case with $\lambda>1$ to avoid repetition. When $\lambda<1$, the GS is a gapless ferromagnet on the $\hat{x}\hat{y}$ plane when $\Omega^{\prime}=0$, as discussed in Sec. IV.1. Nonzero $\Omega^{\prime}$ breaks the degeneracy and polarizes the GS along the $\hat{x}$ direction immediately, meaning that there is no finite critical point between the gapless phase and the $x$-PARA phase. We summarize the analysis above in the phase diagram for $\phi=0$ shown in Fig. 2(a).

We then consider the antiferromagnetic case, i.e. $\phi=\pi$. When $\lambda=1$, via rotating every second spin around the $\hat{z}$ axis by an angle $\pi$, the Hamiltonian can be transformed into an isotropic ferromagnetic Heisenberg chain in staggered field, that is $\hat{H}^{\prime}_{\mathrm{eff}}=-\sum_{j}\vec{\bm{S}}_{j}\cdot\vec{\bm{S}}_{j+1}-\Omega^{\prime}\sum_{j}\left(-1\right)^{j}\hat{S}_{j}^{x}$. The GS remains a gapless LL phase up to some critical point $\Omega^{\prime}_{c}$, beyond which, it is a fully polarized antiferromagnetic phase Alcaraz and Malvezzi (1995); Dmitriev et al. (2002). They translate to a gapless antiferromagnet and the $x$-PARA phase respectively for the original model. When $\lambda<1$, the GS is a gapless antiferromagnet on the $\hat{x}\hat{y}$ plane when $\Omega^{\prime}=0$, as discussed in Sec. IV.1. When the transverse field is turned on, its incompatibility with the antiferromagnetic coupling along the $\hat{x}$ direction makes the GS favor antiferromagnetic order along the $\hat{y}$ direction. To manifest the nature of such a phase, we plot in Fig. 3(c1)-(c3) the energy gaps, the correlation functions and corresponding structure factors for a representative parameter point $(\lambda,\Omega^{\prime})=(0.75,1)$. The GS are twofold degenerate shown in Fig. 3(c1); the correlation function $\langle\hat{S}_{j}^{y}\hat{S}_{l}^{y}\rangle$ displays clear long range antiferromagnetic structure shown in Fig. 3(c2), which is also confirmed by the sharp peak of $Q_{y}(k)$ at $k=\pi$ shown in Fig. 3(c3). Such a antiferromagnetic structure along the $\hat{y}$ direction is accompanied with a finite ferromagnetic correlation along the $\hat{x}$ direction, shown in Fig. 3(c2) and (c3), where $\langle\hat{S}_{j}^{x}\hat{S}_{l}^{x}\rangle$ is finite and $Q_{x}(k)$ has a smaller peak at $k=0$. Anyway, we call such a twofold degenerate state with long range antiferomagnetic correlation in the $\hat{y}$ direction the $y$-AFM phase. As $\Omega^{\prime}$ increases further, the antiferromagnetic correlation breaks down and only the $x$-ferromagnetic correlation survives, and thus the GS enters the $x$-PARA phase. There is a finite critical point separating the $y$-AFM phase and the $x$-PARA phase. Our numerical phase diagram for $\phi=\pi$ is shown in Fig. 2(b), consistent with the result obtained by mean-field approach Dmitriev et al. (2002, 2002) if taking care of the definition of parameters.

IV.3 Interplay between DM interaction and external field

We then explore the interplay between the DM interaction and the transverse field by fixing $\phi=\pi/2$, at which the $x$-$x$ and $y$-$y$ coupling strengths are zero. In this case, the Hamiltonian is reduced to

[TABLE]

When $\lambda<1$, the GS is mainly decided by the DM interaction and the transverse field. The phase diagram when $\phi=\pi/2$ is shown in Fig. 4, in which the $y$-AFM phase (III) is absent, and the $xy$-SP phase (IV) emerges.

The $xy$-SP phase emerges for $0<\Omega^{\prime}<\Omega^{\prime}_{c}$. It differs form the gapless spiral phase at $\Omega^{\prime}=0$. In Fig. 5, we plot the energy gaps, correlation functions and the structure factors for the $xy$-SP phase at a representative parameter point $(\lambda,\Omega^{\prime})=(0.75,0.5)$ with system size $L=295$ and maximum bond dimension $\chi=24$. Although the GS is maintained gapless from the excited states (at least for the first two excited states), see Fig. 5(a), the correlation functions $\langle\hat{S}_{j}^{x}\hat{S}_{l}^{x}\rangle$ and $\langle\hat{S}_{j}^{y}\hat{S}_{l}^{y}\rangle$ show no decaying behaviour even in the long range, see Fig. 5(b). The periodic oscillations of $\langle\hat{S}_{j}^{x}\hat{S}_{l}^{x}\rangle$ and $\langle\hat{S}_{j}^{y}\hat{S}_{l}^{y}\rangle$ are also kept in the $xy$-SP phase, but their periods are not the same anymore, shown in Fig. 5(c), where we see that besides with the shared peak with $Q_{y}(k)$, $Q_{x}(k)$ has another peak at $k=0$. In this way, the $xy$-SP phase also distinguishes itself from the spiral phase in nearest-neighbor spin-flip model induced by SOC. When $\Omega^{\prime}$ increases beyond $\Omega^{\prime}_{c}$, the effect of the transverse field prevails and the GS is in the $x$-PARA phase. In this way, The phase diagram for $\phi=\pi/2$ is shown in Fig. 4.

IV.4 Interplay among anisotropy, DM interaction and external field

Finally, we discuss the general cases when the anisotropic Heisenberg coupling, the DM interaction and the transverse field are nonzero. When $\Omega^{\prime}$ increases, the system may undergo $xy$-SP-$$x-PARA phase transition or transit to the $y$-AFM phase before being in the $x$-PARA phase, which depends on $\phi$. The former is the case when $0<\phi<\pi/2$, within which the $x$-$x$ and $y$-$y$ Heisenberg couplings are ferromagnetic, and only the spiral and $x$-ferromagnetic order may exist. We have checked this argument by calculating the phase diagram by fixing $\phi<\pi/2$ and find there are only quantitative differences from Fig. 4, which is not shown here. When $\pi/2<\phi<\pi$, the $x$-$x$ and $y$-$y$ spin couplings are antiferromagnetic. There might appear $y$-AFM phase when the antiferromagnetic coupling along with the transverse field prevails. Thus, as $\Omega^{\prime}$ increases, the system might undergo $xy$-SP-$$y-AFM-$$x-PARA phase transitions. The phase diagram in the latter case is as in Fig. 6 with fixed $\phi=0.8\pi$.

We now show how we determine the phase boundaries separating the four typical phases according to the behaviour of the order parameters. Other order parameters except the ones related to the four typical phases are zero in all the parameter regime, so we do not show them further. $M_{x}$ is nonzero as far as $\Omega^{\prime}>0$. In the $x$-PARA phase, all the order parameters except $M_{x}$ are zero. In the $z$-FM phase, both $M_{x}$ and $M_{z}$ are non-zero. In the $y$-AFM phase, in addition to non-zero $M_{x}$ and $N_{y}$, the incomplete $y$-AFM structure along with the extra $x$-ferromagnetic structure makes $C_{xy}$ non-zero, but suppressed compared to $M_{x}$ and $N_{y}$. Similarly, in the $xy$-SP phase, apart from $M_{x}$ and $C_{xy}$, there is small amount of $N_{y}$. In Fig. 6, the three red-dashed lines labeled (a), (b) and (c) correspond to $\lambda=1.25$, $0.92$, and $0.75$ respectively. We plot the order parameters as function of $\Omega^{\prime}$ along these lines in Fig. 7. In Fig. 7(a), there is a first-order phase transition between the $z$-FM and the $x$-PARA phase as $\Omega^{\prime}$ increases. The transition point at $\Omega^{\prime}=0.76$ can be determined by the sudden drop of $M_{z}$ to [math]. In Fig. 7(b), with increasing $\Omega^{\prime}$, there is also a first-order phase transition at $\Omega^{\prime}=0.68$, separating the $xy$-SP and $x$-PARA phase. In Fig. 7(c), as $\Omega^{\prime}$ increases, a continuous phase transition firstly takes place at $\Omega^{\prime}=0.92$, dividing the $xy$-SP and $y$-AFM phase, and then, a first-order phase transition happens at $\Omega^{\prime}=1.36$ as the sudden jump of $N_{y}$, dividing the $y$-AFM and $x$-PARA phase.

V summary and discussions

In summary, by employing the variational MPS method, we have studied quantum magnetic phase transitions in the deep Mott insulator of spin-orbit coupled two-level Bose atoms in 1D optical lattice. We find that the introduction of SOC brings a DM interaction and a transverse field in the effective model, which significantly modify the magnetic structures. When the transverse field in absent, the GS is in a gapped $z$-FM phase with ferromagnetic long range correlation in the $\hat{z}$ direction if the ratio of asymmetric interaction strength $\lambda>1$, and a gapless LL phase with algebraic decaying correlations if $\lambda<1$. The gapless LL phase is classified into ferromagnet or antiferromagnet when the DM interaction is zero and spiral phase when the DM interaction in nonzero. When the transverse field is introduced, the three kinds of gapless LL phase is broken to gapped $x$-PARA, $y$-AFM and gapless $xy$-SP phase, with finite long range correlation. From an analysis of the order parameters, we present rich phase diagrams on the parameter spaces. We note that the $xy$-SP phase distinguishes itself from the spiral phase in the nearest-neighbor spin-flip model induced by SOC, which is a LL phase, by the correlation function. Thus, our study provides a complete understanding of the SOC effect on the strongly interacting atoms in optical lattice system.

Interesting extensions of our present study include the magnetic phase transition of the strongly interacting artificial bosonic three-leg ladder, whose single-particle Hamiltonian is realized in Stuhl et al. (2015). The effective model will be a spin-1 chain, the GS of which would be accessible by the presented MPS algorithm. Furthermore, in the fermionic synthetic ladder, when take the nuclear spins into consideration, the interorbital spin-exchange interaction Cappellini et al. (2014); Scazza et al. (2014); Zhang et al. (2014) would couple individual ladders. Thus, the model Eq. (1) can be generalized to multiple two-leg ladders pierced by artificial magnetic flux. The impact of this inter-orbit coupling on the superfluid states has been thoroughly studied in Zhou et al. (2017). It induces exotic vortex states on the nuclear-lattice plane, which competes with the existing phases in the decoupled ladder and leads a rich phase diagram Zhou et al. (2017). How this inter-orbit coupling influences the quantum magnetism remains an open question for future research.

Acknowledgements.

This work was supported by the National Natural Science Foundation of China (Grants No. 11874434 and No. 11574405). Y.K. was partially supported by the International Postdoctoral Exchange Fellowship Program (Grant No. 20180052).