Monte Carlo simulation of a model cuprate

We develop a classical Monte Carlo algorithm based on a quasi-classical approximation for a pseudospin S=1 Hamiltonian in real space to construct a phase diagram of a model cuprate with a high Tc. A model description takes into account both local and nonlocal correlations, Heisenberg spin-exchange interaction, correlated single-particle, and two-particle transport. We formulate a state selection algorithm for a given parameterization of the wave function in order to ensure a uniform distribution of states in the phase space. The simulation results show a qualitative agreement with the experimental phase diagrams.


I. INTRODUCTION
One of the topical problems of the high-T c cuprate physics is the coexistence and competition of antiferromagnetic, superconducting, and charge orderings [1]. Recent accurate measurements of various physical characteristics on thousands of cuprate samples [2] indicate fundamental discrepancies with the ideas based on the canonical Bardeen-Cooper-Schrieffer approach, and rather support the bosonic mechanism of the high-T c cuprates. The study is complicated by the presence of heterogeneity due to dopants or non-isovalent substitution, as well as to the internal electronic tendency to heterogeneity [3]. A large number of theoretical models have been developed to account for the exotic electronic properties of cuprates in the normal state and to reveal the nature of unconventional superconductivity. However, to date, the problem of a consistent theoretical description of the cuprate phase diagram is still far from being solved.
Earlier, we developed a minimal model of cuprates where the CuO 2 planes are considered as lattices of CuO 4 clusters, which are the main element of the crystal and electronic structure of cuprates. The on-site Hilbert space [4,5] is formed by three effective valence states of the cluster: [CuO 4 ] 7− , [CuO 4 ] 6− , and [CuO 4 ] 5− . The very possibility of considering these centers on an equal basis is due to the strong electron-lattice relaxation effects in cuprates [6,7]. The centers differ in the usual spin state: s = 1/2 for the [CuO 4 ] 6− center and s = 0 for the [CuO 4 ] 7− and [CuO 4 ] 5− , respectively, and in the orbital symmetry: B 1g for the ground states of the [CuO 4 ] 6− center, A 1g for the [CuO 4 ] 7− center, and the Zhang-Rice A 1g or more complicated low-lying non-Zhang-Rice states for the [CuO 4 ] 5− center. In these many-electron atomic complexes with strong p−d covalence and strong intra-center correlations, electrons cannot be described within any conventional (quasi)particle approach that addresses the [CuO 4 ] 7−,6−,5− centers within the on-site hole representation |n , n = 0, 1, 2, respectively. Instead of conventional quasiparticle k-momentum description, we make use of a real space on-site S = 1 pseudospin formalism to describe the charge triplets and an effective spinpseudospin Hamiltonian which takes into account both local and nonlocal correlations, single and two-particle transport, as well as Heisenberg spin-exchange interaction.
The pseudospin approach is used for the strongly correlated electron systems [8,9] and for the superconductivity [10] of cuprates for a long time. The pseudospin formalism leads to the possibility of simulation using the well-developed classical Monte Carlo (MC) method for constructing phase diagrams and studying the features of the thermodynamic properties of the system. A similar effective S = 1 spin-charge model for cuprates and its MC implementation were considered in recent papers [11,12].
We organize the article as follows. In Section 2, we present the S = 1 pseudospin formalism and the effective spin-pseudospin Hamiltonian of the model. In Section 3, we introduce quasi-classical approximation for the wave functions and formulate the state selection algorithm. The results of classical MC simulations of our model and their discussion are presented in Section 4.

