From antiferromagnetic order to magnetic textures in the two dimensional Fermi Hubbard model with synthetic spin orbit interaction
FFrom antiferromagnetic order to magnetic textures in the two dimensional FermiHubbard model with synthetic spin orbit interaction
Jiří Minář
Centre for Quantum Technologies, National University of Singapore, Singapore
Benoît Grémaud
Centre for Quantum Technologies, National University of Singapore, SingaporeDepartment of Physics, National University of Singapore, Singapore andLaboratoire Kastler Brossel, Ecole Normale Supérieure CNRS, UPMC; 4 Place Jussieu, 75005 Paris, France
We study the interacting Fermi-Hubbard model in two spatial dimensions with synthetic gaugecoupling of the spin orbit Rashba type, at half-filling. Using real space mean field theory, wenumerically determine the phase as a function of the interaction strength for different values of thegauge field parameters. For a fixed value of the gauge field, we observe that when the strengthof the repulsive interaction is increased, the system enters into an antiferromagnetic phase, thenundergoes a first order phase transition to an non collinear magnetic phase. Depending on the gaugefield parameter, this phase further evolves to the one predicted from the effective Heisenberg modelobtained in the limit of large interaction strength. We explain the presence of the antiferromagneticphase at small interaction from the computation of the spin-spin susceptibility which displays adivergence at low temperatures for the antiferromagnetic ordering. We discuss, how the divergenceis related to the nature of the underlying Fermi surfaces. Finally, the fact that the first order phasetransitions for different gauge field parameters occur at unrelated critical interaction strengths arisesfrom a Hofstadter-like situation, i.e. for different magnetic phases, the mean-field Hamiltonians havedifferent translational symmetries.
PACS numbers: 67.85. d 05.30.Fk 37.10.Jk 71.70.Ej
I. INTRODUCTION
The recent progress of experiments using cold atomicgases , in particular in implementing artificial gaugefields , has open the door to the studies of a wholeclass of model Hamiltonians, some directly inherited fromcondensed matter physics (quantum Hall effects), but,more saliently, some are genuily generating new physi-cal situations, allowing physicists to further develop andtest theoretical ideas, like topological phases , non-abelian particles or mixed dimensional systems .In particular, the experiments involving spinors, made ofeither bosons or fermions in different Zeeman sub-levels,are now able to produce non-abelian gauge-fields, lead-ing to a kinetic term allowing for a modification of theinternal degrees of freedom along the propagation of theparticle . In two-dimensional lattices, the non-abeliangauge fields result in tight-binding Hubbard models withspin-flip hoping terms, i.e. the hoping matrices are notdiagonal anymore in the spin degrees of freedom. Amongall the possibilities, an artificial gauge-field mimicking aspin-orbit coupling term (see below) is probably themost well-known and studied situation, for two reasons:(i) it corresponds to a spatially independent vector po-tential leading to relatively simple analytical treatment,especially in the bulk situation; (ii) it leads to highly non-trivial features like broken time-reversal ground statesand/or magnetic textures with topological properties,like skyrmion crystal .The cases of two components bosons or fermions in thepresence of a spin-orbit coupling have been the objects of recent analytical and numerical studies . In par-ticular, in the case of repulsive interactions, they haveemphasized various magnetic ordering and textures de-picted by the effective spins. One must note that in thelarge interaction strength, and close to half-filling, bothbosonic and fermionic situations can be described by ef-fective and quite similar Heisenberg models involving thespin degrees of freedom only (see below). The case offermions with attractive interaction has also been stud-ied, emphasizing the impact of the spin-flip terms on thepairing states, like the BCS-BEC crossover or the Fulde-Ferrel-Larkin-Ovchinnikov (FFLO) phase . The in-stabilities of attracting bosons with such a spin-orbit cou-pling have also been considered .Even though the Bose-Hubbard model and the Fermi-Hubbard model with a spin-orbit coupling depicts simi-lar phases in the strong (repulsive) interaction limit, thebehavior for small coupling is obviously different: theMott regime of the bosonic models turns into a com-plex superfluid regime whereas fermionic models are ex-pected to be in a Mott insulator state. Still, in the latercase, the evolution of the magnetic ordering from thenon-interacting situation towards the Heisenberg modelregime remains largely unexplored. In the present pa-per, using both a linear response approach and real spacemean-field theory, we emphasize that the Fermi-Hubbardmodel at half-filling, in the presence of a spin-orbit cou-pling, still depicts, at low interaction, an antiferromag-netic phase, corresponding to the Fermi-Hubbard model without spin-orbit coupling. The impact of the gauge-field then results, at intermediate interaction depending a r X i v : . [ c ond - m a t . s t r- e l ] N ov on the strength of the spin-orbit coupling, in a first orderphase transition towards a non collinear magnetic order,which further evolves to magnetic texture at large inter-action, predicted by the effective Heisenberg model.The paper consists in two main parts: in Sec-tion II, we describe the theoretical framework (Fermi-Hubbard model, linear response theory, effective Heisen-berg model). In Section III (i) we show that thenon-interacting energies (band structure) always depicta nesting at the antiferromagnetic order resulting ina diverging spin-spin susceptibility at low temperatureand, thereby, an instability towards an antiferromagneticphase at small interaction; (ii) we provide the numer-ical mean-field results for different values of the spin-orbit coupling and the interaction, showing evidences forthe transition from the antiferromagnetic order to a non-collinear magnetic order. In particular, the skyrmionicphase is shown to already exist for moderate values ofthe interaction. II. METHODSA. model
In this paper we consider a system of spin 1/2 fermionicparticles in 2D square lattice and subjected to a syntheticgauge field (more will be said about the specific form ofthe gauge field later). The lattice spacing is equal tounity, setting both spatial and momentum scales. Thesystem Hamiltonian reads H tot = H kin + H µ + H int , (1)where H kin = − X j , j c † j ,s T s,s j , j c j ,s (2a) H µ = − X j µ n j , + µ n j , (2b) H int = U X j (cid:18) n j , − (cid:19) (cid:18) n j , − (cid:19) , (2c)where T are the tunneling matrices and are taken to betime independent in the following. c † j ,s , c j ,s are the usualfermionic creation and annihilation operators satisfying { c j ,s , c † j ,s } = δ j , j δ s,s , n j ,s = c † j ,s c j ,s is the density opera-tor, µ s the chemical potential, U the interaction strengthand s = 1 , H int = − U X j S j · S j , (3)where S a j = c † j ,s σ as,s c j ,s are the usual spin operators.To be more specific, when needed, we will consider thecase of the gauge fields corresponding to the spin-orbit coupling of the Rashba type, corresponding to the posi-tion independent tunnellings in the x and y directions: T x = t e − iασ y = t (cid:18) cos( α ) − sin( α )sin( α ) cos( α ) (cid:19) (4a) T y = t e iασ x = t (cid:18) cos( α ) i sin( α ) i sin( α ) cos( α ) (cid:19) . (4b)where t denotes the global strength of the hopping am-plitudes and we have used the labelling T x = T j , j +1 x and T y = T j , j +1 y . B. Non-interacting case - U = 0 In the case where the matrices T j , j are position inde-pendent, the Hamiltonian Eq.(2a) +Eq.(2b) is diagonal-ized in the momentum space: H = − X k C † k T k C k , (5)where C † k = ( c † k , , c † k , ) and T k = ( T x e − ik x + T † x e ik x + T y e − ik y + T † y e ik y + ˆ µ )ˆ µ = (cid:18) µ µ (cid:19) . (6)Finally, the matrix T k is diagonalized with a unitarymatrix U k , such that H = X k ,s (cid:15) k ,s d † k ,s d k ,s , (7)where ( d † k , , d † k , ) = ( c † k , , c † k , ) U † k .In the specific case of the spin-orbit coupling (4), oneobtains the following eigenenergies: (cid:15) k , , = − µ − t (cid:20) (cos k x + cos k y ) cos α ± q (sin k x + sin k y ) sin α + h (cid:21) , (8)where µ is the chemical potential and h is the spin im-balance defined through µ , = µ ± h (see Eq.(6)). Thesituation gets simpler in the limit of µ = h = 0, whichcorresponds to half filling for all α . C. Linear Response in small U limit Turning on the interaction, one expects the Fermi liq-uid phase to be unstable towards a magnetically orderedphase. This instability of the system can be captured us-ing the standard linear response theory (see Appendix Afor details), more precisely from the spin-spin suscepti-bility. Let us start with the non interacting Hamiltonianof the form H = H + H ext , where H is given by Eq.(5) and H ext = X j S a j B a j ( t ) = 1 N X q S a − q B a q ( t ) , where B is the external driving force. Using theFourier transform of the fermionic operators c k ,s = √ N P e − i k · j c j ,s , where N is the number of sites, onegets S a j = N P q e i q · j S a q and S a q = P k c † k ,s σ as,s c k + q ,s . The spin-spin susceptibility is diagonal in the momentumspace and reads χ a,b q ( ω ) = − i Z + ∞ d τ e − iωτ (cid:10)(cid:2) S a q ( τ ) , S b − q (0) (cid:3)(cid:11) , (9)where the thermal average and the time evolution aredone using the unperturbed Hamiltonian H . Therefore,the analytic expression for the susceptibility is easily ob-tained by diagonalizing H . After some manipulation,we find, that χ a,b q ( ω ) = 1 N X k ,s,s ( S a k , k + q ) s,s ( S b k + q , k ) s ,s n ( (cid:15) k ,s ) − n ( (cid:15) k + q ,s ) ω + (cid:15) k ,s − (cid:15) k + q ,s + iη → π ) Z d k X s,s ( S a k , k + q ) s,s ( S b k + q , k ) s ,s F s,s ( ω, k , q ) , (10)where n ( (cid:15) ) = (1 + e β(cid:15) ) − is the Fermi function. We haveintroduced F s,s ( ω, k , q ) = n ( (cid:15) k ,s ) − n ( (cid:15) k + q ,s ) ω + (cid:15) k ,s − (cid:15) k + q ,s + iη (11)and S a k , k + q = 12 U k σ a U † k + q , (12)where η is an infinitesimal convergence factor. In theusual situation of a Hamiltonian diagonal in the originalspin-space, the matrices U k are simply the identity andone recovers the standard spin-spin susceptibility.In what precedes, we have derived the susceptibilityEq.(10) for the non-interacting system. However, theinteraction among the fermions affects the spin-spin sus-ceptibility. This can be captured, in the random phaseapproximation (RPA) framework, by deriving, in a self-consistent way, the effective propagator for the spin fluc-tuations. This is equivalent to perform a mean-field ap-proximation to the interacting Hamiltonian .Lets recall the main step: Starting from the interac-tion Hamiltonian Eq.(3), the effect of interaction then amounts to an introduction of an effective driving forcegiven by H effext = H ext + H MFint = X q S b − q ( t )( B b q ) eff ( t ) , (13)where ( B b q ) eff ( t ) = (cid:2) B b q ( t ) − g h S b q i ( t ) (cid:3) and we have in-troduced g = 2 U/ h S q i ( ω ) = M − q ( ω ) χ q ( ω ) B q ( ω ) (14)where M ab q ( ω ) = δ a,b + 2 gχ a,b q ( ω ). M − q ( ω ) χ q ( ω ) istherefore the RPA susceptibility, whose singularitiesin the complex ω plane correspond to the vanishingeigenvalues of M q ( ω ). D. Large U limit - the Heisenberg Hamiltonian Following the method of , one obtains, in the large(repulsive) interaction U limit and at half-filling, the fol-lowing effective Heisenberg Hamiltonian, up to the sec-ond order in the t/U expansion: H = H c + X δ = x,y X X a = x,y,z J aδ S ai S ai + δ + D δ + · ( S i × S i + δ ) + + D δ − · ( S i × S i + δ ) − , (15)where ( S i × S i + δ ) + = ( S y S z , S z S x , S x S y )( S i × S i + δ ) − = − ( S z S y , S x S z , S y S x ) are the "positive" and "negative" part of the vector prod-uct. In the most general case, the coefficients J aδ , D pδ and H c , δ = x, y , a = x, y, z , p = 1 , ,
3, are quadratic func-tions of elements of the tunnelling matrices T . The gen-eral expression can be found in the appendix B. However,the situation will simplify considerably when consideringa specific case of spin orbit coupling of Rashba type, seeEq. (4): J xδ = 4 λ cos (2 α ) J yδ = 4 λJ zδ = 4 λ cos (2 α )for both spatial directions δ = x, y , D xδ + = D xδ − = 0 D yδ + = D yδ − = − λ sin (2 α ) D zδ + = D zδ − = 0for δ = x (tunnelling T x ) and D xδ + = D xδ − = − λ sin (2 α ) D yδ + = D yδ − = 0 D zδ + = D zδ − = 0for δ = y , (the tunnelling T y ), H c = − λ ,λ = t /U . This Hamiltonian is identical with Eq. (2)of , where one has to set their parameter λ = 1. In thefollowing, we take t equal to unity to set the energy scale. III. RESULTS & DISCUSSIONA. Small U limit
Here, we study the small interaction behaviour of thelattice gas at half-filling, µ = h = 0, by means of thesusceptibility, given by Eq.(10). The instability of thenon-interacting ground state is signalled by a divergenceof the DC ( ω = 0) RPA susceptibility (14), correspond-ing to a vanishing eigenvalue of M q (0). The latter corre-sponds to an eigenvalue of χ q (0) having the value − / g .In the following, we evaluate the susceptibility by numer-ical integration of (10) over the first Brillouin zone. Wewill discuss the general q value below, but it turns outthat the mean-field simulations indicate the onset of AFphase for small couplings U , corresponding to a value of q = q π = ( ± π, ± π ). For this specific value of q = q π ,one can see, using Eq.(8), that (cid:15) k + q π ,s = − (cid:15) k , ¯ s , where¯ s means the opposite spin to s . Moreover, in this case,the Fermi energy is E F = 0 and one gets nested Fermisurfaces for any value of the gauge field parameter α - theFermi surfaces together with the susceptibility integrandfunction F s,s (0 , k , ( π, π )), Eq.(11), are plotted in Fig.(1)for different values of α . The coincidence of the maxi-mum of the integrand function F is a specific feature of FIG. 1. (Color online) Plots of Fermi surfaces and suscepti-bility integrand F (see Eq.(11)) in the first BZ for q = ( π, π ).Rows: successive values of α = 0 , π/ , π/
2; Columns: (i)Fermi surfaces. The nesting of Fermi surfaces is indicated ex-plicitly for α = π/
3. For α = π/ F for spins ( s, s )=(2,1),(iii) F for spins ( s, s )=(1,2). The green and red colours in(ii) and (iii) represent the maximum of F . The Fermi surfacesfor given ( s, s ) in (ii) and (iii) are indicated by white lines.The coincidence of Fermi surfaces and the maximum of theintegrand F is not generic, but is specific for q π = ( ± π, ± π ). q π and is responsible for the the onset of AF phase forsufficiently low temperature, i.e. large β values. Indeed,we note, that in the case of nested Fermi surfaces, thesusceptibility possess ln β divergence.Moreover, for α = 0, the presence of Van Hove singulari-ties at the Fermi energy adds a leading divergence ln β .Therefore, we fit the susceptibility with a function χ fit = a ln β + b ln β + c. (16)The results are summarized in Table I and plotted inFig.(2). With respect to the preceding discussion, onecan distinguish two cases - α = 0 and α >
0. We verified,that for α = 0, our results agree with the theoreticalprediction given by Eq.(7) in (first two lines in TableI) χ q π , theor = − π ˜ t ln γ ˜ tπ β + C . (17)For α > α strictly positive), the absolute value of b decreases monotonically as the nesting of Fermi sur-faces decreases. Also, for α > a term becomessignificantly smaller than for the α = 0 case (no ln β divergence). α a × b c σ × π/
12 0.9 -0.065 -0.11 2. π/ π/ π/ π/
12 -0.85 -0.0052 -0.082 1.6 π/ a, b, c of Eq.(16) fordifferent values of the gauge field parameter α . σ is the usualdata variance, defined in the text. The first line correspondsto the theoretical prediction Eq.(17) and agrees well with thenumerical computation of the susceptibility Eq.(10) shown inthe second line. The data variance in the Table I reads σ = n P ni =1 | χ fit ,i − χ i | , where χ i stands for a minimal eigen-value at (given) q π and n = 14 is the number of simu-lated data points. As shown in Fig.(2), the parameter ofthe ln β divergence b increases monotonically (for α > β . This comes from the fact, that the contribu-tion to the susceptibility Eq.(10) is proportional to thearea of the Fermi surface (length of the 1D surface in ourcase), which decreases with increasing α and shrinks tothe Dirac points for α = π/ T = 0, athigher temperature the susceptibility develops minimaat q = q π . An example for a diagonal q , q = ( q d , q d ), q d ∈ [0 , π ] is shown in Fig.(3). Therefore, at this temper-ature, the linear response analysis predicts an instabilitytowards a different magnetic order. In addition, sincethese minima correspond to a finite value of the suscep-tibility χ m , the phase transition is predicted to occur ata finite interaction U = −
34 1 χ m . However, this predic-tion assumes a second order phase transition and at this(large) value of U , another magnetic order might have al-ready appeared. In addition, since the value of the mini-mum of the susceptibility is not much lower than the oneof the AF order, it might explain that, from a numericalpoint of view (finite size...), the onset of those non-AFphases at finite temperature and finite interaction havenot been observed yet. FIG. 2. (Color online) a) Plot of the minimal eigenvalue of χ against β for q = ( π, π ). The data points were calculatedusing Eq.(10) and the lines are the fits Eq.(16). The lowestcurve corresponds to α = 0 and the highest to α = π/ b against the gauge field parameter α > α ). The plot shows a monotonicdecrease of | b | - see text for details. a)b) FIG. 3. (Color online) Susceptibility eigenvalues vs. q on thediagonal of the Brillouin zone, q = ( q d , q d ). An example for α = π/ β = 8 . q = q π . B. Mean Field numerical simulation
To further investigate the properties of the system,we study the properties of the mean-field Hamiltonianground state. More precisely, at finite temperature, weminimize the mean-field free energy F MF = − β ln Z ,where Z is the partition function associated to the mean-field Hamiltonian H MF (See Appendix A for details): H MF = − X j , j c † j ,s T s,s j , j c j ,s − U X j h S j i · S j + 2 U X j h S j i · h S j i . (18)At half-filling, with a repulsive interaction, the total av-erage density is expected to remain fixed to unity, therelevant degrees of freedom being the average values ofthe spin operators S aj . The present mean-field decouplingrespects thus the SU (2) invariance of the interaction andallows for all possible magnetic orderings, in particularthose obtained in the large U limit from the effectiveHeisenberg model . From that point of view, eventhough other mean-field decoupling schemes are possi-ble , the present one is expected to capture qualitativelythe properties of the magnetic phases for different valuesof U and α .The numerical calculation has been performed on a36 ×
36 square lattice with periodic boundary conditionsfor different values of parameters α, U, β . On each lat-tice site, the three components of the spin h S j i are in-dependent mean-field parameters. Since the mean-fieldHamiltonian (18) is quadratic in the fermionic opera-tors, the free energy can be obtained by diagonalizing a2 N × N matrix, where N is the number of lattice sites.Even though the exact structure of the matrix is slightlydifferent from the one obtained in the BCS case ,there is a one to one mapping between the two situa-tions, namely a particle-hole transformation on one of thespecies. From a numerical point of view, the free energyis minimized using a mixed quasi-Newton and conjugategradient method; additional checks (e.g. different initialvalues of the spins) were performed to ensure that theglobal minimum has been reached. Finally, even thoughthe mean-field calculation captures the temperature de-pendence of the spin degrees of freedom, it only describesthe Mott transition, whereas the determination of thetrue critical temperature to a quasi-long range magneticorder, especially in the large interaction limit, amountsto taking into account the effects of terms beyond meanfield .The results are summarized in Table II. Only few val-ues of the gauge field and interaction are presented, sincethe focus of the paper is on the generic evolution of thesystem phase from the non-interacting situation to thelarge interacting limit. For the parameters under consid-eration, we have identified the following phases: • Anti-ferromagnet (AF), q = ( π, π ), • Spiral (SI), q = ( q, q ), TABLE II. (Color online) Summary of the MF results.
Leg-end:
AF (blue) - antiferromagnetic phase, SI (dark red), SI2(orange) - spiral phases, SkX (yellow) - Skyrmion crystal, seetext and Table III for definitions; For fixed values of α , weshow phases which occur as a function of β , the inverse tem-perature, and U , the interaction strength. For small values ofthe interaction, AF phase occurs for all α and for sufficientlyhigh values of β (low temperatures), in agreement with theresults predicted by the linear response theory. For increasingvalues of U , a phase transition occurs towards different phases(SI, SI2, SkX), depending on α . For other values of α , notshown here, we have found a similar scenario. For α ≥ . π ,the logarithmic divergence of the susceptibility with the tem-perature is very slow, such that the AF phase only appearsfor very low values of the temperature and, from a numericalpoint of view, is very difficult to observe with our method.Coexistence of different phases is indicated (either of two welldefined phases or of a well defined phase and an unknownphase (coex)). • Spiral 2 (SI2), q = ( q, π ), • Skyrmion crystal (SkX), q = ( π, ± π/
3) and q is modulo the periodicity of theBZ. As one can see from the Table II, an AF phase al-ways shows up first at low interaction; one can see thatwhen the transition normal-AF occurs for a low interac-tion, the phase only exists for very low temperature, abehaviour similar to the BCS situation. For other valuesof α , not shown here, we have found a similar scenario.For α ≥ . π , the logarithmic divergence of the suscep-tibility with the temperature is very slow, such that theAF phase only appears for very low values of the temper-ature and, from a numerical point of view, is very difficultto observe with our method.For larger interaction values, the AF order furtherevolves to a phase depicting a magnetic texture (spiral...),corresponding to peaks in | S ( q ) | away from the cornerof the BZ. The transition is of first order, since the AForder parameter doesn’t vanish at the transition point, FIG. 4. (Color online)
First row:
Density vector plots of spins h S a j i for different interaction strengths U = 4 . , . , . , ×
36 lattice. The real space z component of the spin S z j is color encoded, where thepink color corresponds to the maxima and the green color to the minima of S z respectively. The color scaling is relative toeach U independently. Second row:
3D plots of the spin density | S ( q ) | in the first Brillouin zone for interaction strengths U = 4 . , . , . , U = 4 . , . ,
10 correspond to AF, SI and SkX phases, which are clearly indicated by peaks of | S ( q ) | . U = 6 . α = 0 . π and β = 50. See also Table II and Appendix C for quantitative details. at which also a magnetic texture appears.Then, for larger interaction strength and dependingon the gauge-field parameter, the magnetic texture canfurther evolve to another phase, to reach finally the spinconfiguration predicted from the Heisenberg-model. It’smore difficult to determine the type of transition fromone phase to the other, but from the numerical data itseems to correspond to a smooth change in the location ofthe peaks of | S ( q ) | , i.e. a crossover within the numericalresolution of the simulation.Finally, one should notice that the critical value U c ( α )of the transition from the AF order is not expected todepict a smooth curve in the U − α plane, since, in themagnetic texture phase, the different mean-field Hamil-tonians H MF ( α ) have different translational symmetries,a situation similar to the Hofstadter-Hubbard model, de-picting fermions (or bosons) in an external magnetic field.Nevertheless, our numerics seems to suggest that thetransition from an AF order to a magnetic texture phaseincreases with the gauge-field parameter α , such that theAF order on the α = 0 line seems to be unstable towardsa magnetic structure phase in the presence of arbitrarysmall spin-orbit coupling.In Fig.(4), an example of a real space spin configu-ration h S a j i together with its Fourier transform | S ( q ) | is shown for α = 0 . π and β = 50. Different phases(AF, SI, SkX for U =4.5, 4.75, 10) or their coexistence(SI + SkX for U = 6.5) are clearly indicated by peaks of | S ( q ) | . As one can see, in the large U limit, one recovers q N ( q ) × N ( q ) N ( q ) · N ( q )AF ( π, π ) 0 0SI ( q, q ) = 0 0SI2 ( q, π ) = 0 0SkX ( π, ± π/
3) non-planar = 0TABLE III. Summary of magnetic phases. Different phasesare defined by the value of q at which | S ( q ) | peaks. Furtherdistinction of magnetic orders given by the values of N · N and N × N is shown. The skyrmion phase (SkX) corre-sponds to a pair of parameters ( N ( q ± ) , N ( q ± )) at inequiv-alent positions q ± in the Brillouin zone, with, in addition,non-collinear N ( q ± ) × N ( q ± ) vectors. See also Eq.(19) andthe text for details. the spin configuration expected from the effective Heisen-berg model, i.e. a 3 × S j · (cid:0) S j + x × S j + y (cid:1) . The magnetic orders canbe parametrized by the peak values of | S ( q ) | . Morespecifically, each peak q gives rise to a spin wave, whichcan be described as h S a j i = N a ( q ) cos q · j + N a ( q ) sin q · j (19)with further distinction of collinear, N ( q ) × N ( q ) = 0,and non collinear, N ( q ) × N ( q ) = 0, orders. The valuesof q , N ( q ) and N ( q ) are summarized in Appendix C. IV. CONCLUSION
In summary, we have studied the quantum phase tran-sitions of the Fermi-Hubbard model in a square latticeat half-filling in the presence of an effective spin-orbitcoupling. We have shown that at small interaction, thesystem always enters an AF order, then undergoes afirst order phase transition to a phase depicting mag-netic texture, and, eventually, reaches (at large interac-tion strength), the magnetic texture predicted by the as-sociated Heisenberg model.In addition to the half-filling situation presented here,a possible study will concern the doped case or theimbalanced population case, where more exotic mag- netic phases are expected to occur, possibly in compe-tition with the non-conventional superconductivity. Oneshould also take into account the effects of terms beyondmean field to determine properly the critical tempera-ture of the transition and to estimate the strength of thequantum fluctuations thus allowing for a better compar-ison with possible experimental results.Apart from the (magnetic) properties of the groundstate, it would be interesting to study the excitationsabove the ground state, in particular to describe thedynamical response of the system to external perturba-tions, like the (sudden) quenches of the interaction orthe gauge-field, which can efficiently be achieved in coldatomic gases experiments. M. Lewenstein et al., Adv. Phys. , 243 (2007). I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. , 885 (2008). W. Ketterle, and M. W. Zwierlein, in
Ultra-cold FermiGases , Proceedings of the International School of PhysicsEnrico Fermi, Varenna, 20-30 June 2006, Course CLXIV,edited by M. Inguscio, W. Ketterle and C. Salomon, IOSPress (Amsterdam), p. 95 (2007). Y. J. Lin, R. L. Compton, A. R. Perry, W. D. Phillips,J. V. Porto, and I. B. Spielman, Physical Review Letters , 130401 (2009). Y .-J. Lin, R.L. Compton, K. Jiménez-García, J.V. Porto,and I.B. Spielman, Nature (London) 462, 628 (2009). Y.-J. Lin, K. Jiménez-García, and I.B. Spielman, Nature(London) 471, 83 (2011). P. Wang, Z.-Q. Yu, Z. Fu, J. Miao, L. Huang, S. Chai, H.Zhai, and J. Zhang, Phys. Rev. Lett. , 095301 (2012). L. W. Cheuk, A. T. Sommer, Z. Hadzibabic, T. Yefsah,W. S. Bakr, and M. W. Zwierlein, Phys. Rev. Lett. 109,095302 (2012). J. Struck, C. Ölschläger, M. Weinberg, P. Hauke, J. Si-monet, A. Eckardt, M. Lewenstein, K. Sengstock, and P.Windpassinger, Phys. Rev. Lett. , 225304 (2012). K. Jiménez-García, L.J. LeBlanc, R. A. Williams, M. C.Beeler, A. R. Perry, and I. B. Spielman, Phys. Rev. Lett. , 225303 (2012). V. Pietilä and M. Möttönen, Phys. Rev. Lett. , 080403(2009). N. Goldman, D.F. Urban, D. Bercioux Phys. Rev. A ,063601 (2011). Michele Burrello and Andrea Trombettoni, Phys. Rev.Lett. , 125304 (2010). Y. Nishida and S. Tan, Phys. Rev. Lett , 170401(2008). Y. Nishida, Phys. Rev. A , 011605(R) (2010). G. Lamporesi et al., Phys. Rev. Lett. , 153202 (2010). W.-M. Huang, K. Irwin and S.-W. Tsai, Phys. Rev. A ,031603(R) (2013). M. Iskin and A. L. Subasi, Phys. Rev. A , 063628 (2010). J. Dalibard, F. Gerbier, G. Juzeli¯unas, and P. Öhberg,Rev. Mod. Phys. 83, 1523 (2011). W. S. Cole, S. Zhang, A. Paramekanti, and N. Trivedi,Phys. Rev. Lett. , 085302 (2012). Daniel Cocks, Peter P. Orth, Stephan Rachel, Michael Buchhold, Karyn Le Hur, and Walter Hofstetter, Phys.Rev. Lett. , 205303 (2012) and references therein. N. Goldman, A. Kubasiak, P. Gaspard, and M. Lewenstein,Phys. Rev. A , 023624 (2009). Z. Cai, X. Zhou, and C. Wu, Phys. Rev. A , 061605(R)(2012). Radić, J. and Di Ciolo, A. and Sun, K. and Galitski, V.Phys. Rev. Lett. , 085303 (2012). M. Sato, Y. Takahashi, and S. Fujimoto, Phys. Rev. Lett. , 020401 (2009). N. Goldman, W. Beugeling, C. Morais-Smith Euro-phys. Lett. , 23003 (2012). Jayantha P. Vyasanakere, Shizhong Zhang, and Vijay B.Shenoy Phys. Rev. B , 014512 (2011) Li Han and C. A. R. Sá de Melo Phys. Rev. A ,011606(R) (2012) A. Kubasiak, P. Massignan, M. Lewenstein, Euro-phys.Lett. , 46004 (2010). M. Iskin and A.L. Subasi, Phys. Rev. A , 043621 (2011). M. Iskin and A.L. Subasi http://arxiv.org/abs/1211.4020(2012). M. Iskin, http://arxiv.org/abs/1304.1473 (2013). K. Riedl, C. Drukier, P. Zalom, and P. Kopietz Phys. Rev.A , 063626 (2013). A. H. MacDonald, S. M. Girvin, and D. Yoshioka, Phys.Rev. B , 9753 (1988). Rear Earth Magnetism , edited by J. Jensen and A. R.Mackintosh (Clarendon Press, Oxford, 1991). Condensed Matter Field Theory , edited by A. Altlandand B. Simons (Cambridge University Press, ADDRESS,2010). David Pekker, Mehrtash Babadi, Rajdeep Sensarma, Niko-laj Zinner, Lode Pollet, Martin W. Zwierlein, and EugeneDemler, Phys. Rev. Lett. , 050402 (2011). H. Shimahara and S. Takada, J. Phys. Soc. Jap. , 1044(1988). In Eq.(16) C = 0 . γ ≈ . t = 1, we have toset the tunneling amplitude ˜ t = 2 in due to the factor2 of difference in the definition of eigenenergies. In addi-tion, for α = 0, the unitary matrices in the definition of S Eq.(12) reduce to the identities and we are left with X k ,s,s ( S a k , k + q ) s,s ( S b k + q , k ) s ,s F ( k ) = X k F ( k ) 14 X s,s σ ass σ bs s In this case, the last sum over elements of the Pauli matri-ces is 0 for a = b and 2 for a = b as can be easily verified.This gives purely diagonal susceptibility, as expected for α = 0. It also yields an additional factor of 2 between thedefinition Eq.(10) and the formula Eq.(17). H.J. Schultz, Phys. Rev. Lett. , 2462 (1990). Y. Dubi and Y. Avishai, Nature , 876 (2007). Y. Chen, Z.D. Wang, F.C. Zhang, C.S. Ting, Phys. Rev. B , 054512 (2009). J.R. Engelbrecht, M. Randeria, and C.A.R. Sá de Melo,Phys. Rev. B , 15153 (1997). S. Sachdev, Rev. Mod. Phys. , 913 (2003). Atom-Photon Interactions , edited by C. C. Tannoudji, J.Dupont-Roc, and G. Grynberg (John Wiley & Sons, NewYork, 1998).
Appendix A: Linear response of the spin systems
From the linear response theory , the expression for the spin-spin susceptibility χ a,b q ( τ ) reads: χ a,b q ( τ ) = − iθ ( τ ) (cid:10)(cid:2) S a q ( τ ) , S b − q (0) (cid:3)(cid:11) . (A1)The frequency domain susceptibility is given by the Fourier transform χ a,b q ( ω ) = R dτ e − iωτ χ a,b q ( τ ).We can now find an explicit analytical expression for χ a,b q ( ω ), given the Hamiltonian Eq.(2a - 2c). Namely we needto evaluate the thermal average (cid:10)(cid:2) S a q ( τ ) , S b − q (0) (cid:3)(cid:11) = Z − Tr (cid:20)(cid:0) S a q ( τ ) S b − q (0) − S b − q (0) S a q ( τ ) (cid:1) e − βH (cid:21) . (A2)In order to evaluate the trace and the time dependence, one needs to diagonalize H . As explained in the main text,this is achieved by going to momentum space and finding 2 × U k ( d k = U k c k ) such that H = X k ,s (cid:15) k ,s d † k ,s d k ,s . (A3)In the diagonal basis, the time evolution of the operators d k ,s is simple, such that one obtains readily the timeevolution of the spin operators: S a q ( τ ) = X k d † k ,s ( S a k , k + q ) s,s d k + q ,s e iτ ( (cid:15) k ,s − (cid:15) k + q ,s ) , (A4)where S a k , k + q = 12 U k σ a U † k + q . (A5)We proceed with the evaluation of the trace Eq.(A2). Lets start with the first part of the commutatorTr (cid:0) S a q ( τ ) S b − q (0)e − βH (cid:1) = X n i h n | S a q ( τ ) | n i h n | S b − q (0) | n i h n | e − βH | n i = X S a q ( τ ) n ,n S b − q (0) n ,n (cid:2) e − βH (cid:3) n ,n , (A6)where the sum runs over complete basis (i.e. momentum and spin) | n i i and we have used the fact, that H is alreadydiagonalized. Plugging the expression Eq.(A4) to Eq.(A6) one findsTr (cid:0) S a q ( τ ) S b − q (0)e − βH (cid:1) == Tr (cid:16) d † k ,s ( S a k , k + q ) s,σ d k + q ,σ d † k ,s ( S b k , k − q ) s ,σ d k − q ,σ e − iτ ( (cid:15) k ,s − (cid:15) k q ,σ ) e − β P (cid:15) p ,σ d † p ,σ d p ,σ (cid:17) = Tr (cid:16) d † k ,s ( S a k , k + q ) s,s d k + q ,s d † k + q ,s ( S b k + q , k ) s ,s d k ,s e − iτ ( (cid:15) k + q ,s − (cid:15) k ,s ) e − β P (cid:15) p ,σ d † p ,σ d p ,σ (cid:17) . (A7)0Following Eq.(A6), we put in the previous k = k + q and σ = s , σ = s . One obtainsTr (cid:0) S a q ( τ ) S b − q (0)e − βH (cid:1) == Z X k ,s,s ( S a k , k + q ) s,s ( S b k + q , k ) s ,s e − iτ ( (cid:15) k + q ,s − (cid:15) k ,s ) e − β(cid:15) k ,s (1 + e − β(cid:15) k ,s )(1 + e − β(cid:15) k + q ,s ) , q = 0= Z X k ,s = s ( S a k , k ) s,s ( S b k , k ) s ,s e − iτ ( (cid:15) k ,s − (cid:15) k ,s ) e − β(cid:15) k ,s (1 + e − β(cid:15) k ,s )(1 + e − β(cid:15) k ,s )+ X k ,s ( S a k , k ) s,s ( S b k , k ) s,s e − iτ ( (cid:15) k ,s − (cid:15) k ,s ) e − β(cid:15) k ,s − β(cid:15) k ,s , q = 0 . (A8)The trace of S b − q (0) S a q ( τ )e − βH is obtained in a similar way and yields a result identical to Eq.(A8) with exchange β(cid:15) k ,s ↔ β(cid:15) k + q ,s (i.e. the energies are exchanged only in the thermal terms including β ). We next notice, thate − β(cid:15) a − e − β(cid:15) b (1 + e − β(cid:15) a )(1 + e − β(cid:15) b ) = e − β(cid:15) a + 1 − − e − β(cid:15) b (1 + e − β(cid:15) a )(1 + e − β(cid:15) b )= 11 + e − β(cid:15) b −
11 + e − β(cid:15) a = (1 − n ( (cid:15) b )) − (1 − n ( (cid:15) a ))= n ( (cid:15) a ) − n ( (cid:15) b ) . (A9)This also holds for the degenerate case ( q = , last line in Eq.(A8)), in which we have directlye − β(cid:15) a − β(cid:15) a − e − β(cid:15) b − β(cid:15) b = n ( (cid:15) a ) − n ( (cid:15) b ) . (A10)We obtain the result for the trace of the commutator (cid:10)(cid:2) S a q ( τ ) , S b − q (0) (cid:3)(cid:11) = X k ,s,s ( S a k , k + q ) s,s ( S b k + q , k ) s ,s ( n ( (cid:15) k ,s ) − n ( (cid:15) k + q ,s ))e − iτ ( (cid:15) k + q ,s − (cid:15) k ,s ) . (A11)The only time dependent factor is the oscillating exponential and we can thus directly compute the time integral inthe definition of the susceptibility Eq.(A1) − i Z ∞−∞ d τ e iωτ θ ( τ )e − iτ ( (cid:15) k + q ,s − (cid:15) k ,s ) = − i Z ∞ d t e iωτ − ητ e − iτ ( (cid:15) k + q ,s − (cid:15) k ,s ) = 1 ω + (cid:15) k ,s − (cid:15) k + q ,s + iη , (A12)where we have added the infinitesimal convergence factor η . Plugging this back to the defining relation for thesusceptibility Eq.(A1), we obtain the final result for the susceptibility of the non interacting system χ a,b q ( ω ) = 1 N X k ,s,s ( S a k , k + q ) s,s ( S b k + q , k ) s ,s n ( (cid:15) k ,s ) − n ( (cid:15) k + q ,s ) ω + (cid:15) k ,s − (cid:15) k + q ,s + iη . (A13) Susceptibility in the interacting MF model
We have derived the susceptibility Eq.(A13) for the non-interacting system subjected to a small external driving force.We will now use this result to find a susceptibility of the interacting system described by the mean-field theory. Letsrecall the interaction Hamiltonian Eq.(3): H int = − U X j S b j ( t ) S b j ( t ) = − U N X q S b − q ( t ) S b q ( t ) MF ≡ − gN X q h S − q ih S q i + h S − q i ( S q − h S q i ) + ( S − q − h S − q i ) h S q i + O (( δS ) )= − gN X q S − q h S q i − h S − q ih S q i = H MFint , (A14)1where we have introduced the coupling strength g = 2 U/
3. In the last line of the preceding equation, the last termdoes not contribute to the spin dynamics and is normalized out in the computation of the spin average values. Wethus drop this term. We obtain the effective external driving Hamiltonian H effext = H ext + H MFint = 1 N X q S b − q ( t ) (cid:2) B b q ( t ) − g h S b q i ( t ) (cid:3) = 1 N X q S b − q ( t )( B b q ) eff ( t ) . (A15)We can thus see, that the inclusion of the MF interaction amounts to replacing the external driving B by thenew effective driving B eff . One then follows the same procedure as in the non interacting case. Since we aremainly interested in the frequency response of the system, we will use the defining formula connecting the frequencycomponents of the spins with the driving through the susceptibility < S a q ( ω ) > = χ a,b q ( ω )( B b q ) eff ( ω ) . (A16)One can easily check, that the frequency components of the effective driving are given by the Fourier transform of itsparts, ( B b q ) eff ( ω ) = Z d t e iωt ( B b q ( t ) − g ¯ S b q ( t ))= B b q ( ω ) − g ¯ S b q ( ω ) . (A17)Therefore, one obtains h S a q i ( ω ) = X b χ a,b q ( ω ) B b q ( ω ) − g h S b q i ( ω ) , (A18)which after a straightforward manipulation yields the equation for the average value of the spin operators h S q i ( ω ) = M − χ B M a,b = δ a,b + 2 gχ a,b q ( ω ) , (A19)where we merely rewrote the equation Eq.(A18) in the symbolic matrix notation. From here, one can obtain theinformation about the critical value of the coupling strength g (and thus U ) from the divergences of the average ofthe spin operators, i.e. when the matrix M becomes singular. Appendix B: Effective Heisenberg model in the large U limit We now restrict our interest to the regime with strong repulsion, high U . In this regime, the ground state of thegrand canonical ensemble has single occupation at each site. Moreover it would cost an energy of order of U toincrease the double occupancy by one. This regime can be described by the method of effective Hamiltonian, which issuitable for systems with well separated energy manifolds . The energy manifolds are separated by U and we wouldlike to evaluate the effect of the coupling between the ground state manifold and the higher lying manifolds. Thiscoupling results in the perturbation of the bare energy levels in the ground state manifold. In this section, we presentthe treatment used in and (page 38).In the following we consider a situation at half filling µ = µ = 0. When we multiply the kinetic HamiltonianEq.(2a) by n j , ¯ s + h j , ¯ s = 1 from the left and by n j , ¯ s + h j , ¯ s = 1 from the right, we obtain H kin ≡ T = T + T − + T , (B1)where T = − X n j , ¯ s c † j ,s T s,s j , j c j ,s n j , ¯ s + h j , ¯ s c † j ,s T s,s j , j c j ,s h j , ¯ s T − = − X h j , ¯ s c † j ,s T s,s j , j c j ,s n j , ¯ s T = − X n j , ¯ s c † j ,s T s,s j , j c j ,s h j , ¯ s , (B2)where n and h denote the particle and hole number operators respectively and ¯ s, ¯ s denote the spin orthogonal to s, s - e.g. ¯ s is spin up for s spin down and vice versa. The sums in Eq.(B2) run over nearest neighbours h jj i and spins2 s, s . Denoting the interaction Hamiltonian Eq.(2c) as V , one can easily verify the following commutation relations[ V, T ] = 0[ V, T − ] = − U T − [ V, T ] = U T , (B3)which can be summarized as [ V, T m ] = mU T m . (B4)We are now ready to proceed with the effective Hamiltonian derivation. We wish to rewrite the current Hamiltonian H = V + T as H eff = e iS H e − iS = H + [ iS, H ] + 12! [ iS, [ iS, H ]] + . . . = ∞ X n =0 n ! [ iS, H ] ( n ) , (B5)where S is some Hermitian matrix (we require the transformation to be unitary). [ iS, H ] ( n ) denotes the n -times nestedcommutator and [ iS, H ] (0) = H . The matrix S can be determined in the following way. Lets write S as S = λS + λ S + . . . = ∞ X n =1 λ n S n , (B6)where λ = 1 /U is our small parameter around which we will do our perturbative expansion. The elements of thematrix S k can be determined in an iterative way by requiring, that after the unitary transformation up to the order k in the expansion Eq.(B5) all terms bringing the energy out of the ground state manifold have to vanish. We alsodenote S ( k ) = k X n =1 λ n S n (B7)Important thing to note is that the energy ratio between V and T is of order λ − . It is thus more transparent torewrite the Hamiltonian as H = V + λ ˜ T , where V and ˜ T = λ − T are now contributions of the same order. With thisnotation, the commutator Eq.(B4) can be written as λ (cid:2) V, ˜ T m (cid:3) = m ˜ T m . (B8)In the first order in λ we have from Eq.(B5) H eff = H + [ iλS , H ] = V + λ ( ˜ T + ˜ T − + ˜ T ) + λ [ iS , V ] + . . . . (B9)The terms which bring one from the ground state manifold are the terms changing the double occupancy, i.e. ˜ T − and ˜ T . In order to cancel these terms with the commutator, one gets iS = λ ( ˜ T − ˜ T − ) . (B10)We will now generalize this procedure in an iterative way to arbitrarily high order in λ . Lets define a Hamiltonian oforder k + 1 as H ( k +1) = e iS ( k ) H e − iS ( k ) = ∞ X n =0 n ! h iS ( k ) , H i ( n ) . (B11)As a matter of example, lets take k = 1 H (2) = H + λ [ iS , H ] + λ λ [ iS , [ iS , H ] + O ( λ )= H + λ [ iS , V ] + λ (cid:2) iS , ˜ T (cid:3) + λ iS , [ iS , V ] + O ( λ ) . (B12)3The last two terms in the second line are of the same order λ since S ∝ V ∝ λ − . The idea is now to use thisHamiltonian to find the explicit form of the next elements in the expansion of iS , namely the element iS . This canbe done as follows. The third order Hamiltonian can be written as H (3) = H (2)0 + H (2)ch + λ [ iS , V ] + O ( λ ) , (B13)where we have decomposed the second order Hamiltonian into a part which does not change the double occupancy H (2)0 and the one which changes the double occupancy, H (2)ch . We require, that the lowest order term of iS cancelsthe double occupancy changing term. This is the general procedure for an arbitrary order expansion, so that H ( k +1) = H ( k )0 + H ( k )ch + λ k [ iS k , V ] + O ( λ k +1 ) . (B14)In order to find the explicit form of iS k , we will need the following relationships. Lets denote the product of tunnellingoperators as ˜ T k ( m ) ≡ ˜ T m . . . ˜ T m k . (B15)The commutator with V then reads λ (cid:2) V, ˜ T k ( m ) (cid:3) = M k ( m ) ˜ T k ( m ) , (B16)where M k ( m ) = P i m i . Looking at Eq.(B14), one can write the last two terms as H ( k )ch + λ k [ iS k , V ] = λ k − X m C k ( m ) ˜ T k ( m ) + λ k [ iS k , V ] = 0 . (B17)The last equality can be achieved by noting that λ (cid:2) V, ˜ T k ( m ) (cid:3) · λ k − C k ( m ) M k ( m ) = λ k − C k ( m ) ˜ T k ( m ) , (B18)so that we finally obtain iS k = λ k X m C k ( m ) M k ( m ) ˜ T k ( m ) . (B19)Note that since M k ( m ) = 0 for double occupancy changing T k , we can safely divide by it.We have implemented the above described iterative procedure in Mathematica for general tunnellings T s,s j , j .Up to the second order in the 1 /U expansion, one obtains the following effective Heisenberg Hamiltonian ( δ = x, y are the spatial directions in the 2D lattice) H = H c + X δ = x,y X X a = x,y,z J aδ S ai S ai + δ + D δ + · ( S i × S i + δ ) + + D δ − · ( S i × S i + δ ) − , where ( S i × S i + δ ) + = ( S y S z , S z S x , S x S y )( S i × S i + δ ) − = − ( S z S y , S x S z , S y S x )are the "positive" and "negative" part of the vector product. In the most general case, the coefficients read: J xδ = 2 λ (cid:16) T , δ T , δ ∗ + T , δ T , δ ∗ + T , δ T , δ ∗ + T , δ T , δ ∗ (cid:17) J yδ = 2 λ (cid:16) T , δ T , δ ∗ − T , δ T , δ ∗ − T , δ T , δ ∗ + T , δ T , δ ∗ (cid:17) J zδ = 2 λ (cid:16) T , δ T , δ ∗ − T , δ T , δ ∗ − T , δ T , δ ∗ + T , δ T , δ ∗ (cid:17) D δ + = 2 iλ (cid:16) − T , δ T , δ ∗ + T , δ T , δ ∗ + T , δ T , δ ∗ − T , δ T , δ ∗ (cid:17) D δ + = 2 λ (cid:16) T , δ T , δ − ∗ T , δ T , δ ∗ + T , δ T , δ ∗ − T , δ T , δ ∗ (cid:17) D δ + = 2 iλ (cid:16) T , δ T , δ ∗ + T , δ T , δ ∗ − T , δ T , δ ∗ − T , δ T , δ ∗ (cid:17) D δ − = 2 iλ (cid:16) − T , δ T , δ ∗ + T , δ T , δ ∗ + T , δ T , δ ∗ − T , δ T , δ ∗ (cid:17) D δ − = 2 λ (cid:16) − T , δ T , δ ∗ − T , δ T , δ ∗ + T , δ T , δ ∗ + T , δ T , δ ∗ (cid:17) D δ − = 2 iλ (cid:16) T , δ T , δ ∗ − T , δ T , δ ∗ + T , δ T , δ ∗ − T , δ T , δ ∗ (cid:17) H c = 2 λ (cid:16) − T , δ T , δ ∗ − T , δ T , δ ∗ − T , δ T , δ ∗ − T , δ T , δ ∗ (cid:17) Appendix C: Magnetic orders
In Table IV we list the details of the magnetic orders parametrized by Eq.(19) and summarized in Table III. Thetables have the following format { α, U, β } , Phase q x /π q y /π N x N x N y N y N z N z { α, U, β } = { } , AF q x /π q y /π N x N x N y N y N z N z − . × − − . × − . × − − . × − { α, U, β } = { } , SI2 q x /π q y /π N x N x N y N y N z N z − − . × − . × − − . × − − . × − − . × − . × − − . × − − . × − { α, U, β } = { } , AF q x /π q y /π N x N x N y N y N z N z − . × − − . × − . × − − . × − { α, U, β } = { } , SI q x /π q y /π N x N x N y N y N z N z −
56 59 − . × − − . × − . × − . × − − . × − . × −
356 59 − . × − . × − − . × − . × − − . × − − . × − { α, U, β } = { } , SI q x /π q y /π N x N x N y N y N z N z −
12 89 . × − . × − − . × − − . × − . × − − . × −
112 89 . × − − . × − . × − − . × − . × − . × − { α, U, β } = { } , SI q x /π q y /π N x N x N y N y N z N z −
12 89 . × − . × − − . × − − . × − . × − − . × −
112 89 . × − − . × − . × − − . × − . × − . × − { α, U, β } = { } , AF q x /π q y /π N x N x N y N y N z N z − . × − − . × − . × − − . × − { α, U, β } = { } , SI q x /π q y /π N x N x N y N y N z N z −
718 79 − . × − . × − . × − − . × − . × − . × − − . × − . × − − . × − . × − − . × − − . × − { α, U, β } = { } , AF q x /π q y /π N x N x N y N y N z N z − − . × − − . × − { α, U, β } = { } , SI q x /π q y /π N x N x N y N y N z N z
229 1318 − . × − − . × − − . × − − . × − − . × − . × − − . × − . × − − . × − . × − . × − . × − { α, U, β } = { } , SkX q x /π q y /π N x N x N y N y N z N z − − . × − . × − . × − . × − − − . × − − . × − . × − − . × − − . × − . × − . × − . × − − . × − . × − . × − . × − TABLE IV. Magnetic order parameters N i = ( N xi , N yi , N zi ), i = 1 , q = ( q x , q y ),together with the parameters { α, U, β } for phases given in Table II. In the table, we provide the data for magnetic orderparameters only for q y >
0, since the values for q y < q → − qq