Ising t-J model close to half filling: A Monte Carlo study
aa r X i v : . [ c ond - m a t . s t r- e l ] O c t Ising t – J model close to half filling: A Monte Carlostudy M M Ma´ska , M Mierzejewski , A Ferraz , E A Kochetov Department of Theoretical Physics, Institute of Physics, University of Silesia,40–007 Katowice, Poland International Center of Condensed Matter Physics, Universidade de Brasilia, CaixaPostal 04667, 70910–900 Brasilia, DF, Brazil Laboratory of Theoretical Physics, Joint Institute for Nuclear Research, 141980Dubna, Russia
Abstract.
Within the recently proposed doped-carrier representation of theprojected lattice electron operators we derive a full Ising version of the t – J model.This model possesses the global discrete Z symmetry as a maximal spin symmetry ofthe Hamiltonian at any values of the coupling constants, t and J . In contrast, in thespin anisotropic limit of the t – J model, usually referred to as the t – J z model, the global SU (2) invariance is fully restored at J z = 0, so that only the spin-spin interaction hasin that model the true Ising form. We discuss a relationship between those two modelsand the standard isotropic t – J model. We show that the low-energy quasiparticlesin all three models share the qualitatively similar properties at low doping and smallvalues of J/t . The main advantage of the proposed Ising t – J model over the t – J z oneis that the former allows for the unbiased Monte Carlo calculations on large clustersof up to 10 sites. Within this model we discuss in detail the destruction of theantiferromagnetic (AF) order by doping as well as the interplay between the AF orderand hole mobility. We also discuss the effect of the exchange interaction and that of thenext nearest neighbour hoppings on the destruction of the AF order at finite doping.We show that the short-range AF order is observed in a wide range of temperaturesand dopings, much beyond the boundaries of the AF phase. We explicitly demonstratethat the local no double occupancy constraint plays the dominant role in destroyingthe magnetic order at finite doping. Finally, a role of inhomogeneities is discussed.PACS numbers: 71.10.Fd, 71.10.Pm, 74.25.Jb sing t – J model close to half filling: A Monte Carlo study
1. Introduction
Since the discovery of the high–temperature superconductors (HTSC), many theoreticalinvestigations have been focused on the study of the doping evolution from theantiferromagnetically ordered Mott insulator to the BCS-type superconductor. It iscommonly assumed that this evolution is the key ingredient for the understanding of thephysics behind high–temperature superconductivity and it can adequately be addressedin terms of the two–dimensional (2 d ) t - J lattice model. This model is believed to capturethe essential low-energy physics of doped Mott insulators driven by strong electroncorrelations. Although, it is generally accepted that strong electron correlations playan important role close to half filling, it remains unclear to what extent HTSC can bedescribed entirely in terms of this minimal electronic model.One of the reasons for this lack of clarity is that, despite its simplicity, the exactproperties of the t − J model, apart from a few limiting cases, are still unknown. Thevast majority of the results away from half filling have been obtained for one or two holesintroduced into the AF background. This problem has been thoroughly analyzed withthe help of various analytical and numerical approaches [1, 2, 3, 5, 4, 6, 7, 8, 9, 10, 11].Results for larger dopings are less comprehensive. Here, one of the important problemsconcerns the robustness of the AF order against doping. Most of the theoreticalapproaches predict, that the long range AF order persists up to much larger dopingsthan observed in cuprates [12, 13]. It has recently been suggested that including intoconsideration the nearest neighbour hopping may help in resolving this discrepancy [14].However, a reliable and well controlled analytical treatment of the isotropic t – J modelposes a severe technical problem: it is very hard to analytically deal with the localno double occupancy (NDO) constraint. On the other hand, because of this constraint,numerical treatment is available only on rather small lattice clusters, which immediatelyrises the problem of the finite-size effects and, consequently, of the thermodynamicalstability of the results obtained.An interesting question then arises: is there available a modification of the isotropic t − J model that allows for unbiased numerical treatment and at the same time capturesthe essential properties of the original t − J model at least for a certain parameter range?Some of the computational difficulties related with the t – J model may be overcome bymeans of investigations of its anisotropic limit, i.e., the t – J z model. Numerical workhas suggested that the t − J and t − J z models have many similar properties [15]. Inparticular, the stripes, pairing, and Nagaoka states found in the t − J z model are verysimilar to those of the isotropic t − J model [16, 17, 18, 19, 20]. In general, despite somesignificant differences between both the models [20, 21], at the low-energy scale of order J ( ≪ t ) it is reasonable to consider the background spin configuration to be frozen withrespect to the hole dynamics time scale. In this case the properties of the low-energyquasiparticle excitations in the t − J model are at least qualitatively similar to thosein the anisotropic t − J z model. Although this model is more amenable to numericalcalculations, again only rather small lattice clusters are allowed. sing t – J model close to half filling: A Monte Carlo study t − J model in which the t -termpossesses the global discrete Z spin symmetry rather than the SU(2) one results in amore tractable though still nontrivial model. However, it is not clear how such a modelcan be derived directly in terms of the Gutzwiller projected lattice electron operators,since those operators transform themselves in a fundamental representation of SU(2). Inthis paper we use the recently proposed doped-particle representation of the projectedelectron operators to derive the full Ising version of the t − J model. We believe that theresults of our study capture some essential low-energy properties of the isotropic t − J model since the strongly correlated nature of the problem is preserved.The doped-particle representation of the t - J model is especially suited forinvestigations of the underdoped regime [22, 23, 24]. Some crucial points of this approachare recalled at the beginning of the following section. Although this is a slave–particleformulation, it differs from other similar approaches in that the NDO constraint isidentically fulfilled at half–filling. Since the holes are the only charge degrees of freedompresent in the system, the number of charge carriers is very small close to half-filling.Therefore, this approach is particularly useful for the description of the underdopedregime, whereas its applicability to the description of strongly overdoped cuprates ismuch more involved.Until now, the doped particle representation has been analyzed only within mean–field approximations [22, 23]. Our aim is to go beyond this approximation in such away, that the slave–particle constraint is exactly taken into account. We demonstratethat this can be achieved for very large clusters of the order of 10 lattice sites, providedthe SU(2) symmetry is broken down to Z not only in the spin–spin interaction termbut also in the hopping term. Since the effectiveness of the resulting calculations isindependent of the translational invariance, this approach also allows one to investigatethe role of inhomogeneities which are expected to play an important role in the cupratecompounds [25].The paper is organized as follows. Section 2 describes in detail the derivation of thetheoretical model. Section 3 comprises the results of the numerical calculations. Section4 lists our conclusions.
2. Model and Approach t – J model We start with the t – J Hamiltonian on a square lattice [26] H tJ = − X ijσ t ij ˜ c † iσ ˜ c jσ + J X h ij i ( Q i Q j −
14 ˜ n i ˜ n j ) , (1)where ˜ c iσ = P c iσ P = c iσ (1 − n i, − σ ) is the projected electron operator (to exclude theon-site double occupancy), Q i = P σ,σ ′ ˜ c † iσ τ σσ ′ ˜ c iσ ′ , τ = 3 / , is the electron spin operatorand ˜ n i = P n i P = n i ↑ + n i ↓ − n i ↑ n i ↓ . Hamiltonian (1) contains a kinetic term with thehopping integrals t ij and a potential J describing the strength of the nearest neighbour sing t – J model close to half filling: A Monte Carlo study P = Q i (1 − n iσ n i − σ ) projects out the doubly occupied states | ↑↓i thereby reducingthe quantum Hilbert space to a product of 3-dimensional spaces spanned by the states | i i , | ↑i i and | ↓i i . Physically this modification of the original Hilbert space results instrong electron correlation effects. The crucial local no double occupancy constraintis rigorously incorporated into Eq. (1). However, this is achieved at the expense ofintroducing the constrained electron operators, ˜ c † iσ , that obey much more complicatedcommutation relations than the conventional ”unconstrained” fermion operators. Itshould be stressed that it is precisely close to half filling where the Gutzwiller projectionis of a crucial importance: the projected electron operator ˜ c † iσ in this regime significantlydiffers from the bare electron operator c † iσ (right at half filling ˜ c † iσ = 0).A natural question then arises: is it possible to rewrite the t - J Hamiltonian in termsof the conventional fermion and spin operators in such a way that the NDO constraintfor the lattice electrons transforms itself into the one that can be treated, close to halffilling, in a controlled way? Recently, it has been shown that the t – J Hamiltonian canindeed be represented in that form [23, 22, 24].For the reader’s convenience we sketch below the main points of this scheme. Thebasic idea behind this approach is to assign fermion operators to doped carriers (holes,for example) rather than to the lattice electrons. The t - J Hamiltonian is expressed thenin terms of the lattice spin operators, S i , and doped carrier operators represented byspin-1 / d iσ .To accommodate these new operators one obviously needs to enlarge the originalonsite Hilbert space of quantum states. This enlarged space is characterized by the statevectors | σa i with σ = ⇑ , ⇓ labeling the spin projection of the lattice spins and a = 0 , ↑ , ↓ labeling the dopon states (double occupancy is not allowed). In this way the enlargedHilbert space becomes H enli = {| ⇑ i i , | ⇓ i i , | ⇑↓i i , | ⇓↑i i , | ⇑↑i i , | ⇓↓i i } , (2)while in the original Hilbert space we can either have one electron with spin σ or avacancy: H = {| ↑i i , | ↓i i , | i i } . (3)The following mapping between the two spaces is then defined [22]: | ↑i i ↔ | ⇑ i i , | ↓i i ↔ | ⇓ i i , (4) | i i ↔ | ⇑↓i i − | ⇓↑i i √ . (5)The remaining states in the enlarged Hilbert space, ( | ⇑↓i i + | ⇓↑i i ) / √ | ⇑↑i i , | ⇓↓i i are unphysical and should therefore be removed in actual calculations. In this mapping,a vacancy in the electronic system corresponds to a singlet pair of a lattice spin and adopon, whereas the presence of an electron is related to the absence of a dopon. sing t – J model close to half filling: A Monte Carlo study t - J Hamiltonian then reads [24] H t − J = X ijσ t ij ˜ d † iσ ˜ d jσ + J X h ij i [( S i + M i ) ( S j + M j ) − (cid:16) − ˜ n di (cid:17) (cid:16) − ˜ n dj (cid:17)(cid:21) , (6)with ˜ d iσ = d iσ (1 − d † i, − σ d i, − σ ) being a projected dopon operator and ˜ n di = P σ ˜ d † iσ ˜ d iσ .The application of H t − J in this form should be accompanied by the implementationof the constraint to eliminate the unphysical states, S i M i + 34 ˜ n di = 0 , (7)where M i = P σ,σ ′ ˜ d † iσ τ σσ ′ ˜ d iσ ′ stands for the dopon spin operator so that Q i = S i + M i . Note the important factor of 2 in front of the first term in Eq. (6). It originates fromthe fact that the vacancies are represented in this theory by the spin-dopon singletsgiven by Eq. (5). The projected lattice electron operators can be explicitly expressedin terms of the projected dopon operators. For example,˜ c † i ↑ = √ P phi ˜ d i ↓ P phi = 1 √ (cid:20)(cid:18)
12 + S zi (cid:19) ˜ d i ↓ − S + i ˜ d i ↑ (cid:21) , (8)where P phi = 1 − ( S i M i + ˜ n di ) is the projection operator which eliminates the unphysicalstates from the i th site.The above representation of the t – J Hamiltonian is particularly useful for thedescription of strongly underdoped cuprates. Close to half filling n di ≪
1, so that onecan safely drop the tilde sign off the projected dopon operators. This is due to the factthat in the low doping regime the probability for the realization of a doubly occupieddopon state is indeed very low. Despite that, the NDO constraint (7) must be imposedto eliminate the unphysical degrees of freedom that are present in this formalism atany finite doping. Note, however, that at half–filling the left hand side of Eq. (7)vanishes, and thus, in contrast to the original NDO constraint for the lattice electrons,this equation turns into a trivial identity. ‡ Additionally, in this regime one can neglectboth the hole–hole interactions represented by the M i M j term, and the ˜ n di ˜ n dj couplings.Note also that the insulating phase in this representation is directly associated with theabsence of charged particles. t – J z and the Ising t − J models So far, the doped–particle representation of the t - J model has been analyzed only withinthe mean–field approximations [23, 22, 24]. Our aim is to go beyond the mean–fieldanalysis and to carry out calculations for large enough systems, to make sure that thefinite–size effects are truly negligible. Let us start with the anisotropic limiting case ofthe t – J model with the spins polarized only along the z-component, namely, the t – J z ‡ The original local NDO constraint for the lattice electrons, P σ c † iσ c iσ ≤ , right at half filling reads P σ c † iσ c iσ = 1. sing t – J model close to half filling: A Monte Carlo study t – J z model can be considered as a limiting case of the t – J model (1) whichhas an Ising rather than a Heisenberg spin interaction: H t − J z = − X ijσ t ij ˜ c † iσ ˜ c jσ + J z X h ij i (cid:18) Q zi Q zj −
14 ˜ n i ˜ n j (cid:19) , (9)Here ˜ c iσ represent the Gutzwiller-projected electron operators. The global continuousspin SU(2) symmetry of the t – J model now reduces to the global discrete Z symmetryof the t – J z model. Although Q zi Q zj interaction possesses discrete Z symmetry, theoriginal SU(2) symmetry of all other terms of the Hamiltonian is preserved. Therefore,the symmetry of the t – J z model depends on whether J z is zero or finite. Namely, for J z = 0 the SU(2) symmetry is restored again. In contrast, in the full t − J model boththe t - and J -terms possess the same SU(2) symmetry.It is therefore natural to seek for a representation of the full Ising version of the t − J model in which the symmetry of the model does not depend on the values ofthe model parameters. Such a representation can straightforwardly be derived withinthe dopant–particle formulation of the t – J model. The physical consequences as wellas computational advantages of such an approach will be discussed in the subsequentsections.To proceed, we start right from the original t – J Hamiltonian described in terms ofthe lattice electrons given by Eq. (1), where we now put Q + i = Q − i = 0 at the operatorlevel. We then have for the physical electron projected operators:˜ c i ↑ = P phi ˜ d † i ↓ P phi = ( 12 + S zi ) ˜ d † i ↓ , ˜ c i ↓ = P phi ˜ d † i ↑ P phi = ( 12 − S zi ) ˜ d † i ↑ , (10)where the projection operator now reads P phi = 1 − (2 S zi M zi + ˜ n di / Q + i = ( Q − i ) † = ˜ c † i ↑ ˜ c i ↓ ≡ . The kinetic t -term built out of the physical electron operators given by Eqs. (10)possesses the global Z symmetry rather than the SU(2) one.Accordingly, the underlying onsite Hilbert space rearranges itself in the followingway. The operators ˜ c i ↓ , ˜ c † i ↓ act on the Hilbert space H ↓ = {| ⇓ , i , | ⇓ , ↑i} . Theseoperators do not mix up any other states. Operator ˜ c i ↓ destroys the spin-down electronand creates a vacancy. This vacancy is described by the state | ⇓ , ↑i . The similarconsideration holds for the ˜ c i ↑ operators. Now, however, the vacancy is described bythe state | ⇑ , ↓i . Those two vacancy states are related by the Z transformation. Theoperator ( Q zi ) = (1 − ˜ n di ) produces zero upon acting on the both. The physicalHilbert state is therefore a direct sum H ph = H ↑ ⊕ H ↓ . Under the Z transformation( ↑↔↓ , S zi → − S zi ) we get H ↑ ↔ H ↓ , which results in H ph → H ph . § § In the isotropic t – J model these two 2 d spaces merge into a 3 d SU(2) invariant physical space, where sing t – J model close to half filling: A Monte Carlo study Z symmetry whereas the global continuous SU(2) symmetry iscompletely lost. Close to half–filling this Hamiltonian reduces to the form, H Isingt − J ≡ H t − J | Z = X ijσ t ij ˜ d † iσ ˜ d jσ + J X h ij i (cid:20)(cid:18) S zi S zj − (cid:19) + S zi M zj + S zj M zi i , (11)which should be accompanied by the constraint,2 S zi M zi + 12 ˜ n di = 0 . (12)In order to emphasize the difference between the t – J z and the present model we dubthe latter the Ising t − J or for short the t − J | Z model.The factor of 2 which is presented in the hopping term of the isotropic t – J Hamiltonian (6) drops out from the hopping term in the Ising t − J Hamiltonian. Thisoccurs because of the fact that the equations˜ c † i ↑ = √ P phi ˜ d i ↓ P phi , P phi = 1 − (cid:18) S i M i + 34 ˜ n di (cid:19) which are valid for both the t – J and t – J z models, are in the t − J | Z model replaced bythe following ones, ˜ c † i ↑ = P phi ˜ d i ↓ P phi , P phi = 1 − (cid:18) S zi M zi + 12 ˜ n di (cid:19) . In practical calculations, the NDO constraint (12) can be taken into account withthe help of a Lagrange multiplier. In order to do that we introduce an additional termto the Hamiltonian, λ X i (cid:18) S zi M zi + 12 n di (cid:19) = λ X i (cid:20)(cid:18)
12 + S zi (cid:19) d † i ↑ d i ↑ + (cid:18) − S zi (cid:19) d † i ↓ d i ↓ (cid:21) . (13)Notice that the operator ( + S zi ) d † i ↑ d i ↑ + ( − S zi ) d † i ↓ d i ↓ produces eigenvalues 0 and 1,when acting on the onsite physical and unphysical states, respectively. Because of this,the global Lagrange multiplier λ → ∞ enforces the NDO constraint locally. The doubledopon occupancy of an arbitrary site results in an appearance of an unphysical stateand hence enhances the total energy by λ . Therefore, in the large- λ limit all unphysicalstates are automatically eliminated and we can in this limit safely remove the tilde signoff the d operators. In the following section we show that this constraint is of crucialimportance for the description of the AF order at finite doping. the vacancy is just an antisymmetric linear combination given by the SU(2) spin singlet (5). Thesymmetric combination splits off, since it represents an unphysical spin-triplet state. sing t – J model close to half filling: A Monte Carlo study The total Hamiltonian takes the form H λt − J | Z = H ↑ + H ↓ + J X h ij i S zi S zj + const , (14)with H ↑ = X ij t ij d † i ↑ d j ↑ + X i d † i ↑ d i ↑ λ (cid:18)
12 + S zi (cid:19) + J X h j i i S zj , (15) H ↓ = X ij t ij d † i ↓ d j ↓ + X i d † i ↓ d i ↓ λ (cid:18) − S zi (cid:19) − J X h j i i S zj , (16)where h j i i denotes neighbouring sites of a given site i . We have neglected the hole–holeinteraction in H λt − J | Z , which is perfectly justified in the low doping regime. In order toverify this approximation we have carried out additional calculations with the hole-holeinteraction being taken into account in the mean–field approximation. The differenceis negligible and therefore we do not present them here. Note, that the interactionstrength in H λt − J | Z is exactly the same as in the standard formulation of the t – J model.Absence of any renormalization of the model parameters originates from the fact thatthe projection procedure is explicitly built in Eq. (14), provided λ → ∞ . Despite itscomplexity, with the Monte Carlo (MC) method, one can investigate the Hamiltonian(14) for very large systems without any approximation. Since [ S zi , H λt − J | Z ] = 0 the spindegrees of freedom can be analyzed within the classical Metropolis algorithm. However,since the effective Hamiltonian (14) includes both fermionic as well as classical degreesof freedom, this algorithm needs to be modified. The procedure is as follows:(i) an initial configuration of { S zi } is generated;(ii) the Hamiltonians (15, 16) are diagonalized and the free energy F of the fermionicsubsystem in the canonical ensemble is determined;(iii) two sites with opposite spins S z are randomly chosen; then, both the spins areflipped;(iv) step (ii) is repeated, determining new value of the free energy F ′ ;(v) if F ′ < F or exp [( F − F ′ ) /kT ] > x , where x is a random number from the interval[0; 1), the new { S zi } configuration is accepted, added to the ensemble and theprocedure goes to step (iii), otherwise it goes directly to step (iii).It is the Metropolis algorithm, but with the internal energy in statistical weightsreplaced by the free energy of the fermionic subsystem. A detailed description of thisapproach can be found in Ref. [27]. Concurrently with the MC simulation an iterateprocedure calculating the distribution of holes is carried out in a self–consistent way. sing t – J model close to half filling: A Monte Carlo study ×
20 systems with periodicboundary conditions. However, in order to check the influence of finite–size effects wehave carried out calculations on clusters of up to 1600 lattice sites and with averagingover the boundary conditions [28]. This problem is discussed at the end of the nextsection.In order to eliminate the unphysical states in the Monte Carlo simulations, we havetaken λ = 100 t . Therefore, the Lagrange multiplier is by far the largest energy scale inthe system, which guaranties the single occupancy of each lattice site. The simulationshave been carried out in the canonical ensemble, what allows for accurate control ofthe doping level. We assume an absence of the ferromagnetic order. Namely, we take P i S zi = P i M zi = 0, what is reflected in the third point of the MC procedure.The aim of the simulations is to determine how the antiferromagnetic order andthe spectral properties are affected by doping. We start our discussion with the dopingdependence of the spin–spin correlation function for the projected physical electronoperators g ( r ) = 4 N X i X j h ( S zi + M zi )( S zj + M zj ) i× exp [ i K · ( R i − R j )] ¯ δ ( r − | R i − R j | ) , (17)where K = ( π, π ) and ¯ δ ( x ) = ( | x | ≤ . a, , with a being the lattice constant. h . . . i in Eq. (17) means an average over thespin configurations generated in MC run. g ( r ) allows one to distinguish betweenthe long range order (LRO), when it remains finite for arbitrary r , quasi long rangeorder (QLRO), when g ( r ) decays algebraically, and a short range order (SRO), when g ( r ) decays exponentially. Calculations of the spin–spin correlation function will beaccompanied by results obtained for a static spin–structure factor, defined as S ( q ) = 1 N X ij e i q ( R i − R j ) h ( S zi + M zi )( S zj + M zj ) i . (18)The third quantity that we use in the following discussion, is the hole spectral functiongiven by A ( k , ω ) = − π Im G (cid:16) k , ω + i + (cid:17) , (19)where G ( k , z ) = X i X j exp { i k ( R i − R j ) }× hG σ ( R i , R j , z ) (cid:20) − s ( σ ) S zi (cid:21) (cid:20) − s ( σ ) S zj (cid:21) i , (20)with s ( ↑ ) = 1 and s ( ↓ ) = −
1. Here, similarly to Eq. (17), h . . . i indicates averaging overspin configurations and G σ ( R i , R j , z ) = n [ z − H σ ] − o ij (21) sing t – J model close to half filling: A Monte Carlo study { S zi } . The presenceof factors − s ( σ ) S zi in Eq. (20) follows from Eq. (10). Note that the spin–spincorrelation function, the spin–structure factor and the spectral function are defined forphysical electron operators, ˜ c i .
3. Numerical results
As discussed in the preceding sections the derived representation differs from thestandard t – J z Hamiltonian in that the SU(2) symmetry is broken also for J = 0. Inorder to visualize the physical consequences of this difference we start with calculationsfor the one–hole case. In this regime large clusters have been analyzed numerically bothfor the t – J z and t – J models. J/t -4-3-2-10 E /t present results t-J z model t-J model Figure 1.
One–hole energy as a function of
J/t for a 8 × T = 0. Thedata labeled as t − J z and t − J have been taken from Refs. [29] and [11], respectively.In the latter case, it is the energy of a hole with momentum ( π/ , π/ × t − J | Z and t – J z systems. Here, the dashed line showsresults presented in Ref. [30] for the t – J z model. Since for J → P i S zi = 0 is now relaxed. In the main panel of Fig. 1 we show the one–hole energy calculated at T = 0 fora 8 × ˜ n i ˜ n j term. Here, we compare our data with exact resultsobtained for a 50–site t – J z cluster [29] as well as with recent exact results for a bulk t – J system [11]. In the inset of Fig. 1 we compare exact results obtained for 4 × t − J | Z and t – J z clusters [30] for a wider range of the J/t ratio. In two limiting cases
J/t → J/t → ∞ the one–hole energy obtained in our approach is the same as in the t – J z model. It can be explained in the following way. For J = 0, the ferromagnetic Nagaokastate becomes a ground state in both the approaches. In this case the propagation of ahole is not perturbed by the magnetic order and the one–hole energy equals − t . Themain difference between t – J z and t − J | Z approaches consists in the symmetry of the sing t – J model close to half filling: A Monte Carlo study J/t ≫
1, when the system propertiesare determined predominantly by the same spin–spin interaction. In the regime ofintermediate J the differences are most pronounced. The one–hole energy obtained forthe t – J model in this regime is in between the results obtained for t – J z and t − J | Z approaches. Here, the differences between t – J and t – J z models are comparable to thosebetween t – J and t − J | Z ones. δ S ( π , π ) J = 0.2 tJ = 0.4 t C) g ( r ) r g ( r ) δ = 0.02 δ = 0.04 δ = 0.06 δ = 0.08 δ = 0.02 δ = 0.04 δ = 0.06 δ = 0.08 A)B)
Figure 2.
Panels A) and B) show g ( r ) calculated for kT = 0 . t . J = 0 . t (A)and J = 0 . t (B). The curves from the top to the bottom have been obtained for δ = 0 . , . , . , .
08. Solid (dashed) lines show results obtained for t ′ = t ′′ = 0( t ′ = − . t and t ′′ = 0 . t ). Panel C) shows doping dependence of the static spin–structure factor S ( π, π ) at kT = 0 . t for t ′ = − . t and t ′′ = 0 . t . Next, we investigate how the doping affects the antiferromagnetic order. Theprevious studies of the t – J model clearly indicate that the long range hopping amplitudessignificantly modify the bandwidth and the dispersion of the quasiparticles [31].Recent Green’s function Monte Carlo calculations demonstrate that the next nearest sing t – J model close to half filling: A Monte Carlo study t – J model withonly the nearest neighbour hopping, antiferromagnetic correlations persist up to holeconcentrations much larger than the ones observed in HTSC materials. One may expectthat in the absence of the transverse spin–spin interaction the robustness of LRO shouldbe even more pronounced. Moreover, the previous analysis of the one– and two–holespectra in the t – J z model [20] has shown that for t/J <
5, half of the one–hole bandwidth does not exceed 0 . t . Therefore, the intra–sublattice hopping should be a sourceof an important contribution to the kinetic energy even for small values of t ′ and t ′′ .The significance of the long range hopping for the AF order is demonstrated in Fig.2, where we compare the spin–spin correlation functions calculated with and without t ′ and t ′′ . In all the figures showing g ( r ) we use logarithmic scale for the vertical axis.Therefore, for LRO, QLRO and SRO, g ( r ) should be represented asymptotically by aconstant function, logarithmic function and a straight line, respectively. In the following δ denotes the average concentration of holes.One can see that the influence of the long range hopping depends, even qualitatively,on the value of the exchange coupling J . For a small value of J , the AF order is enhanced when hoppings to second and third nearest neighbours are allowed (for J = 0 . t see panelA) in Fig. 2). On the other hand, for bigger values of J , these hoppings reduce the AForder (for J = 0 . t see panel B) in Fig. 2). Such a behaviour could be explained asfollows. In the presence of only nearest neighbour hopping there is a strong competitionbetween the energy of spin–spin interaction and the hole kinetic energy. It results fromthe fact, that in this case only inter sublattice hopping is allowed. From Eq. (13)one can then infer that it is possible only in regions where the AF order is absent.Then, nonzero t ′ and t ′′ allow for intra sublattice hopping, thereby leading to gainingof the kinetic energy without destroying the AF order. This mechanism is effective for J ≤ . t . On the other hand, for t ′ = t ′′ = 0 and large J , holes are almost localizedand, therefore, only weakly frustrate the AF state. The intra–sublattice hopping allowsfor the propagation of holes, which effectively reduces the AF LRO.The doping induced destruction of the AF LRO can directly be seen in Fig.2C, where we present spin–structure factor obtained for t ′ = − . t and t ′′ = 0 . t .This quantity is important in that it is directly accessible in, e.g., neutron scatteringexperiments. The maximal doping for which the AF state still exists strongly dependson the magnitude of the exchange interaction. This result contrasts with the recentlyreported Green’s function Monte Carlo study of the t − J model [14], where the AFLRO vanishes at δ = 0 . δ = 0 .
13 for J = 0 . t and J = 0 . t , respectively. Inour approach the experimental data for the critical doping in HTSC can be reproducedprovided J < . t .In order to illustrate the interplay between the AF order and the mobility of holeswe have calculated the hole spectral functions A ( k , ω ) (see Figs. 3 and 4). For t ′ = t ′′ = 0 and small δ one can see almost localized particles with very small dispersion.Similar situation occurs in the t – J z model, but it is not the case for the t − J one, where sing t – J model close to half filling: A Monte Carlo study -6-4-2024 / ( , ) (0, ) (0,0) ( , )-6-4-2024 / ( , ) (0, ) (0,0) ( , )-6-4-2024 / ( , ) (0, ) (0,0) ( , ) -6-4-2024 / ( , ) (0, ) (0,0) ( , )-6-4-2024 / ( , ) (0, ) (0,0) ( , )-6-4-2024 / ( , ) (0, ) (0,0) ( , ) J=0.2t J=0.4t δ = . δ = . δ = . Figure 3.
Spectral functions A ( k , ω ) calculated for kT = 0 . t along the mainsymmetry lines of the Brillouin zone with t ′ = t ′′ = 0. the the spin–flip term can undo the defects generated by a moving hole and hence itallows for its much higher mobility. In the present approach holes become mobile whendoping increases, i.e., when the AF background disappears. The onset of mobile holesis accomplished through a gradual transfer of the spectral weight from the vicinity ofthe almost localized level. Close to half–filling the most significant transfer takes placein states with k = (0 ,
0) and k = ( π, π ). However, even for relatively large doping thespectral functions remain broad for all the momenta.For non–zero t ′ and t ′′ there are mobile holes even for small doping, but the spectralfunctions still remain very broad. Also in this case, doping is responsible for significantmodification of the dispersion relation of holes. In Fig. 4 we compare A ( k , ω ) calculatedfor δ = 0 . , .., .
14 with t ′ = − . t and t ′′ = 0 . t . Along with the destruction ofthe AF LRO, there is an increasing contribution of nearest neighbour hopping to thehole kinetic energy. For δ = 0 .
02 the peaks in spectral functions can be fitted bythe dispersion relation with t = 0, whereas for δ = 0 .
14 the AF correlations hardlyinfluence the nearest neighbour hopping. Note that such a substantial modification ofthe dispersion relation may change the topology of the Fermi surface. Doping affectsnot only the effective dispersion relation, but also frustrates the AF background. The sing t – J model close to half filling: A Monte Carlo study -6-4-2024 / ( , ) (0, ) (0,0) ( , ) -6-4-2024 / ( , ) (0, ) (0,0) ( , )-6-4-2024 / ( , ) (0, ) (0,0) ( , ) -6-4-2024 / ( , ) (0, ) (0,0) ( , ) δ = 0.02 δ = 0.06δ = 0.14δ = 0.10 Figure 4.
Spectral functions A ( k , ω ) calculated along the main symmetry lines of theBrillouin zone for J = 0 . t , kT = 0 . t and various dopings. t ′ = − . t and t ′′ = 0 . t have been assumed. latter effect is responsible for strong broadening of the spectral functions that is visiblein Fig. 4. Comparison of Fig. 2 with Figs. 3 and 4 demonstrates that the mobility ofholes and destruction of the AF LRO are mutually connected with each other in thesense that mobility affects AF LRO and, vice versa , AF order affects the mobility ofholes.We now turn our attention to the temperature dependence of the spin–spincorrelation function and the spectral properties of holes. It is known that theN´eel temperature drops rapidly when the parent compounds of the high–temperaturesuperconductors are doped with holes. Similar behaviour can be inferred from Fig. 5,where temperature dependence of g ( r ) is presented for different doping levels. One cannote that doping strongly reduces the LRO, whereas its influence on the SRO is muchweaker. In particular, the nearest neighbour correlation functions g (1) calculated for δ = 0 .
02 and δ = 0 .
06 are qualitatively and quantitatively close to each other. Theseresults suggest that, the AF SRO should be observed in a wide range of temperaturesand dopings, much beyond the boundaries of the AF phase. It remains in agreement withrecent experiments on high–temperature superconductors suggesting that with doping,the long-range N´eel order gives way to short-range order with a progressively shortercorrelation length. As a result, at optimal doping the static spin correlation length isno more than two or three lattice spacing [32].The discussed above correlation function and the spin–structure factor describethe background composed of the localized spins, which, as mentioned in the precedingparagraphs, is up to some degree affected by the motion of doped holes. Therefore, the sing t – J model close to half filling: A Monte Carlo study g ( r ) r = 1 r = 2 r = 10 k B T/t g ( r ) r = 1 r = 2 r = 10 k B T/t S ( π , π ) k B T/t S ( π , π ) Figure 5. g ( r ) as a function of temperature for J = 0 . t and δ = 0 .
02 (upper panel)and δ = 0 .
06 (lower panel). The lines from the top to the bottom show g (1), g (2) and g (10). t ′ = − . t and t ′′ = 0 . t have been assumed. reduction of the antiferromagnetic correlations has to be observed also in the dynamicsof the carriers. It is shown in Fig. 6, where we demonstrate the temperature dependenceof the spectral functions. When the temperature increases, the number of spin defectsin the N´eel state increases as well, and this enables the nearest neighbour hopping,thereby allowing holes to lower their kinetic energy. This mechanism leads to trappingof holes in the regions of broken antiferromagnetic bonds and forming ferromagneticspin polarons, where the hole hopping does not frustrate the spin background. Thecontribution of the nearest neighbour hopping becomes visible in the spectral functions,where the increase of the temperature causes a significant broadening of the spectrallines. Similarly to the spectral functions obtained for a single hole in the t – t ′ – t ′′ – J model[33], the width of the peaks in the A ( k , ω ) is too small, when compared to the resultsof the angle–resolved spectroscopy (ARPES) measurements [34] on Sr CuO Cl . It hasrecently been argued, that strong electron–phonon interaction [35, 36] may explain thevery broad peaks observed in the insulating copper oxides [37, 38].In the present approach, holes interact with spins through the intersite interactionof strength J as well as through the onsite constraint with the Lagrange multiplier λ , sing t – J model close to half filling: A Monte Carlo study -6-4-2024 / ( , ) (0, ) (0,0) ( , ) -6-4-2024 / ( , ) (0, ) (0,0) ( , )-6-4-2024 / ( , ) (0, ) (0,0) ( , ) -6-4-2024 / ( , ) (0, ) (0,0) ( , )-6-4-2024 / ( , ) (0, ) (0,0) ( , ) -6-4-2024 / ( , ) (0, ) (0,0) ( , ) δ = 0.02 k T = . t δ = 0.06 k T = . t k T = . t Figure 6.
Spectral functions A ( k , ω ) calculated along the main symmetry lines of theBrillouin zone for J = 0 . t . t ′ = − . t and t ′′ = 0 . t have been assumed. and both of these interactions can be responsible for the destruction of the AF LRO.In order to determine the underlying mechanism we have carried out simulations takinginto account one of these interactions at a time. The results shown in the upper panelof Fig. 7 clearly demonstrate that the constraint plays the dominating role in thedestruction of the AF state. When both these interactions are taken into account theLRO is completely destroyed and only AF SRO can be observed for δ = 0 .
08. However,if we ignore the constraint LRO is restored. The results obtained without the intersitespin–hole interaction ( S zi M zj terms) are almost indistinguishable from the ones obtainedwith both the interactions. Note, that for λ = 0 the NDO constraint is completelyneglected. This illustrates the importance of the constraint for a realistic descriptionof the AF order at finite doping. Additionally, in order to determine the value of theLagrange multiplier for which the results converge and the NDO constraint is fulfilledwe have calculated g ( r ) for a wide range of λ . The results are presented in the lowerpanel of Fig. 7. One can see that the assumed value λ = 100 t is large enough to enforcethe constraint. After explaining the role of the farther neighbour hopping we restrictthe following analysis to the case of t ′ = − . t and t ′′ = 0 . t . sing t – J model close to half filling: A Monte Carlo study g ( r ) full Hamiltonianwithout S ⋅ M term λ = 0 r g ( r ) λ = 0 λ = 3 λ = 10 λ = 30 λ = 100 λ = 300 Figure 7. g ( r ) for J = 0 . t, kT = 0 . t, δ = 0 . t ′ = − . t and t ′′ = 0 . t . In theupper panel the topmost curve has been calculated with λ = 0. The other curves showresults obtained for λ = 100 t with and without the spin–hole exchange interaction. Thelower panel shows g ( r ) for different values of λ . Note, that the results for λ = 100 t and λ = 300 t are almost indistinguishable. Since our method works for systems with broken translational invariance, it is temptingto apply it as well to inhomogeneous systems. It has recently been shown with the help ofscanning tunneling spectroscopy that nanoscale electronic inhomogeneity is an inherentfeature of many groups of high–temperature superconductors. By a direct probing of thelocal density of states, these methods reveal strong spatial modulation of the energy gapin the superconducting Bi–based compounds [39]. Very recently, STM experiments haveshown a strong correlation between position of the dopant atoms and all manifestationsof the nanoscale electronic disorder [40, 41]. Thus, these experiments proved essentiallythat the impurities were the source of the inhomogeneity. On the other hand, theyrevealed a very important feature: there is a positive correlation between the magnitudeof a gap and the position of an out–of–plane oxygen atoms [40]. These are the atomswhich have been doped into the insulating parent compound in order to introduce holesto the CuO planes. The gap–impurity correlation has been explained as a result of theinhomogeneity–enhanced exchange interaction in the t – J model [42]. Assuming thatpurely electronic models contain the essential physics of cuprates, the same interactionis responsible for both superconductivity and the AF order. Therefore, inhomogeneitymay affect the AF state as well. Additionally, localization of holes by the electrostaticpotential of out–of–plane oxygen atoms may also affect the AF order, since the hopping sing t – J model close to half filling: A Monte Carlo study H λt − J | Z → H λt − J | Z + X iσ ε i d † iσ d iσ , (22)where ε i = V , if there is an out–of–plane oxygen atomabove site i ,0, otherwise.Since in HTSC each doped oxygen atom introduces one hole in the CuO plane, wehave carried out calculations for the number of impurities equal to the number of holes.Technically, for each Monte Carlo simulation we generate a random configuration of theout–of–plane oxygen atoms and keep it frozen during the whole run. In that way bothholes and localized spins feel a quenched disorder.In the upper panel of Fig. 8 we show a comparison of the correlation function g ( r )calculated in the presence of the diagonal disorder and without it. ¿From this figureone sees that the influence of the diagonal disorder is almost negligible, at least fora small–to–moderate values of the potential V . For larger values of V , the presenceof negatively charged out–of–plane oxygen atoms reduces the hole mobility resultingin a visible enhancement of the spin–spin correlation function. It is worthwhile toemphasize that the NDO constraint becomes very important in the presence of thediagonal disorder despite the low concentration of holes. In a homogeneous systemswith δ ≪ δ . Since the negatively charged out–of–plane oxygen atoms locallyenhance the hole concentration, neglecting of the NDO may significantly modify theresults [42].Now we turn to the influence of the inhomogeneity–induced enhancement of theexchange interaction. Following the results of Ref. [42] we assume J to a be site–dependent quantity: J ij = J (1 + η ij ) , (23)where η ij = η >
0, if there is an out–of–plane oxygen atomabove site i or j ,0, otherwise.In contrast with the superconducting gap [42], the AF order is hardly modified by thismechanism. This can be clearly inferred from the lower panel in Fig. 8. The regime fora magnetic ordering predicted by many calculations in the t – J model extends to muchlarger dopings than observed in cuprates and this discrepancy is sometimes attributedto the inhomogeneities, which are neglected in many theoretical approaches (see thediscussion in Ref. [14]). Although, inhomogeneities are expected to play an important sing t – J model close to half filling: A Monte Carlo study g ( r ) δ = 0.02, V = 0 δ = 0.02, V = - t δ = 0.02, V = -5 t δ = 0.06, V = 0 δ = 0.06, V = - t r g ( r ) δ = 0.02, η = 0 δ = 0.02, η = 0.1 t δ = 0.06, η = 0 δ = 0.06, η = 0.1 t Figure 8. g ( r ) for J = 0 . t in the presence of the out–of–plane oxygen atoms. Theupper panel demonstrates the influence of the diagonal disorder, whereas the lowerpanel shows effects coming from the site–dependent exchange interaction. The modelparameters ( δ, V, η ) are given in the legend. role in high–temperature superconductors [25], our results indicate that their influencein the AF ordering is rather limited. In particular, we expect that the inhomogeneitiesintroduced by the out–of–plane oxygen atoms cannot explain the above discrepancy. Since our analysis has been carried out on finite clusters, it is necessary to check towhat extent the results are affected by the finite size effects. One of the measures of thesignificance of the finite size effects is a sensitivity to the boundary conditions. Therefore,we have calculated the correlation functions and the spectral functions for systems withdifferent boundary conditions and compared them to those which have been obtainedwith periodic boundary conditions. Here, we have used a method known as averagingover boundary conditions (ABC) [28]. Each time a particular hole jumps out of thecluster, it is mapped back into the cluster with wave function with a different phase.Then the results are averaged over these phases; thereby the reciprocal space is probedin a much greater number of points than in the case of periodic boundary conditions. In sing t – J model close to half filling: A Monte Carlo study ×
40 cluster. Fig. 9 shows a comparison of spectral functionsobtained for 20 ×
20 and 40 ×
40 clusters with periodic boundary conditions as well asfor 20 ×
20 cluster with ABC. Since the false–colour plots of these spectral functions arevery similar to each other, we present their energy dependence for a few selected pointsof the Brillouin zone. One can see from this figure that the coherent part of the spectral -4 -3 -2 -1 0 1 20.00.10.2 A [ a . u . ] -3 -2 -1 0 1 2 3 ω /t A [ a . u . ] -4 -3 -2 -1 0 1 20.00.10.2 A [ a . u . ] × 20 × 40 ABC k =(0, 0) k =(π/2, π/2) k =(π, π) Figure 9.
Spectral functions calculated for selected points of the Brillouin zone. Theseresults have been obtained on 20 ×
20 (solid line) and 40 ×
40 (dashed line) clusters withperiodic boundary conditions and on 20 ×
20 cluster with ABC (dotted line). Since thepositions of the coherent peaks obtained with help of these three approaches are almostindistinguishable, we have cut the vertical axis in such a way that the incoherent partsare more pronounced. These results have been obtained for J = 0 . t , kT = 0 . t , and δ = 0 . functions is almost exactly the same in these three cases. The low–intensity parts also sing t – J model close to half filling: A Monte Carlo study g ( r ). Fig. 10 shows g ( r )determined for the same three systems, for which the spectral functions are presentedin Fig. 9. Despite minor quantitative differences the overall character of all three r g ( r )
20 × 2040 × 40
ABC
Figure 10.
Correlation function g ( r ) for the same systems as in Fig. 9. Also,the same parameters have been used. The insets show examples of snapshots of thespin configurations for 40 ×
40 and 20 ×
20 clusters [black (white) square corresponds to S zi = ( S zi = − )]. correlation functions is the same. Though we have not carried out a systematic finitesize scaling, the similarity of both the spectral functions and the correlation functionsconstitutes a significant indication that our results are valid also in the thermodynamiclimit.
4. Summary
We have developed a doped-carrier representation of the Ising t – J model. In thisformulation, the system is described in terms of fermions interacting with staticlocalized spins. Although it is a slave–particle approach, in contrast with many similarapproaches, the local NDO constraint is taken into account exactly. The proposedHamiltonian has the global Z symmetry at any values of the parameters, J and t . Thismodel is of an interest in itself since it represents a simple though nontrivial electronsystem which captures the physics of strong electron correlations. The issue of how thesecorrelations affect the magnetic ordering of the lattice spins is thouroughly investigatedin the present work. Besides, this model may provide at least for some values of themodel parameters a guess supported by unbiased numerical calculations regarding theactual low-energy behaviour of the quasiparticle excitations in a more realistic isotropic t − J model.In particular, we have calculated the one–hole energy and compared our results withthose obtained for t – J z and t – J models. We have found that the one–hole energy is the sing t – J model close to half filling: A Monte Carlo study t – J z model in two limiting cases, J → J → ∞ . For intermediate J the one–hole energy obtained for the t – J model is in between the results obtainedfor t – J z and t − J | Z approaches and the differences between t – J and t – J z models arecomparable to those between t – J and t − J | Z ones.The main advantage of the present approach consists in that it can be applied forvery large systems and the computational effort increases much slower with the size ofthe system than, e.g., for exact diagonalization. Moreover, it works for arbitrary valueof the coupling J and for arbitrary doping level. In particular, in the small– J regimethe exact diagonalization and Quantum Monte Carlo methods give rather poor results.This is because of the fact that in this regime the size of the defects generated by movingholes is comparable to or larger than the size of cluster the calculations can be carriedon. Since the clusters in our approach are much larger, this problem is less significant.Additionally, our method does not require translational invariance of the system. Thisfeature is especially important in the context of the recent experimental results, whichclearly indicate the presence of inhomogeneities in curates. It could also be applicableto optical lattices, where the translational symmetry is broken by a trap.Using the proposed approach we have found that the AF SRO persists fortemperatures and dopings which are much beyond the boundaries of the AF LRO phase,what is in agreement with recent experiments on the high–temperature superconductors.We have also demonstrated that the AF LRO depends on the exchange interaction J .It concerns the transition temperature as well as the maximal doping at which the AFLRO vanishes. We explicitly demonstrate that the local no double occupancy constraintplays the dominant role in destroying the magnetic order at finite doping.Finally, we have shown that the inhomogeneities induced by the out–of–planeoxygen atoms have a rather limited influence on the spin–spin correlations functions,at least in the underdoped regime. Although, localization of holes by their electrostaticpotential stabilizes the AF LRO, this mechanism becomes important only for a relativelystrong diagonal disorder. Such a limited influence of inhomogeneities on the AF orderis closely related with the NDO constraint. Note that exactly at half filling the diagonaldisorder does not influence the system, provided the NDO is properly taken into account. Acknowledgments
This work has been supported by Bogolyubov–Infeld program and by the Polish Ministryof Education and Science under Grant No. 1 P03B 071 30.
References [1] Schmitt-Rink S, Varma and Ruckenstein A E 1988
Phys. Rev. Lett. , 2793[2] Shraiman and Siggia E D 1988 Phys. Rev. Lett. , 740[3] Trugman S A 1988 Phys. Rev. B , 1597[4] Wr´obel P and Eder R 1998 Phys. Rev. B , 15160[5] Dagotto E, Joynt R, Moreo A, Bacci S and Gagliano E 1990 Phys. Rev. B , 9049 sing t – J model close to half filling: A Monte Carlo study [6] Prelovsˇek P, Sega I and Bonˇca J 1990 Phys. Rev. B , 10706[7] Martinez G and Horsch P 1991 Phys. Rev. B , 317[8] Poilblanc D, Ziman T, Schulz H J and Dagotto E 1993 Phys. Rev. B , 14267[9] Leung P W and Gooding R J 1995 Phys. Rev. B , R15711[10] Brunner M, Assaad F F and Muramatsu A 2000 Phys. Rev. B , 15480[11] Bonˇca J, Maekawa S and Tohyama T 2007 Phys. Rev. B , 035121[12] Lugas M, Spanu L, Becca F and Sorella S 2006 Phys. Rev. B , 165122[13] Morinari T 2004 JMMM , 188[14] Spanu L, Lugas M, Becca F and Sorella S 2008 Phys. Rev. B , 024510[15] Dagotto E 1994 Rev. Mod. Phys. , 763[16] Chernyshev A L, White S R and Castro Neto A H 2002 Phys. Rev. B , 214527[17] Chernyshev A L, Castro Neto A H, and White S R 2005 Phys. Rev. Lett. , 036407[18] Riera J A 2001 Phys. Rev. B , 104520[19] Martins G B, Gazza C and Dagotto E 2000 Phys. Rev. B , 13926[20] Chernyshev A L and Leung P W 1999 Phys. Rev. B , 1592[21] Elrick B M and Jacobs A E 1995 Phys. Rev. B , 10369[22] Ribeiro T C and Wen X–G 2005 Phys. Rev. Lett. , 057001; Ribeiro T C and Wen X–G 2006 Phys. Rev. B , 155113[23] Ferraz A, Kochetov E and Mierzejewski M 2006 Phys. Rev. B , 064516[24] Pepino R T, Ferraz A and Kochetov E 2008 Phys. Rev. B , 035130[25] Dagotto E 2005 Science , 257[26] Chao K A, Spa lek J and Ole´s A M 1978
Phys. Rev. B , 3453[27] Ma´ska M M and Czajka K 2006 Phys. Rev. B , 035109[28] Poilblanc D 1991 Phys. Rev. B , 9562; Gros C 1992 Z. Phys. B , 359; 1996 Phys. Rev. B ,6865; Czajka K and Ma´ska M 2004 Acta Phys. Pol. A , 703[29] Riera J and Dagotto E 1993 Phys. Rev B. Phys. Rev. B , 10977[31] Ba la J, Ole´s A M and Zaanen J 1995 Phys. Rev. B , 4597; Xiang T and Wheatley J M 1996 Phys. Rev. B , R12 653; Kyung B and Ferrell R A 1996 Phys. Rev. B , 10125; Lee T K andShih C T 1997 Phys. Rev. B , 5983; Lee T K, Ho C–M and Nagaosa N 2003 Phys. Rev. Lett. , 067001; Danmascelli A, Shen Z–X and Hussain Z 2003 Rev. Mod. Phys. , 473[32] Regnault L P, Bourges P, Burlet P 1998 Neutron Scattering in Layered Copper-OxideSuperconductors , edited by A. Furrer (Kluwer, Amsterdam); Kastner M A, Birgeneau R J,Shiran G and Endoh Y 1998
Rev. Mod. Phys. , 897; Aeppli G, Mason T E, Hayden S M,Mook H A and Kulda J 1997 Science , 1432; Curro N J, Imai T, Slichter C P and DabrowskiB 1997
Phys. Rev. B , 877[33] van den Brink J and Sushkov O P 1998 Phys. Rev. B , 3518[34] Wells B O, Shen Z–X, Matsuura A, King D M, Kastner M A, Greven M and Birgeneau R J 1995 Phys. Rev. Lett. , 964; Pothuizen J J M, Eder R, Hien N T, Matoba M, Menovsky A A andSawatzky G A 1997 Phys. Rev. Lett. , 717[35] Cataudella V, De Filippis G, Mishchenko A S and Nagaosa N 2007 Phys. Rev. Lett. , 226402[36] Bonˇca J, Maekawa S, Tohyama T and Prelovˇsek P 2008 Phys. Rev. B , 054519[37] Shen K M, Ronning F, Lu D H, Lee W S, Ingle N J C, Meevasana W, Baumberger F, DamascelliA, Armitage N P, Miller L L, Kohsaka Y, Azuma M, Takano M, Takagi H and Shen Z–X 2004 Phys. Rev. Lett. , 267002[38] Shen K M, Ronning F, Meevasana W, Lu D H, Ingle N J C, Baumberger F, Lee W S, Miller L L,Kohsaka Y, Azuma M, Takano M, Takagi H and Shen Z–X 2007 Phys. Rev. B , 075115[39] Pan S H, O’Neal J P, Badzey R L, Chamon C, Ding H, Engelbrecht J R, Wang Z, Eisaki H, UchidaS, Guptak A K, Ngk K–W, Hudson E W, Lang K M, and Davis J C 2001 Nature , 282[40] McElroy K, Lee J, Slezak J A, Lee D–H, Eisaki H, Uchida S and Davis J C 2005
Science ,1048 sing t – J model close to half filling: A Monte Carlo study [41] Mashima H, Fukuo N, Matsumoto Y, Kinoda G, Kondo T, Ikuta H, Hitosugi T and Hasegawa T2006 Phys. Rev. B , 060502(R)[42] Ma´ska M M, ´Sled´z ˙Z, Czajka K and Mierzejewski M 2007 Phys. Rev. Lett. , 147006[43] Czajka K and Ma´ska M M 2007 phys. stat. sol. (b)244