II. MODEL
A minimal model to describe the charge degree of freedom in cuprates [4,5] implies that for the CuO 4 centers in CuO 2 plane the on-site Hilbert space reduced to a charge triplet formed by the three many-electron valence states [CuO] 7−,6−,5− 4 (nominally Cu 1+,2+,3+ ). These states can be considered to be the components of the S = 1 pseudospin triplet with projections M S = −1, 0, +1. Effective pseudospin Hamiltonian of the model cuprate with the addition of the Heisenberg spin-spin exchange coupling of the s = 1/2 [CuO] 6− 4 (Cu 2+ ) centers can be written as follows: (1) Here, the first term describes the on-site and inter-site nearest-neighbour density-density correlations, respectively, so that ∆ = U/2, U being the correlation parameter, and V > 0. The sums run over the sites of a 2D square lattice, ij means the nearest neighbors. The second term is the antiferromagnetic (J > 0) Heisenberg exchange coupling for the CuO 6− 4 centers, where σ = P 0 s/s operators take into account the on-site spin density P 0 = 1 − S 2 z , and s is the spin s = 1/2 operator. The third term has the following form: where the transfer integrals t p , t n , t pn describe three types of the correlated "one-particle" transport. P and N operators are the combinations of the pseudospin S = 1 operators [4]: where the transfer integral t b describes the two-particle ("composite boson") transport [4]. The last term with the chemical potential µ allows to account the constraint for the total charge density n: The mean-field approximation (MFA) for a model with the Hamiltonian (1) allowed us to find [13] the critical temperatures equations for the antiferromagnetic (AFM) ordering, charge ordering (CO), superconducting ordering (SC), and for the "metal" phase (M). The MFA phase diagram [14,15] of the model (1) demonstrates the possibility of correctly describing the features of phase diagrams typical of cuprates.
In the quasi-classical approximation, we write the on-site wave function of the charge triplet as follows where |+1 and |−1 are the orbital wave functions of [ and the coefficients can be written in the following form: Here The matrices of pseudospin operators S z and S ± in a basis {|+1 , |0 , |−1 } are written as Using equations (7-9) we can write average values for all operators in the Hamiltonian (1) on the ith site: where the x-, y-, and z-components of vector are listed in curly brackets. This allows us to obtain the energy for a model (1) in the quasi-classical approximation in the following form: The constraint (6) can be written as:

III. STATE SELECTION ALGORITHM
An average value of the spin s = 1/2 operator σ|s|σ = 1 2 {sin η cos χ, sin η sin χ, cos η} (20) maps the spin states to the unit sphere. It is well known, that the uniform distribution of randomly generated points over the unit sphere is given by the following state selection algorithm: where γ 1,2 are random numbers in the [0, 1] range. The (χ, η)-histogram is shown in Fig.1, left panel, and this produces the flat (χ, z)-histogram shown in Fig.1, right panel. This state selection algorithm is based on the well-known lemmas [16] of probability theory: is the marginal density, and f 2 (x 2 |x 1 ) is the conditional density. Let ξ 1 and ξ 2 be continuous random variables that satisfy the system of equations Then ξ 1 and ξ 2 have the joint density function f (x 1 , x 2 ).
The states with coefficients (9) corresponds to a point in the octant of the unit sphere. We use the Metropolis algorithm for a system with conservation of the total charge. The charge at the site, n i , is related to the parameters of the wave function by the expression We require when the states of sites 1 and 2 change simultaneously, the total charge of the pair is preserved, n 1 + n 2 = n 1 + n 2 = 2n, and the points representing states uniformly fill the allowed area in the octant. The state selection algorithm based on the Lemmas 1 and 2 consists of the following steps: 1. caclulation of n 1 , where −1 + n + |n| ≤ n 1 ≤ 1 + n − |n|, from equation where γ ∈ [0, 1] is the uniformly distributed random value, is the complete elliptic integral of the third kind, K(m) is the complete elliptic integral of the first kind; 2. calculation of the value n 2 = 2n − n 1 ; 3. calculation of cos θ i 2 from equation where γ i ∈ [0, 1], i = 1, 2, are the uniformly distributed random values, sn (x, m) is the Jacobi function. If n i = 0, we take cos θ i 2 = γ i . 4. calculation of cos φ i from equation If n i = 0 and cos θ i 2 = 1, φ i is a random uniformly distributed quantity, 0 ≤ φ i ≤ π. The distributions in the case n = 0 for angles (φ, θ) and for states on the octant of unit sphere are shown in Fig.2, left and right panels, respectively.

IV. RESULTS
In MC simulation, we calculated the structure factors where A l and B m are the on-site operators and the summation is performed over all sites of the square lattice. To determine the type of ordering, we monitored the following structure factors: • F (π,π) (σ, σ) for antiferromagnetic (AFM) order, • F (π,π) (S z , S z ) for the charge order (CO), • F (0,0) (S 2 + , S 2 − ) for the superconducting order (SC), • F (0,0) (P + , P ) for the "metal" phase (M ).
The results of the MC simulation for the doping dependencies of the main structure factors near the ground state, T /t b = 0.05, of the model cuprate are presented in Fig. 3, left panel. The critical temperatures for the AFM, CO, and SC phases were determined from the jump in the structure factor from zero to a certain finite value. Fig. 3, right panel, shows the reconstruction from the MC simulations of the T -x phase diagram. The insert shows the typical for the hole doped cuprates [17].
The obtained phase diagram for model cuprate with the Hamiltonian (1) reproduces some most important features of real phase diagrams: with given parameters, ∆ = 0.8, V = 0.625, J = 1, t p = 0.35, t n = 0, t pn = −0.24, (all in units of the t b ), near the "parent" composition x = 0, we get AFM ordering, which is replaced with increasing x by the SC ordering, that coexists with the CO ordering in the form of phase separation. The found SC ordering exits in the finite region of doping, and it is replaced by the "metal" phase at x ≥ 0.3. The main problems in the implementation of our modeling, such as inhomogeneous phase states and the associated difficulties in identifying them, are predetermined by the enhanced role of fluctuations in low-dimensional systems. But at the same time, the obtained phase diagrams show promising possibilities to describe the coexistence and competition of various phase orders in cuprates.