Quantum incommensurate Skyrmion crystals and Commensurate to In-commensurate transitions in cold atoms and materials with spin orbit couplings in a Zeeman field

In this work, we study strongly interacting spinor atoms in a lattice subject to a 2 dimensional (2d) anisotropic Rashba type of spin orbital coupling (SOC) and an Zeeman field. We find the interplay between the Zeeman field and the SOC provides a new platform to host rich and novel classes of quantum commensurate and in-commensurate phases, excitations and phase transitions. These commensurate phases include two collinear states at low and high Zeeman field, two co-planar canted states at Mirror reflected SOC parameters respectively. Most importantly, there are non-coplanar incommensurate Skyrmion (IC-SkX) crystal phases surrounded by the 4 commensurate phases. New excitation spectra above all the 5 phases, especially on the IC-SKX phase are computed. Three different classes of quantum commensurate to in-commensurate transitions from the IC-SKX to its 4 neighboring commensurate phases are identified. Finite temperature behaviors and transitions are discussed. The critical temperatures of all the phases can be raised above that reachable by current cold atom cooling techniques simply by tuning the number of atoms $ N $ per site. In view of recent impressive experimental advances in generating 2d SOC for cold atoms in optical lattices, these new many-body phenomena can be explored in the current and near future cold atom experiments. Applications to various materials such as MnSi, Fe$_{0.5}$Co$_{0.5}$Si, especially the complex incommensurate magnetic ordering in Li$_2$IrO$_3$ are given.

