Phases of attractive spin-imbalanced fermions in square lattices
PPhases of attractive spin-imbalanced fermions in square lattices
Simone Chiesa and Shiwei Zhang Department of Physics, College of William & Mary, Williamsburg, VA 23188, USA
We determine the relative stability of different ground-state phases of spin-imbalanced popula-tions of attractive fermions in square lattices. The ground state is determined within Hartree-Fock-Bogoliubov theory, with care taken to remove finite size effects. The phases are systematicallycharacterized by the symmetry of the order parameter and their real- and momentum-space struc-tures. For quarter- to half-filled lattices, where the Fermi surfaces are most distorted from theirspherical counterpart in the continuum, we predominantly find unidirectional Larkin-Ovchinikov-type phases. We discuss the effect of commensuration between the ordering wave vector and thedensity imbalance, and describe the mechanism of Fermi surface reconstruction and pairing for var-ious orders. A robust supersolid phase exists when the ordering wave vector is diagonally directed.Charge and pairing order coexist, rather than competing, and are responsible for the opening ofthe gap on different portions of the Fermi surface. A variational determination of the correct pairmomentum of the Larkin-Ovchinikov phases shows that phase separation does not occur in theconsidered regime of density and magnetization.
PACS numbers:
INTRODUCTION
There has been a surge of interest in the possibility ofrealizing unconventional fermionic superfluids using coldatomic gases. Amongst the many possibilities offered bythe highly tunable Hamiltonians available in cold-atomsexperiments, the simplest remains that of unequal pop-ulations of two hyperfine states in the presence of at-tractive interaction. The theoretical study of such sys-tems dates back to Fulde and Ferrel (FF) [1] and Larkinand Ovchinikov (LO) [2] who independently suggestedthat the mismatch between the Fermi surfaces of the twospecies could result in the formation of a condensate offinite-momentum pairs. Atomic gases offer a direct routeto the realization of FFLO phases, circumventing most ofthe difficulties of solid state systems, thanks to the possi-bility of controlling independently the density of the twospecies, the absence of disorder and, most importantly,the ability to engineer strong interactions. In spite ofthis, the existence of an FFLO phase in three dimensionshas been argued to be confined to a small range of inter-action strengths and polarizations[3], and detection hasremained elusive.The possibility of using optical lattices has been sug-gested by several authors[4, 5] as a key ingredient to ob-serve FFLO-type states. The best empirical indicationthat this may be the case is provided by experiments onstrongly-correlated-electrons materials and the fact that,when doped, these systems show a tendency toward for-mation of inhomogeneities in the form of spin, chargeand, possibly, pairing density waves. The relevance ofthese experiments to the properties of attractive fermionsin optical lattices stems from the belief that, in bothcases, the essential physics can be captured by a one-band Hubbard model: with an on-site repulsive inter-action for many of the electronic systems and an on-site attraction in an optical lattice. The attractive and repul-sive cases are mapped into each other by a particle-holetransformation[4] and the presence of spin-texture in thedoped repulsive case translates into the occurrence of amodulated superfluid in an imbalanced population of at-tractive fermions, i.e. an FFLO phase. This is reinforcedby recent quantum Monte Carlo results [6] on the two-dimensional repulsive model showing spin-density waveswith long wavelength modulation.Despite this mapping and several works addressingthe existence of a possible FFLO phase [7–9], the na-ture of the ground state phases in a spin-imbalancedtwo-dimensional optical lattice remains largely undeter-mined. On the one hand, information on the repulsivemodel is entirely confined to the case of unpolarized sys-tems, which maps into the attractive case at half-filling: n ↑ + n ↓ = 1; for the case of imbalanced fermionic pop-ulation, one is interested in the more general case of apolarized system and arbitrary density. On the otherhand, works addressing the physics in the lattice haveeither focused on the single plane-wave form of the or-der parameter[7, 8] or on selected states[9, 10] (e.g. infixed size supercells) because of the challenge of remov-ing large finite-size effects. Such restrictions can bias theresult and lead to, for example, an incorrect pair momen-tum. The accurate determination of the spatial structureof the order parameter is also indispensable to addressingthe issue of phase separation.In this work we establish the correct form that FFLOphases have in the thermodynamic limit on a square lat-tice, and show that a proper determination of the leadingpairing wave vectors in the ordered state leads to a char-acterization of different physical regimes based on theproperties of the nodal (excess) particles. Small to mod-erate interaction strengths are considered, where mean-field theory is expected to capture the correct physics. In a r X i v : . [ c ond - m a t . qu a n t - g a s ] O c t regime of density and polarization where the presence ofa lattice alters most dramatically the shape of the Fermisurfaces, we find unidirectional order and the existenceof three distinct phases: 1) a nodal metallic state char-acterized by the presence of Fermi arcs, 2) a nodal bandinsulator where the densities of excess particles and nodallines are equal and 3) a charge density wave that resultsin a robust supersolid phase obtained when the order-ing wave vector is directed along the diagonal direction.Once the variational search for the optimal pairing wavevector includes the proper symmetry of the order param-eter, the Hartree-Fock-Bogoliubov ground state does notphase separate in the parameter regime of interest. HARTREE-FOCK-BOGOLIUBOV THEORY
Results presented in this work are obtained usingHartree-Fock-Bogoliubov theory so that modulations incharge, spin and pairing are all handled on the same foot-ing. The starting Hamiltonian reads H = − t (cid:88) (cid:104) ij (cid:105) σ c † iσ c jσ − U (cid:88) i (cid:0) n i ↑ n i ↓ − µn i − h m i (cid:1) , (1)where c iσ are fermionic annihilation operators of spin σ on site i , n iσ = c † iσ c iσ , n i = n i ↑ + n i ↓ and m i = n i ↑ − n i ↓ .In order to accommodate the inhomogeneities, the cal-culations are performed on supercells whose shape isdictated by the symmetry of the targeted phase. Thesupercell is characterized by two basis vectors, L and L , whose components are integers. Once the super-cell shape is chosen, we define Bloch states as c j ( k ) ∝ (cid:80) L c j + L exp (cid:2) ik · L (cid:3) where L = n L + n L , and k is avector that varies freely within the super-lattice first Bril-louin zone. Then, using these states and the mean-fieldapproximation, the Hamiltonian decouples into a sum of k -dependent pieces, H = (cid:80) k H ( k ), of the form H ( k ) = [ c †↑ c ↓ ] (cid:20) H ↑ ( k ) ∆∆ † − H T ↓ ( G − k ) (cid:21) [ c ↑ c †↓ ] T , (2)where c ↑ and c ↓ represent an array (row) of operators, { c i ↑ ( k ) } and { c i ↓ ( G − k ) } respectively, with index i run-ning over the N sites of the supercell, and G is definedbelow. H and ∆ are N × N matrices with elements[ H σ ( k )] ij = − t ij ( k ) + δ ij ( D iσ − µ − s σ h/ ∆ ] ij = δ ij ∆ i . (3)In the above, t ij ( k ) = (cid:80) L exp( ik · L ) t i,j + L , s ↑ / ↓ = ± D iσ , ∆ i , µ and h are determined by the requirementthat the Free energy F = (cid:104) H (cid:105) − T S is a minimum for thetarget average densities n σ . This amounts to imposing U Order parameter0.080.090.10 2.0 2.2 2.4 | q | /π U Energy ( t × − ) 0.00.51.01.52.0 2.5 3.0 3.5Kin. En. SpiralLinearChkbd FIG. 1: Left Panel: local order parameter, max i |(cid:104) c i ↑ c i ↓ (cid:105)| ,and leading wavevector (inset) versus U for (cid:104) c i ↑ c i ↓ (cid:105) ∝ e iq x · r i (Spiral), (cid:104) c i ↑ c i ↓ (cid:105) ∝ cos( q x · r i ) (Linear), (cid:104) c i ↑ c i ↓ (cid:105) ∝ cos( q x · r i ) + cos( q y r i ) (Chkbd) with q x = | q | (1 ,
0) and q y = | q | (0 , n = 0 .
95 and m/n = 0 .
1. Kinetic energies are in units of t . the following self-consistency conditions D i − σ = − U (cid:90) dk (cid:104) c † iσ ( k ) c iσ ( k ) (cid:105) , ∆ i = − U (cid:90) dk (cid:104) c i ↓ ( k ) c i ↑ ( k ) (cid:105) ,n σ = N − (cid:88) i (cid:90) dk (cid:104) c † iσ ( k ) c iσ ( k ) (cid:105) , (4)with expectation values evaluated in the supercell. Notethat, although we target specific densities, the mean-fieldapproach works in the grand-canonical ensemble.The variable G in Eq. (2) is a vector such that θ = G · L gives the twist angle of the pairing order parameter un-der translation by L . For collinear phases, θ is either0 or π , giving periodic or anti-periodic boundary condi-tions on (cid:104) c † i ↑ c † i ↓ (cid:105) . Note that uniform phases with a spiralorder parameter (FF type) amount to the choice of a one-site cell and a G equal to the wave-vector of the spiralmodulation. Thus H ( k ) is specified by a 2 × SYMMETRY OF THE ORDER PARAMETERDetermination
Given a set of symmetry-equivalent q vectors and as-suming a continuous phase transition, linear responsetheory can be used to show that the onset of instabil-ities of the form ∆ i = (cid:88) q ∆ q e iq · r i , (5) -0.3-0.2-0.10.00.10.20.3 (a) Pairingamplitude-0.3-0.2-0.10.00.10.20.30 10 20 30 40site position(b) -0.30-0.25-0.20-0.15-0.10-0.050.00 0.960.981.001.021.04(a) charge fi spin ‹ -0.30-0.25-0.20-0.15-0.10-0.050.000 10 20 30 400.920.940.960.98site position(b) (a) A › (w) A fl (w) -3.0 -1.5 0.0 1.5 3.0 w / t (b) (a) |(cid:209) n k › ||(cid:209) n k fl | (b)0.91.00.8 0.9 1.0density l A0.91.00.8 0.9 1.0densityB
FIG. 2: Local properties and gradient of the momentum distribution for n ↑ − n ↓ = − .
05 at two densities: (A) n = 0 . n = 0 .
93 (bottom). In both cases, the finite-momentum of the pair, q , is in the (1 , λ = mπ/q as a function of density; case A belongs to the commensurate regimewith unit density of excess-spin particles per node, while case B is in the incommensurate regime. The second column reportthe spin, (cid:104) n i ↑ − n i ↓ (cid:105) , and charge, (cid:104) n i ↑ + n i ↓ (cid:105) , density profiles in the direction of q . The third column is the local density ofstates measured at site 0. In the fourth column, each panel shows two halves of the Fermi surface, for the ↑ (minority) and ↓ (majority) spins, respectively. The arrows in the figure indicate the pairing construction and applies to every pair across theFS, leading to one common q . must happen at exactly the same value of U , regard-less of the choice of ∆ q . Below such value, mean fieldtheory is guaranteed to return a normal, spin-polarizedFermi liquid. In order to determine the correct form oforder parameter, we proceed as follows. We first de-termine U c and the associated non-zero wave-vector q c using the single plane-wave form as this allows for aquick exploration of phase space. We find that q c is di-rected along any of the four, symmetry-equivalent direc-tions ( ± , , ±
1) and, therefore, any linear combi-nation of the four associated plane waves is a candidateground state order parameter just above U c . To resolvewhich one leads to the largest lowering of energy, we pro-ceed by solving the mean-field equations for the threecases corresponding to spiral (∆ i ∝ exp( iq x · r i )), uni-directional (∆ i ∝ cos( q x · r i )) and checkerboard (∆ i ∝ cos( q x · r i )+cos( q y · r i )) pairing density wave and track theevolution of | q | as a function of U . Explicit calculations inlarge simulations cells on the repulsive model (with sim-ulated annealing starting with random initial fields[11])have shown that instabilities involving q -vectors in differ-ent, non-equivalent directions, which the above approachwould miss, are unlikely to occur in the range of U con-sidered here. Dependence on density and polarization
Mean-field results[11] on the repulsive model and theparticle-hole transformation relating the attractive andrepulsive models imply the existence of the followingproperties at half-filling, i.e., when the average parti-cle density is precisely one fermion per site: 1) a crit-ical U exists such that, above it, the system developsa phase with an inhomogeneous order parameter 2) thepairing order parameter is characterized by a wave vector | q | = mπ where m = n ↑ − n ↓ is the average magnetiza-tion 3) both fermionic species have a gap in their singleparticle spectrum 4) as U grows larger there is a transi-tion from a pair-density-wave in the (1 ,
0) state to one inthe (1 ,
1) direction 5) order can be arbitrarily broken intoa charge or a pairing instability or a combination of thetwo. This last point is a consequence of the possibilityof breaking spin symmetry in any of the three equivalentdirections in the repulsive case. It implies that chargeand pairing orders compete, in the sense that the largerone is, the smaller the other must be.Although there have been studies on several aspects ofthe physics away from half-filling[10], the correct leadingpairing wave vector (which requires a scan through super- m/n n = 0 . n = 0 . n = 0 . . U c = 1 . U c = 1 . U c = 2 . q c = 0 . π q c = 0 . π q c = 0 . π Linear; Linear Linear; Linear Chkbd; Linear0 . U c = 3 . U c = 3 . U c = 2 . q c = 0 . π q c = 0 . π q c = 0 . π Linear; Linear Linear; Linear Chkbd; ChkbdTABLE I: Critical parameter at the onset of order for a fewdensities ( n ) and polarizations ( m/n ). The wavevector is di-rected along (1,0). The two entries in the last line in eachtable cell describe, respectively, the type of order just above U c and at U = 4 (see Fig.1’s caption for detail). cell sizes) has not been determined. This has preventeda characterization of the nature of such phases. As a re-sult the order of possible transitions and the related pos-sibility of phase separation have not been resolved. Weaddress these questions here using the strategy outlinedin the previous subsection. A representative example isgiven in Figure 1 for the case with m ≡ n ↑ − n ↓ = − . n = 0 . n = 0 . , .
75 and 0 . m/n = 0 . .
4, so as to have a rather complete picture closeand away from half-filling, at small and large polariza-tion and in an interaction range that extends from U c upto U = 4.The results are summarized in Table I. They are consis-tent with earlier results addressing the physics of FFLOphases in the continuum, which found checkerboard orderclose to U c thanks to an expansion in powers of the orderparameter[12], and show that for n ≥ .
75 lattice effectsare strong enough to return unidirectional order inde-pendently of polarization or proximity to U c , as knownto happen at and around n = 1 .
0. Given the difficulty inobserving unconventional pairing in the continuum, thelatter is the density regime where an experimental real-ization of the FFLO state could be more feasible. Wewill therefore focus on it in the following.
CHARACTER OF THE NODAL PHASESMetal and band insulator at weak coupling
Figure 2 characterizes the FFLO phase at U = 3 t and for m = 0 .
05 in two qualitatively different densityregimes. In particular, for n > .
95, the wave vector isprecisely determined by the magnetization via the samerelation holding at half-filling, q = mπ . This commen-surate regime is characterized by a density of one excess particle per node (of the order parameter) and conse-quent band-insulating behavior along the node. This ismost clearly seen in the local density of states at the node,with both species showing a gap at the Fermi energy.Correspondingly the gradient of the momentum distri-bution shows no sharp lines indicative of the existenceof a Fermi surface. The commensurate regime here has,however, some important distinctions from half-filling.First, the interchangeability between pairing and chargeorders is broken as soon as the average density deviatesfrom n = 1, and pairing emerges as the dominant order.Second, the density is not perfectly uniform, as it is forthe purely superfluid phase at half-filling. The density isinstead characterized by a weak modulation that reflectsthe different degrees of localization of the Andreev statesfor the two spin species. The density profile shows weakpeaks at the nodes indicating that the majority-speciesnodal states have stronger localization. A similar phasewas found in the context of a two-dimensional array oftubes[13].In the second density regime, n < .
95, the majorityspin species develop a finite density of states at the nodesof the order parameter. Fermi arcs appear in |∇ n k ↓ | in the form of sharp lines. The arc in Fig. 2.B is ofperfectly one-dimensional nature, indicating a completedecoupling between the metallic states living at differ-ent nodes. At larger polarization, the arcs will morestrongly resemble the underlying Fermi surface of thenon-interacting species. Finally, the density profile showsminima at the nodes, rather than the maxima observedin case A, as a direct consequence of the gradual emp-tying of the majority-spin Andreev bands visible in thelocal density of states. This is a potentially importantexperimental signature as it allows the characterizationof the nature of the nodal phase, metallic or insulating,via a static local property instead of a collective prop-erty such as the distance between nodes or the presenceof Fermi arcs. Supersolid phase at intermediate coupling
Next we consider the phases at larger interactionstrengths U . Results on the repulsive model[11] indi-cate that the ordering wave vector switches to the (1 , m and U > t . We focuson U = 4 t for an extensive and systematic examinationof the properties of these phases. Indeed unidirectionalLO states along the diagonal direction are found; in ad-dition, we find that the system can accommodate bothpairing and charge orders as non-competing instabilitiesaway from half-filling. The real- and momentum-spaceproperties of the phases are shown in Fig. 3. In contrastto the (1,0)-direction states, pairing here is achieved bya more radical reconstruction of the Fermi surface. Thenon-interacting surface of the majority spin is stretched, n › +n fl =0.95 ‹ charge fi pairing |(cid:209) n k › ||(cid:209) n k fl | n › +n fl =0.85 n › +n fl =0.75 FIG. 3: Left panels: Evolution of charge (cid:104) n i ↑ + n i ↓ (cid:105) and pair-ing order as density is reduced, for m = 0 .
05 along the (1 , q , is in the ( − , while that of the minority spin is compressed, along the(-1,1) direction so that the reconstructed surfaces become“nested” via k → − k + q as illustrated in the middle rowof Fig.3. The charge-density wave that develops at largerdensity is a consequence of the ( π, π ) nesting along the(1 , q , chargeorder must happen at ( π, π ) − q and one can dial the formof the order parameter interpolating between a purely su- n ! + n " n ! n " ! n " VCVI DIDC VCVI
FIG. 4:
U/t = 3 (left) and
U/t = 4 (right) T = 0 phase dia-gram. The color axis reports values of λ = mπ/ | q | . Trianglesrepresent the commensurate phases ( λ = 1). Up-triangleshave q ∝ (1 ,
1) (DC=Diagonal commensurate) while down-triangles have q ∝ (1 ,
0) (VC=Vertical Commensurate). Cir-cles and square are Vertical and Diagonal Incommensuratephases (VI and DI). perfluid and a purely charge-density wave state. Awayfrom half-filling, however, pairing is characterized by apolarization dependent wave-vector q while charge orderappears at ( π, π ) and cannot be dialed away. PHASE DIAGRAM
Our results are summarized in the phase diagram ofFig. 4 for U = 3 t and U = 4 t . We found no evidence ofdiagonal supersolid order at U = 3 t with a region of sta-bility for the commensurate phase limited to small polar-ization and density close to one. At U = 4 t , the diagonalwave vector is almost always commensurate with m , andpairing and charge order coexist in all but a small fractionof the phase diagram where the diagonal phase is the cor-rect ground state. This suggests that the charge-densitywave play a small but important role in stabilizing diag-onal order.Because we have accurately determined the pair mo-menta and the correct symmetry of the order param-eter as a function of magnetization and density, wecan now tackle the issue of phase separation in spin-imbalanced systems on a square lattice within Hartree-Fock-Bogoliubov theory. To do this we check the convex-ity of the energy by diagonalizing P = (cid:20) ∂µ∂n ∂µ∂m∂h∂n ∂h∂m (cid:21) (6)for the same set of densities and magnetizations reportedin Fig. 4. The values of µ , h and their derivatives aredetermined numerically from the Helmholtz free energyusing the grid in Fig.4. As a result, data in Fig. 5 havenumerical uncertainty that we estimated to be compa-rable to the symbol size. These data, summarizing thecase of U = 3 t , do not suggest any tendency towardphase separation. Although the smallest eigenvalue of P approaches 0 at small magnetization we do not interpretthis as a signal of incipient phase separation. Because thedistance between nodal lines is inversely proportional tothe magnetization, the ground state at small magnetiza-tions is characterized by the presence of essentially non-interacting domain walls where the excess spin particlesreside. The ground state energy is then simply deter-mined by the energy to create one such wall multipliedby the walls density. Because the latter is proportional tothe magnetization, the energy displays linear behavior in m and causes ∂h/∂m to become vanishingly small. This,in turn, is responsible for the vanishing behavior of oneof the eigenvalues. These results show that, contrary toprevious conclusion based on restricted searches[14] ontwo-dimensional lattices and contrary to the widely ac-cepted scenario in the dilute continuum, there is no phaseseparation in the true Hartree-Fock-Bogoliubov groundstate on a lattice in these parameter regimes.Similar results hold for the phases at U = 4 t when thedifferent order parameter patterns are separately consid-ered. Obviously, in the global phase diagram, the dif-ferent symmetry between diagonal and vertical phasesimplies that the transition is discontinuous and accom-panied by phase separation. We have not attempted anin-depth study of how this affects the situation at U = 4 t .However, a simpler analysis based on considering phaseseparation in two phases having the same density and dif-ferent magnetization leads to a narrow coexistence region(∆ m = 0 .
02 at n = 1) and makes it sensible to assumethat the broad feature of the phase diagram are, in fact,robust. SUMMARY
We have determined the type of phases that an imbal-anced population of two fermionic species with attrac-tive interactions support in the two-dimensional squarelattice. Their real- and momentum-space properties arequantitatively characterized. We find that unidirectionalLO states are the most stable mean field solutions ina large range of parameter regimes, and checkerboardorder is only clearly favored at small density and largepolarization. At lower U , the finite-momentum pairingwave-vector is along (1 , ,
1) direction. Our Z ee m a n fi e l d ( h ) n ↑ − n ↓ s m a ll e s t e i g e n v a l u e o f P n ↑ − n ↓ n=0.70n=0.72n=0.76n=0.80n=0.84n=1.00 FIG. 5: Left: Zeeman field, h , as a function of magnetization, m , for different densities. At small magnetizations the Zee-man field reaches a plateau as explained in the main body ofthe paper. Right: Smallest eigenvalue of the matrix P . Inthe density/magnetization regime considered this eigenvalueremains positive (although very small in correspondence ofthe Zeeman field plateau) indicating no region of phase sepa-ration. Data are for U = − t . results suggest that, besides time of flight or Bragg spec-troscopy, the different local phases we discussed are alsoidentifiable by their local density profiles. This, in turn,should help their real-space characterization in the pres-ence of a confining potential. Our results conclusivelyshow that, within Hartree-Fock-Bogoliubov, there is nophase separation in the presence of a lattice for the con-sidered density and magnetization regimes. Acknowledgments
We acknowledge support from NSF (Grant no.DMR-1006217), DOE (de-sc0008627), and ARO (Grantno. 56693-PH) and computational support from the Cen-ter for Piezoelectrics by Design and by DOE leadershipcomputing through an INCITE grant. [1] P. Fulde and R. A. Ferrell, Physical Review , A550(1964), ISSN 0031-899X, URL http://link.aps.org/doi/10.1103/PhysRev.135.A550 .[2] A. Larkin and I. Ovchinnikov, Soviet Physics JETP , 762 (1965), URL .[3] D. E. Sheehy and L. Radzihovsky, Physical Review Let-ters , 060401 (2006), URL http://arxiv.org/abs/cond-mat/0508430 .[4] A. Moreo and D. Scalapino, Physical Review Letters ,216402 (2007), ISSN 0031-9007, URL http://link.aps.org/doi/10.1103/PhysRevLett.98.216402 . [5] M. J. Wolak, B. Gr´emaud, R. T. Scalettar, and G. G.Batrouni, Physical Review A p. 15 (2012), 1206.5050,URL http://arxiv.org/abs/1206.5050 .[6] C.-C. Chang and S. Zhang, Phys. Rev. Lett. ,116402 (2010), URL http://link.aps.org/doi/10.1103/PhysRevLett.104.116402 .[7] T. Koponen, J. Kinnunen, J.-P. Martikainen,L. M. Jensen, and P. T¨orm¨a, New Journal ofPhysics , 179 (2006), ISSN 1367-2630, URL http://stacks.iop.org/1367-2630/8/i=9/a=179?key=crossref.58b7381cef82e44d246df7a855d2b7dd .[8] T. K. Koponen, T. Paananen, J.-P. Martikainen,and P. T¨orm¨a, Physical Review Letters , 120403(2007), URL http://link.aps.org/doi/10.1103/PhysRevLett.99.120403 .[9] Y. L. Loh and N. Trivedi, arXiv09070679v1 cond-matquantgas , 4 (2009), URL http://arxiv.org/abs/0907.0679 . [10] Q. Wang, H.-Y. Chen, C.-R. Hu, and C. S. Ting, Phys.Rev. Lett. , 117006 (2006), URL http://link.aps.org/doi/10.1103/PhysRevLett.96.117006 .[11] J. Xu, C.-C. Chang, E. J. Walter, and S. Zhang, Journalof physics. Condensed matter : an Institute of Physicsjournal , 505601 (2011), ISSN 1361-648X, URL .[12] H. Shimahara, Journal of the Physical Society of Japan , 736 (1998), URL http://arxiv.org/abs/cond-mat/9711017 .[13] M. M. Parish, S. K. Baur, E. J. Mueller, and D. A. Huse,Phys. Rev. Lett. , 250403 (2007), URL http://link.aps.org/doi/10.1103/PhysRevLett.99.250403 .[14] P. A. Igoshev, M. A. Timirgazin, A. A. Katanin,A. K. Arzhnikov, and V. Y. Irkhin, Phys. Rev. B , 094407 (2010), URL http://link.aps.org/doi/10.1103/PhysRevB.81.094407http://link.aps.org/doi/10.1103/PhysRevB.81.094407