Introduction.-Quantum phases and quantum phase transitions in quantum spin systems have been an important and vigorous research field in material science for many decades [1][2][3]. It was widely believed and also partially established that due to the tremendous tunability of all the parameters, ultracold neutral fermions or spinor bosonic atoms loaded on optical lattices can provide unprecedented experimental quantum simulations and manipulations of these quantum spin systems [4,5]. However, so far, there still few theoretical works studying the effects of various kinds of spin-orbit couplings (SOC) on these quantum spin systems [6]. Recently the investigation and control of SOC have become the subjects of intensive research in both condensed matter and cold atom systems after the discovery of the topological insulators [7,8]. In the condensed matter side, there are increasing number of new quantum materials with significant SOC, including several new 5d transition metal oxides and heterostructures of transition metal systems [9]. In the cold atom side, there are very recent remarkable advances to generate SOC in optical lattices [5,[10][11][12][13][14][15][16][17][18]. It becomes topical and important to investigate possible dramatic effects of SOC in various quantum spin systems from both sides [9,[19][20][21][22][23].
In a recent work [22], we studied interacting spinor bosons at integer fillings loaded in a square optical lattice in the presence of non-Abelian gauge fields. In the strong coupling limit, it leads to a Rotated Ferromagnetic Heisenberg model (RFHM) which is a new class of quantum spin models with a strong SOC. We introduced Wilson loops to characterize frustrations and gauge equivalent class. For a special equivalent class, we identified a new spin-orbital entangled commensurate ground state: the Z − x state at the h = 0 axis in Fig.1. It supports not only commensurate magnons, but also a new gapped elementary excitation: in-commensurate magnons with two gap minima continuously tuned by the SOC strength. At low temperatures, these magnons lead to dramatic effects in all the thermodynamic quantities and spin correlation functions. The existence of the incommensurate magnons above a commensurate phase is a salient feature of the RFHM. They indicate the short-ranged in-commensurate order embedded in a long-range ordered commensurate ground state. Under the changes of some parameters such as the gauge parameters (α, β) at generic equivalent classes, or spin-anisotropy interaction, or external Zeeman field, they are the seeds to drive various new class of quantum phase transitions. In this paper, we study the RFHM in an external Zeeman field. We find that the Zeeman field drives the Z − x state at h = 0 to several new phases through new universality class of quantum phase transitions shown in Fig.1. The Z − x state at h = 0 remains stable only upto a low critical field h c1 . A ferromagnetic phase exists above a upper critical field h c2 , two canted co-planar states stand near β = 0 and β = π/2 respectively. Most importantly, a non-coplanar incommensurate Skyrmion ( IC-SkX ) crystal phases are surrounded by the 4 phases. The IC-SkX state reduces to a commensurate SkX state with a 2×4 spin-orbital structure at β = π/4. We determine the orbital and spin structures and the novel excitation spectra of all the collinear, coplanar and non-coplanar quantum phases. The QPT from the Z − x state to the 2 × 4 at β = π/4 is in the same universality class of zero density Mott to superfluid transition with the exponents z = 2, ν = 1/2, η = 0, that from the canted state to the IC-SkX is a new type of bosonic Lifshitz transition due to SOC [24] with the anisotropic dynamic exponents (z x = 1, z y = 3). We also outline the physical pictures of the FM to the canted state with z = 1 and the FM to the IC-SkX transitions with z = 2. New finite temperature transitions above all the phase are analyzed. Experimental implementations in both original and U (1) basis and various experimental detections of these new phenomena are presented. Implications to condensed matter systems with or without SOC are given. Our results show that quantum spin systems with SOC hosts new classes of quantum phases, excitations and quantum phase transitions. Some technical details are left in the Supplementary Materials (SM).
The RFHM with the two rotating matrices R(x, 2α), R(ẑ, 2β) where (α = π/2, β) in a Zeeman field Ω is [25]: where Ω is the Raman laser induced Zeeman field [10][11][12][13][14][15][16][17][18] along theẑ direction [31]. As shown in [22], the RFH at h = 0 has the translational symmetry, the time reversal T , the three spinorbital coupled Z 2 symmetries P x , P y , P z and most importantly, a hidden U (1) symmetry generated by U = e iφ i (−1) i S z i . The h breaks the T , P x , P y symmetries, but still keeps the translation, P z and the hidden U (1) symmetry. It can be shown that under the local rotation S i = R(x, π)R(ŷ, πn 2 ) S i , followed by the T , (β, h) → (π/2 − β, h). In the following, we use the dimensionless Zeeman field h = Ω 2SJ . Z − x phase and QPT at the low critical field h c1 .-It was shown in [22] that at h = 0, the Z − x state is the exact ground state with an excitation gap. It remains the exact ground state at a small h until the gap closes at h c1 . It only breaks the translational symmetry to two sites per unit cell, but keeps the P z and the U (1) symmetry.
Following [22], we perform SWE on Eqn.1 to leading order in 1/S. Introducing the two Holstein-Primakoff (HP) bosons for the two A/B sublattices respectively and introducing a unitary transformation, one can put the Hamiltonian in the diagonal form ( See SM): where E 0 = −2N JS 2 is the ground state energy ( the same as that at h = 0 ), k belongs to the reduced Brillioun zone (RBZ) and the excitation spectrum ω ± (k) = 1 − 1 2 cos 2β cos k y ± 1 2 cos 2 k x + (sin 2β sin k y − h) 2 .  Fig.2b. The excitation spectra, the quantum phase transitions between ( among ) these phases with the dynamic exponents z = 1, z = 2 and the anisotropic one (zx = 1, zy = 3) ( abbreviated by z = 3 in the figure to save space ) and finite temperature phase transitions above all the phases are discussed in the text.
As shown in [22], at h = 0 ( horizontal axis in Fig.1 ), the gap minimum of the lower branch ω − (k) is at (0, 0) when 0 < β < β 1 ( called C − C 0 magnons ), at two minima (0, ±k 0 y ) when β 1 < β < β 1 ( called C − IC magnons ), at (0, π) when β 2 < β < π ( called C − C π magnons ). However, any H > 0 will turn all these magnons into the C − IC magnons located only at one minimum k 0 = (0, 0 < k 0 y < π) where k 0 y = k 0 y (β, H) depends on both β and H, its constant contour was shown in Fig.1a in SM. Expanding ω − (k) near the minimum k = k 0 + q leads to: where the ∆ Z (β, h) is the gap and m Z,x (β, h), m Z,y (β, h) are their two effective masses. The lower critical Zeeman field h c1 is determined by ∆ Z (β, h c1 ) = 0 where the dynamic exponent z = 2. Its expression is given in the SM and shown in Fig.1 The two effective masses remain non critical (SM). The condensation of the C −IC magnons indicates a transition from the Z − x state into a IC-SkX phase. The nature of this transition and the IC-SkX phase will be explored from h c1 < h < h c2 below.
FM phase and QPT at the upper critical field h c2 .-At a strong Zeeman field H ≫ J > 0 , the spin is a FM state ( subject to quantum fluctuations ) as shown in Fig.1. It breaks no symmetry of the Hamiltonian. One only need to introduce one HP boson. After performing a Bogliubov transformation, we obtain: where ω k = (h − cos 2β cos k y ) 2 − cos 2 k x − sin 2β sin k y .
We find that the ω k always has two degenerate minima located at (0, k 0 y ) and (π, k 0 y ) where 0 ≤ k 0 y ≤ π. The upper critical field h c2 , determined by the vanishing gap at the two minima, takes a piece-wise form: where the two different pieces indicate transitions to two different states. At β = β 1 , β 2 , the two expressions coincide [32].
When 0 < β < β 1 , expanding ω k around the k 0 y = 0 leads to: where cos 2β, c F = sin 2β where z = 1. At β = 0, c F = 0, Eqn.6 recovers the relativistic form. The simultaneous condensations of the C magnons at k 0 = (0, 0) and (π, 0) indicates a transition from the FM state into a canted phase shown in Fig.1. The nature of the canted phase will be explored from h L < h < h c2 below. Increasing β along h = h c2 where ∆ F 0 = 0, the slope v F,y (β) − c F (β) of the dispersion Eqn.6 along q y > 0 decreases. At β = β 1 , the slope vanishes. As shown below the dynamic exponent along q y directions becomes anisotropic (z x = 1, z y = 3). This is a multi-critical point with (z x = 1, z y = 3) of the three phases: FM ( collinear ), canted ( Co-planar ) and IC-SkX ( Non-Coplanar ) phase.
When β 1 < β < β 2 , expanding ω k around the minima (0, 0 < k 0 y < π), we obtain a similar form as Eqn.S5: where , m F,y = sin 2 2β/m F,x , z = 2. At β = β 1 , h c2 = 1+cos 2β, we find m F,x = 0 and m F,y = ∞ which match the anisotropic (z x = 1, z y = 3) dynamic behaviors in Eqn.6. The condensation of the C − IC magnons indicates a transition from the FM state into a IC-SkX phase shown in Fig.1. The nature of this IC-SkX phase and the transition will be explored from h c1 < h < h c2 below.
Canted phase and QPTs at h c2 and the left critical field h L .-When β is near the Abelian case β = 0, we first determine the simplest classical state to be a FM state in the XZ plane ( Fig.1 ) whose tilted angle with the z axis is: where h c2 = 2 cos 2 β matches the upper critical field achieved from h > h c2 Eqn.S18. Note that Eqn.8 is reached by a classical condition, while Eqn.S18 is by a quantum condition to the order of 1/S. This match indicates h c2 does not receive quantum corrections at least to the order of 1/S. The most general form of the classical ground state can be obtained by applying the U (1) symmetry operator [33] on the simplest FM state in the XZ plane denoted as Q 1 = (0, 0): ψ(φ) = e iφ i (−1) x S z (0, 0) which leads to the classical canted state: It contains the two ordering wavevectors Q 1 = (0, 0) and Q 2 = (π, 0) which match those of the condensed C-magnons reached from h > h c2 . Setting φ = 0, π recovers the two FM states in Eqn.8 which does not break the translational symmetry. While setting φ = π/2 leads to a canted ( co-planar ) state in the Y Z plane. When φ = 0, π, Eqn.9 breaks the translational symmetry, the P z and the U (1) symmetry. It is the U (1) symmetry breaking which leads to a gapless Goldstone mode. However, as shown below, due to the rather special symmetries and symmetry breaking of the SOC systems, the Goldstone mode takes very peculiar form. It is most convenient to perform SWE on the simplest FM state in the XZ plane Eqn.8. We first make a global rotation R y (θ) to align the spin quantized axis along the Z axis then only need to introduce one HP boson to perform the SWE H = H 0 + H 1 + H 2 + · · · . We obtain H 0 = −2N JS 2 cos 2 β + h 2 4 cos 2 β , H 1 = 0 which is dictated by the correct classical ground state Eqn.8 and We find the Goldstone mode φ is located at k 0 = (π, 0) and takes the rather peculiar form: where v g,x , v g,y , c g = h tan β are listed in the SM. At β = 0, c g = 0, Eqn.S23 recovers the relativistic form.   Fixing 0 < β < β 1 , as one increases to h c2 (β) from below, a roton minimum develops at (0, 0) as shown in Fig.2a. Its spectrum takes a similar form as Eqn. 6: At h = h c2 = 1 + cos 2β, 0 < β < β 1 , we find the Goldstone mode at (π, 0) in Eqn.S23 and the roton mode at (0, 0) in Eqn.S25 achieved from below [34] h ≤ h c2 coincide with Eqn.6 achieved from above h ≥ h c2 , namely: √ cos 2β and c g = c r = c F = sin 2β. This indicates the transition from the FM to the canted state maybe a second order transition with z = 1.
(b) Bosonic Lifshitz type of transition from the canted state to the IC-SkX phase at the left critical field h L .
At fixed h < h c2 (β 1 ), as the SOC strength increases, the slope of the Goldstone mode v g,y − c g in Eqn.S23 along q y > 0 decreases. Setting the slope vanishing, we obtain h L : which is shown in Fig.1. Combining Eq.(S18) and Eq.(13) and setting h c2 = h L , we obtain β = β 1 and (Golden ratio). As shown in the SM, near h L , by expanding the Goldstone mode Eqn.S23 to higher orders, we find ω(q x = 0, q y ) = (v g,y − c g )q y + c 3 q 3 y + · · · where c 3 > 0 is given in the SM. When v g,y − c g > 0, the minimum position is at q 0 y = 0, so it is in the canted state. When v g,y − c g < 0, the minimum position is at q 0 y = ( cg −vg,y c3 ) 1/2 , it is in the IC-SkX state with the orbital order at (π, q 0 y ). Indeed, this infinitesimal small orbital order connects the one at h c2 , β = β + 1 smoothly to the one at h c1 , β = 0 + . This is a bosonic type of Lifshitz transition [24], however, with the odd power of terms such as q y , q 3 y , .... which is due to the SOC. So it is completely new class of bosonic type of Lifshitz transition with the anisotropic dynamic exponent (z x = 1, z y = 3).
One can also show that ∂hL ∂β β=0 = ∂hc1 ∂β β=0 = 2. However, h L is always above h c1 , so there is always a narrow window of IC-SkX phase sandwiched between the collinear Z − x phase and the co-planar canted phase. There is NO direct transition between the two. This is consistent with the contour (0, k y 0 → 0) in the β → 0 limit from h < h c1 in Fig.1.
2 × 4 Skyrmion crystal phase at β = π/4 and QPTs at h c1 and h c2 .-As shown in [22], due to the explicit U (1) symmetry in the U (1) basis ( Eqn.18 ), the anomalous spin correlation functions vanish. So the most general form of classical state at β = π/4 can be most conveniently determined in the U (1) basis: the condensation at (0, π/2) of the C − IC magnons at β = π/4 directly implies the resulting state takes the form Transforming back to the original basis, then acting on it by the U (1) symmetry operators leads to: which breaks the translational, P z and the U (1) symmetry. Without losing any generality, one can pick up φ = 0 which breaks the U (1) symmetry, leading to a gapless mode whose form will be determined in Eqn.16. Using Eqn.S29, we find [35] a oscillating Skyrmion den- where i, j, k are taken as three lattice points around a square. At φ = 0, it reduces to We have also taken a general 4 × 4 structure, then numerically minimize the classic ground state energy with respect to the 32 parameters (16 polar angles + 16 azimuthal angles). Numerical results always find the configuration described by Eqn.S29. Minimization of the classical energy leads to the two independent polar angles θ A , θ B in the two sublattices shown in Fig.2b.
From Eqn.S29, we make suitably chosen to align the spin quantization axis along the Z axis. Then one need only introduce two HP bosons a/b for the two sublattices A/B respectively and perform a Bogoliubov transformation to obtain: where ω ± (k) = 1 The expressions of C k , D k are given in the SM.
As shown in the SM, D 2 k k=0 = 0, therefore ω − (k = 0) = 0 which is the gapless Goldstone mode φ at k = 0. Performing a long wavelength expansion around the Γ point leading to the Goldstone mode in the long wavelength limit: where v G,x , v G,y are given in the SM.
In fact, as shown in Fig.2b, putting θ A = 0, θ B = π and θ A = θ B = 0, one can also push the calculations to the Z −x state at h < h c1 and the FM state at h > h c2 , but in different gauges than those used in the previous sections. As expected, the minimum positions of excitations shift at different gauges. The critical behaviors are shown in Fig.3. As shown in the Fig.3b, as h → h − c2 , there is also roton mode developing at (0, π) which takes the relativistic form: where k = q + (0, π) and ∆ R = 4 Pushing the expansion to k 4 in both Eqn. 16 and Eqn.17 ( see SM), we find the effective masses of both the Goldstone mode and the roton mode coincide with the m F,x , m F,y Eqn.7 achieved from the FM state, so z = 2. Similarly, at Pushing the expansion to k 4 in Eqn.16( see SM), we find the effective masses of the Goldstone mode coincide with the m Z,x , m Z,y achieved from the Z − x state Eqn.S5, so z = 2.
For any state |ψ L on the left β ≤ π/4, one can get the state on the right |ψ R = T e i π 2 σx e i π 2 yσz |ψ L . Indeed, applying the operation on the state Eqn.9 on the left side leads to the state on the right hand side S z = cos θ, S + = [S − ] † = − sin θ(−1) y e −i(−1) x φ which contains two ordering wavevectors Q 1 = (0, π) and Q 2 = (π, π). Setting φ = π gives the state shown on the right in Fig.1. Applying the operation on the state Eqn.S29 leads back to itself as expected. The right critical field h R is given by Eqn.13 by setting β → π/2 − β which is a reflected image with respect to β = π/4 in Fig.1. Therefor, one only need to focus on the left half of Fig.1 with 0 < β < π/4.
Quantum phase transitions at T = 0 .-When β = π/4, h ≤ h c1 , from Eqn.S5 with the critical behavior ∆ Z ∼ (h c1 −h) 1 and when β = π/4, h ≥ h c1 , from Eqn. 16 with the critical behavior v x (h c1 ) = 0 and v y (h c1 ) = 0 shown in Fig.3a, we conclude the Z − x to the 2 × 4-SkX phase at β = π/4 driven by the condensation of the magnons at (0, π/2) is in the same universality class of zero density Mott to superfluid transition with z = 2, ν = 1/2, η = 0 with possible logarithmic corrections at its upper critical dimension d u = 2 [37]. In addition to the U (1) symmetry breaking, there are also translational symmetry breaking along the y direction in the 2×4 SkX phase, The connection between the spin-orbital orders of both phases and the complex order parameter ψ ∼ e iφ in the effective long wavelength action is precisely given by Eqn.S29. Indeed, when h < h c1 , ψ = 0, the system is in the Z − x state. When h > h c1 , ψ = 0, the system is in the 2 × 4-SkX state.
In the transition from the 2 × 4-SkX phase to the FM phase from h < h c2 shown in Fig.3b, in addition to the Goldstone mode at (0, 0), there is also a roton mode at (0, π) with the spectrum Eqn.17. The two critical modes with z = 2 are decoupled in the linear SWE, but will be coupled in higher orders. In principle, one need to go beyond any orders of the SWE to understand the quantum critical behaviors. When approaching the transition from above h > h c2 , there are two degenerate minima at (0, 0) and (0, π). Following the methods in [43][44][45][46], one can construct a quantum Ginzburg-Landau action in terms of the two order parameters at (0, 0) and (0, π) respecting the T ran, P z , U (1) symmetries in the FM state to study this new class of QPT.
If β = π/4, following a fixed β route, the minimum positions of the dispersions in Eqn.S5 will continuously change from both h < h c1 , 0 < β < π/2 or h > h c2 , β 1 < β < β 2 , so it is more convenient to fol-low the constant contour k y 0 = k y 0 (β, H) in Fig.1a in SM to study the universality class of the transition. In fact, we expect to connect any 0 < k y 0 < π contour at β 1 < β < β 2 , h = h c2 smoothly to 0 < β < π/2, h = h c1 . It stands for the IC-SkX with orbital orders at (0, k y 0 ) and (π, k y 0 ). The form of the Goldstone mode Eqn.16 will also change. The transition from the Z − x to the IC-SkX state at every point on the h c1 line holds a different class of transition. Similarly, the transition from the FM to the IC-SkX state at every point on the line h = h c2 , β 1 < β < β 2 holds a different class of transition. The line β = π/4 enjoys the special symmetry.
In the transition from the canted phase to the FM phase from below h L < h < h c2 shown in Fig.2a, in additional to the Goldstone mode at (π, 0), there is also the roton mode Eqn.S25 developing at (0, 0) which touches zero at h c2 . The two critical modes with z = 1 are decoupled in the linear SWE, but will be coupled in higher orders. In principle, one need to go beyond any orders of the SWE to understand the quantum critical behaviors. When approaching the transition from above h > h c2 , 0 < β < β 1 , the spectrum Eqn.6 has two degenerate minima at (0, k y 0 ) and (π, k y 0 ) where the k y 0 continuously decreases until hitting zero at h c2 . As said above, one can study the transition from the FM to the IC-SkX from above h > h c2 , β 1 < β < β 2 by following a constant contour k y 0 = k y 0 (β, H). Here there is no such contour to follow. At β = 0, it corresponds to a FHM in a staggered (π, 0) Zeeman field in the U (1) basis ( Eqn.18 ), c = 0, k y 0 = 0 is kept along the contour. We expect the transition from the FM to the canted state along the whole line 0 < β < β 1 , h = h c2 is different from that at β = 0.
As shown in the previous section, the transition from the canted state to the IC-SkX state at h = h L is a bosonic Lifshitz type of transition. The orbital order of the IC-SkX near the boundary are at (π, k y 0 ) and (0, k y 0 ) with k y 0 → 0 + . The corresponding spin order is given by Eqn.S29 if replacing iπ/2y by ik y 0 y. Indeed, the h = h + L connects the k y 0 = ( cg −vg,y c3 At the two multi-critical points β = β 1 , β 2 and h = h c2 where the three phases: FM ( collinear ), canted ( Co-planar ) and IC-SkX ( Non-Coplanar ) phase meet, the gap ∆ and the slope v y − c = 0 vanish at the same time. The nature of this novel multi-critical points with (z x = 1, z y = 3) remain to be determined.
Finite Temperature phase transitions-As argued in [22], there is only one finite temperature phase transition above the Z −x phase. The canted phase and the IC-SkX phase also breaks the continuous U (1) symmetry, so any finite T > 0 will transfer it to an algebraic order.
In the canted phase, from Eqn.9, one can see that at any T > 0, the Goldstone mode fluctuations Eqn.S23 lead to S + = 0, so the transverse spin correlation functions also become algebraic order at the two ordering wavevectors Q 1 = (0, 0) and Q 2 = (π, 0). So there is only one finite temperature phase transition T * KT above the canted phase to destroy the algebraic order. In view of the rather peculiar anisotropic form of the Goldstone mode [47] in Eqn.S23, it remains interesting to study if the T * KT falls into the conventional Kosterlize-Thouless (KT) one.
In the 2 × 4 SkX phase, from Eqn.S29, one can see that at any T > 0, the Goldstone mode fluctuations Eqn.16 also lead to S + = 0, so the transverse spin correlation functions also become algebraic orders at the two ordering wavevectors Q 1 = (0, π/2) and Q 2 = (π, π/2). So there are two finite temperature phase transitions above the 2 × 4 SkX state: one KT transition T KT in the transverse spin sector to destroy the algebraic order and then another Ising ( Z 2 ) transition in the longitudinal spin sector T 2 to destroy the A and B sublattice symmetry breaking. We expect T KT < T 2 . At all the quantum phase transition boundaries in Fig.1 The quantum phases, excitation spectra and finite temperature transitions at a general fractional value k y 0 = p/q ( p and q are prime numbers ) or any irrational value k y 0 along a constant contour k y 0 = k y 0 (β, H) in Fig.1 remains an outstanding problem and is under the current investigations [58].
Following [22], when T ≫ J, h c2 , we may also perform a high temperature expansion in terms of J/T and h/T where mixing terms in J/T and h/T are expected. The connections between the Wilson loop and specific heat can also be established.
Experimental realizations and detections-In the original basis Eqn.1, the gauge field configuration is achieved by putting π/2σ x in the x− bond, βσ z in the y− bond. Then the Raman laser induced Zeeman field Ω is alongẑ direction. The β = 0 case has been realized in previous experiments [5,[10][11][12][13][14][15][16][17][18]. As pointed out in [16], the βσ z in the y− bond can be achieved by adding spin-flip Raman lasers or by driving the spin-flip transition with RF or microwave fields.
Following [22], one can work out the thermodynamic quantities such as magnetization, uniform and staggered susceptibilities, specific heat and Wilson ratio at the low temperatures in all the 5 phases in Fig.1. Of course, due to the Goldstone modes in the canted and the SkX phases, the specific heat takes the power law C v ∼ T 2 . Similarly, one can work out various kinds of spin correlation functions at the low and high. temperatures in all the 5 phases in the Fig.1.
As shown in [22], the U (1) basis may be more experimentally easily realized. In the U (1) basis, the RFHM in the Zeeman field Eqn.1 becomes where H U(1) is the RFHM in the U (1) basis at Ω = 0 given in [22] and the Zeeman field becomes a staggered one ( Fig.4). Then applying the R(x, n 1 π) on all the states shown in Fig.1 leads to the corresponding states in the U (1) basis. The thermodynamic quantities are gauge invariant, so are the same in both basis. But the spin-correlation functions are gauge dependent, need to be re-evaluated in the U (1) basis at both low and high temperatures. As shown in [22], all these physical quantities can be detected by atom or light Bragg spectroscopies [38,39], specific heat measurements [40,41] and the In-Situ measurements [42].
To detect the canted phase, at T = 0, the transverse Bragg spectroscopy will display sharp peaks at Q 1 = (0, 0) and Q 2 = (π, 0). However at 0 < T < T * KT , the transverse peaks at Q 1 and Q 2 will be replaced by power law singularities. At T > T * KT , the power law singularities disappear.
To detect the 2 × 4 SkX state, at T = 0, the elastic longitudinal Bragg spectroscopy will display a sharp peak at (π, 0) the transverse Bragg spectroscopy will display sharp peaks at Q 1 = (0, π/2) and Q 2 = (π, π/2). However at 0 < T < T KT , the transverse peaks at Q 1 and Q 2 will be replaced by power law singularities, the longitudinal peak remains sharp. At T KT < T < T 2 , the power law singularities disappear, but the longitudinal peak remains sharp. When T > T 2 , the longitudinal peak disappears. Similarly, a general IC-SkX phase can be detected by the incommensurate peaks at (0, k 0 y ) and (π, k 0 y ) in the transverse Bragg spectroscopy. The roton minimum dropping shown in Fig.2a and Fig.3b can also be monitored through in-elastic Bragg spectroscopies.
Conclusions and Perspectives.-Our main results on new class of quantum phases and phase transitions of RFH in a Zeeman field are presented in Fig.1. From down to up driven by the Zeeman field h, the IC-SkX intermediate phase sandwiches between the Z−x state at h < h c1 and the FM state at h > h c2 . From left to right driven by the SOC strength β, the IC-SkX intermediate phase sandwiches between the canted state at h L < h < h c2 to another canted state at h R < h < h c2 . At β = π/4, the IC-SkX reduces to a 2 × 4 commensurate SkX. It displays spin-orbital correlated collinear, co-planar ( canted), non-coplanar ( Skyrmion crystal ) phases which are due to SOC in a bipartite lattice. They can be contrasted with collinear, spiral, non-coplanar phases found in geometrically frustrated lattices [1]. It was shown in [22] that β = π/4 falls into the most frustrated regime where the Wilson loop W R = −1, the Dzyaloshinskii-Moriya (DM) term dominates. So the 2 × 4 commensurate SkX may be realized in some material with a strong DM interaction. Indeed, a 2d skyrmion lattice has been observed in between h c1 = 50mT and h c2 = 70mT in some chiral magnets [50] M nSi or a thin film of F e 0.5 Co 0.5 Si [50].
It is interesting to compare the roton mode inside the canted phase Fig.2a and the 2 × 4 SkX phase Fig.3b with the roton minima in continuous superfluid systems and also lattice systems without SOC. In the superfluid helium system, the roton inside the superfluid state indicates the short -ranged solid order inside the off diagonal long-ranged SF order [51,52]. As the pressure increases, the roton minimum drops and signals a first order transition to a solid order ( or a putative supersolid order ). In 3d, the roton sphere is a 2d continuous manifold, so its dropping before touching zero signals a first order transition. Similarly, in a 2d electron-hole semi-conductor bilayer (EHBL) system [53] or 2d bilayer quantum Hall (BLQH) systems [54], the roton circle inside the exciton superfluid is a 1d continuous manifold, so its dropping before touching zero signals a first order transition. In boson Hubbard models in various lattices [43-46, 52, 55, 56], the roton dropping indicates a transition from the SF to a supersolid. The most astonishing difference is that in the all the previous ( exciton ) or lattice superfluid systems, the roton drives to a lattice symmetry breaking state, here the roton plays just the opposite role: restoring the lattice symmetry.
It is stimulating to compare the C-IC transition from the canted state to the IC-SkX at h = h L to the C-IC transition described by a 2 = 1 dimensional Pokrovsky-Talapov (PT) model in the Bilayer quantum Hall systems at total filling factor ν T = 1 [57]. The SOC parameter β here plays a similar role as the in-plane Zeeman field B || in the BLQH. At the fixed h < h c2 (β 1 ) = √ 5+1 2 , setting h = h L in Eqn.13 leads to a critical β c (h). Then as the β increases from 0 to β c , the tilted angle of the canted state θ in Eqn.8 increases from cos θ 0 = h/2 accordingly until reaching its maximum at cos θ c = h L /h c2 < 1 evaluated at β c . Then after β > β c , the system can not follow the SOC anymore, just give up the follow and fall into a IC-SkX state with the phase gradient along the y axis ∂ y φ = k y 0 = ( cg −vg,y c3 ) 1/2 . There are KT transition on both sides of the C-IC transition. However, in the BLQH, there is only KT transition on the IC side through dislocations in the domain wall structure of the IC phase.
It is instructive to compare the QPT at h = h L , h R with the quantum dimer model [28] which describes the transition from one valence bond solid to another VBS state. It was known that near the solvable RK point [59], in the height representation, the transition can be described by the low energy effective quantum bosonic Lifshitz action L RK = 1 2 [(∂ τ h) 2 +K 2 (∇h) 2 +K 4 (∇h) 4 +· · · ] where K 2 = 0 at the RK point with z = 2 [28]. However, the canted state to the IC-SkX transition in our RFH+H is described by an quantum bosonic Lifshitz type action with the anisotropic dynamic exponents (z x = 1, z y = 3) and the appearance of linear q y and cubic q 3 y . It is the SOC which breaks the rotational and inversion symmetry leading to the appearance of these odd terms. The exact action form and the nature of the transition and its scaling functions at finite temperatures will be presented in a future publication. It is possible gravity dual [60] may also be interesting to work out. Of course, the underlying physical quantity in the quantum dimer model is the dimer density on a bond ( spin singlet bond ), while here it is the quantum spin. RFHM in generic equivalent class [25] in a Zeeman field and RAFH in a Zeeman field will also be presented in separate publications. In short, Quantum spin systems with SOC display rich and novel class of quantum phases and phase transitions. The results achieved in this paper just reveals a tip of an iceberg.
We thank I. Bloch for discussions on current experimental achievements. In this Supplementary material, we provide some details to the results achieved in the main text. We also split it into four sections to discuss the Z − x state below h c1 , the FM state above h c2 , the canted state on the left side h L < h < h c2 and the 2 × 4 SkX state at β = π/4. The canted state on the right side h R < h < h c2 is reflected to that in the left by the reflection symmetry (β, h) → (π/2 − β, h).
The unitary transformation leading to the excitation spectrum in the Z − xstate: ω ± (k) = 1 − 1 2 cos 2β cos k y ± 1 2 cos 2 k x + (sin 2β sin k y − h) 2 (S1) is: where the matrix elements are given by: Minima location (0, k 0 y ) of ω ± (k) is one of the roots of following quartic equation sin 2 2β sin 4 k y − 2h sin 2β sin 3 k y + (1 + h 2 − sin 2 2β − sin 4 2β) sin 2 k y + 2h sin 3 2β sin k y − h 2 sin 2 2β = 0 (S4) It turns out that there is always one and only one physical root. The constant contour of k 0 y is shown in Fig.S1a. Expanding ω − (k) near the minimum k = k 0 + q leads to: where the ∆ Z (β, h) is the gap of the C − IC magnons and m Z,x (β, h), m Z,y = m Z,y (β, h) are their two effective masses. The lower critical magnetic field h c1 is determined by ∆ Z (β, h c1 ) = 0 and is given by: When β = π/4, it simplifies to: The coefficient is nonzero for β = 0, π/2, thus we obtain ∆ ∼ (h c1 − h) 1 whose slope is given in Fig.S1b. The values of the two effective masses at h c1 are shown in Fig.S1c.

II. DETERMINATIONS OF THREE SEGMENTS OF hc2 FROM HIGH FIELD SWE
The crucial difference of the upper critical field h c2 from the lower critical field h c1 is that one has to split h c2 into 3 different segments ( namely, piece-wise ). It indicates transitions to 3 different class of states: two canted states and one IC-SkX states shown in Fig.1.
The spin wave dispersion in the FM state above h c2 is Due to the reflection symmetry about β = π/4, we only need to focus on 0 < β < π/4. It is easy to see that there are always two degenerate minima located at k x = 0, π. The minimization in k y leads to: 0 = ∂ω k ∂k y = cos 2β sin k y (h − cos 2β cos k y ) (h − cos 2β cos k y ) 2 − 1 − sin 2β cos k y (S11) If Eq.(S11) has a real solution k 0 y , then plugging it back into Eq.(S10) leads to: ω min = sin k 0 y (h cos 2β − cos k 0 y ) sin 2β cos k 0 y (S12) The gap vanishing condition is: In fact, Eq.(S11) is a quartic equation of cos k y sin 2 2β cos 2 k y + (cos 2 2β − cos 2 k y )(h − cos 2β cos k y ) 2 = 0 (S14) If there exists one root with cos k 0 y ≤ 1 in Eq.(S14), substituting Eq.(S13) into Eq.(S14) leads to If all positive roots of Eq.(S14) require cos k 0 y > 1, then the minimum is located at cos k 0 y = 1, namely, k 0 y = 0. Then substituting cos k 0 y = 1 into Eq.(S10) leads to: which is also the condition ensuring a real spectrum.
Combining the two piece-wise h c2 equations leads to: which is shown in Fig.1. After extending to β ∈ (0, π/2), we obtain: The two different piece-wise forms of h c2 in the regime I and II indicates transitions to two different states: canted state and IC-SkX state respectively with the dynamic exponents z = 1 and z = 2 respectively.

III. GOLDSTONE AND ROTON MODE IN THE CANTED STATE
The spin wave excitation spectrum in the canted phase is cos k x + sin 2 β 1 − h 2 4 cos 4 β cos k y , The excitation spectrum is always gapless at k 0 = (π, 0).
Expanding around k = k 0 + q, we obtain the Goldstone mode Eqn.M10 As h increases to h c2 , 0 < β < β 1 , there is also a roton minimum developing at (0, 0) shown in Fig.S1a: At h = h c2 = 1 + cos 2β, we find the Goldstone mode at (π, 0) and the roton mode at (0, 0) achieved from below h ≤ h c2 , 0 < β < β 1 coincide with those achieved from above h ≥ h c2 , 0 < β < β 1 , namely: v g,x = v r,x = v F,x = 1, v g,y = v r,y = v F,y = √ cos 2β and c g = c r = c F = sin 2β. This indicates the transition from the FM to the canted state maybe a second order transition with z = 1.
Now we look at the Lifshitz type of transition at h = h L from the canted state to the IC-SkX state. We need to perform higher-order gradient expansion around (π, 0) in the canted state near h = h L to see the nature of the transition: where Note that v xx is not positive-define, but it will not cause any instabilities due to v x > 0. Both v y and v yy are positive. Without losing the physics of the transition from the canted state to the IC-SkX, setting q x = 0 simplifies Eqn.S26 to: Because the transition first happens in q y > 0 direction, one can see that when v g,y − c g > 0, the minimum position is at q 0 y = 0, so it is in the canted state. When v g,y − c g < 0, the minimum position is at q 0 y = ( cg−vg,y c3 ) 1/2 where c 3 = vyy 2vg,y + c ′ , then it is in the IC-SkX state with the orbital order at (π, q 0 y ). Indeed, this infinitesimal small orbital order connects the one at h c2 , β = β + 1 smoothly to the one at h c1 , β = 0 + . This is a bosonic type of Lifshitz transition, however, with the odd power of terms such as q y , q 3 y , .... which is due to the SOC. So it is completely new class of bosonic type of Lifshitz transition with the anisotropic dynamic exponents z x = 1, z y = 3.

IV. GOLDSTONE AND ROTON MODE IN THE 2 × 4 SKYRMION CRYSTAL STATE
The classical state of the 2 × 4 SkX state in Eqn.M12 can be written as: S i = S(sin θ A cos(φ A − i y π/2), sin θ A sin(φ A − i y π/2), cos θ A ), i ∈ A; S j = S(sin θ B cos(φ B + j y π/2), sin θ B sin(φ B + j y π/2), cos θ B ), j ∈ B; (S29) Using Eqn.S29, we find the classic ground state energy is Because sin θ A sin θ B > 0, the minimization requires φ A + φ B = 0. For θ A and θ B , the minimization requires Solving the equation leads to the two polar angles shown in Fig.S2b. Performing the SWE on the classical state leads to the excitation spectrum in Eqn.M13 We now evaluate D 2 k at the Γ = (0, 0) point, Note that Eqn.S31 leads to (A 0 + 1)(B 0 − 1) − 1 = 0, therefore D 2 k k=0 = 0. It indicates ω − (k = 0) which is the gapless Goldstone mode at k = 0. Now we perform a long wavelength expansion around the Γ point, where we have introduced C 0 = C k k=0 = A 2 0 + B 2 0 − 2 cos(θ A + θ B ) + 2(A 0 cos 2 θ A − B 0 ) + 2 cos 2 θ A > 0 D x = (A 0 + cos 2θ A )(B 0 − 1) − cos 2 (θ A + θ B ) thus we can extract the Goldstone mode in the long wavelength limit: where its velocity v 2 G,x = Dx C0 , v 2 G,y = Dy C0 are shown in Fig.S2b. In fact, as shown in Fig.S2b, putting θ A = 0, θ B = π and θ A = θ B = 0, one can also push the calculations to the the Z − x state at h < h c1 and the FM state at h > h c2 , but in different gauges than those used in the previous sections. As expected, the minimum positions of excitations may shift at different gauges. The gaps along the whole line β = π/4 are shown in Fig.S2a.