A three-state model with competing antiferromagnetic and pairing interactions

Motivated by the rich phase diagram of the high temperature superconductors, we introduce a toy model with three state variables which can be interpreted as two state particles and holes. The Hamiltonian has a term which favors antiferromagnetism and an additional competing interaction which favors bonding between pairs of antiparallel spins mediated by holes. For low concentration of holes the dominant interaction between particles has antiferromagnetic character, leading to an antiferromagnetic phase in the temperature-hole concentration phase diagram, qualitatively similar to the antiferromagnetic phase of doped Mott insulators. For growing concentration of holes antiferromagnetic order is weakend and a phase with a different kind of order mediated by holes appears. This last phase has the form of a dome in the T-hole concentration plane. The whole phase diagram resembles those of some families of high $T_c$ superconductors. We compute the phase diagram in the mean field approximation and characterize the different phase transitions through Monte Carlo simulations.


I. INTRODUCTION
One of the frontiers of modern condensed matter physics is the description of phase transitions in complex many body systems that challenge the well known theories of classical fluids or the well established quantum theories of Fermi or Bose liquids. In particular, when two or more competing interactions are present, usually strong correlations play a significant role in the physics at low temperatures and the phase diagram can show a rich variety of different phases, with different symmetries and thermodynamic properties. Competing interactions are common in nature, an important example are magnetic systems in which different types of magnetic interactions usually coexist like exchange, dipolar, Dzyaloshinsky-Moriya and Zeeman [1][2][3][4][5]. Each of these terms usually lead to a particular kind of order at low temperatures, but when at least two of them are simultaneously present a much richer behaviour appears. Competition between interactions can lead to pattern formation and phases with complex symmetries [6][7][8][9]. High temperature superconductors appear as another example of extrememly complex and fascinating systems. Since their discovery in the 80's thousands of papers have been published on the subject, but after 30 years of very intense research both the experimental and theoretical situations are still not satisfactory [10][11][12][13][14]. From the theoretical perspective, perhaps the biggest open question is the nature of the microscopic pairing mechanism which leads electrons to form Cooper pairs in the high T c superconductor families of compounds [10,15]. It is accepted that the BCS theory, so successful at explaining superconductivity in the "usual" superconductors [16,17], does not explain the superconductiong behavior of the cuprates or the iron based families. In the cuprates or iron based supercon-ductors strong electronic correlations seem to be at work forcing the theoretical description of the relevant physics to go beyond the classic BCS theory. One of the striking effects of the strong electronic correlations is the simultaneous relevance and consequent competition between different degrees of freedom, leading to possibly different kind of orders associated to electron density , spin and orbital degrees of freedom. It is not even completely clear when these different orders compete with each other or instead in some sense "catalyse" the appearence of superconductivity. Sometimes this is refereed to as "intertwinned orders" [12,18]. From this simple perspective it is immediately clear that the phase diagrams of these systems are extremely rich and complex. In fact, both experimental and theoretical work usually focus on one or few relevant observables in limited regions of parameter space. An overall understanding of the phase diagram, collecting all or most of the pieces at work is still not available.
The more well known families of high T c superconductors are the "cuprates" [12,13] and the iron based ones or "pnictides" [11,19]. Both families show similar thermal properties, although not identical. In both cases superconductivity arises by doping with electrons or holes the so called parent compounds. The undoped parent compounds are insulators at low temperatures, showing antiferromagnetic order at half filling. As long as the doping is introduced in the system antiferromagnetism gradually gets weaker and eventually dissapears. Beyond some degree of doping, which depends on the compound, a superconducting phase appears at very low temperatures. By increasing the doping the superconducting phase extends to higher temperatures until a maximum T c is obtained at an "optimal" doping level. Beyond optimal doping the critial superconducting temperature goes down and arXiv:1810.11373v1 [cond-mat.supr-con] 26 Oct 2018 eventually goes to zero at another critical doping level. This is the well known "superconducting dome". In general, in the cuprate superconductors antiferromagnetism dissapears for a doping level less than that at which superconductivity appears, the antiferromagnetic phase and the superconducting dome do not intersect in the temperature-doping plane. On the other side, in the iron based superconductors the antiferromagnetic and superconducting phases usually appear superposed, leading to a region of coexistence between antiferromagnetism and superconductivity. In the same temperature-doping plane of the phase diagram other phases, associated with other degrees of freedom and different symmetries usually appear, like modulated electronic and spin orders, called stripes and nematic phases [14,20,21]. Structural transitions in the crystal symmetry of the compounds are also observed probably being relevant for the emergence of superconductivity or other types of order seen in the phase diagram. The common presence of quenched disorder in the samples also leads to freezing of degrees of freedom and spin glass like behavior in some cases. From this crude exposition of the thermal phenomenology of high T c superconductors one can immediately conclude that obtaining a complete phase diagram is a formidable task.
From a statistical mechanics point of view it is usually very helpful to develop a simplified model which captures some of the relevant properties of a real system. This approach has proved to be extremely successful in the context of critical phenomena after the appearance of the concept of universality. In a nutshell, universality means that the behavior of the order parameters and the associated responses to the conjugate fields of any system near a continuous phase transition depend only on a few relevant charateristics of the system, basically the symmetries of the order parameter and the dimensionality of space. In particular, microscopic details of the interactions are not relevant regarding the critical behaviour. A well known example is the Blume-Emery-Griffiths (BEG) model [22] which, despite its extreme simplicity, com-pletely describes the phase diagram global topology of 3 He and 4 He mixtures, including the right order of the phase transitions. Extensions of the BEG model has also been proposed to describe some aspects of superconductivity [23]. With this in mind and motivated by the rich complexity of the phase diagram of high T c superconductors, in this paper we introduce a very simplified model which captures the topology of the antiferromagnetic and superconductor phases of those systems. To our knowledge, these two phases have only been obtained from a single model Hamiltonian in very few cases e.g. [24,25] and references therein (see also [23]), and because of the complexity of the models involved and the different nature of the order parameters a detailed analysis of the phase transitions have not been done yet. Then, although the present model is clearly too simplified to describe the complex physics of the cuprates or iron based superconductors, we nevertheless think it can be useful to think on the universal or robust thermal properties which can lead to the particular phase diagrams of such systems.

II. THE MODEL
The model is defined by a set of three-state classical variables S i,j = 0, ±1, i, j = 1 . . . L, interacting in a square lattice of N = L × L sites, where S i,j = 0 represents a hole and the two other states S i,j = ±1 may be thought of as representing the spin states of an electron. The Hamiltonian has two terms which compete with each other: a nearest-neighbour interaction of strength J A which favors antiferromagnetic order, and a three-site interaction mediated by a hole, of strength J B . If the concentration of holes is zero this interaction is not present and the system shows only an antiferromagnetic phase of Néel type at low temperatures. The presence of a hole favors another kind of antiferromagnetic order, this time between third neighbours along the two principal directions of the square lattice with a hole in between. The Hamiltonian is given by This order parameter is able to capture emergent order induced by the second interaction term in (1).

III. MEAN FIELD PHASE DIAGRAM
The grand partition free energy density is given by: where H = H − µM and is the total number of particles and µ the corresponding chemical potential. A mean field free energy Φ can be obtained using the Bogoliubov inequality: where H 0 is a reference non-interacting variational Hamiltonian with and its partition function. H 0 assumes the general form: where µ i,j and η i,j are local variational parameters. and A. Antiferromagnetic solution We consider first the antiferromagnetic solution η i,j = ±η, according to which sublattice the site of coordinates i, j belongs, and µ i,j = µ. This implies a constant density: and a staggered manetization S i,j 0 = ±m with m = 2 e βµ sinh (βη) 1 + 2 e βµ cosh (βη) = 2 sinh (βη) e −βµ + 2 cosh (βη) . (13) It is easy to obtain Z i,j 0 = 1 1−ρ and therefore f 0 = 1 β ln(1− ρ).
Then, the variational free energy density of the antiferromagnetic phase reads: Solving for the stationary condition on the variational parameter η, ∂φ AF ∂η = 0 gives: Inserting this result in (12) and (13) we obtain two coupled equations for the order parameters ρ and m which must be solved self-consistently. From (15) m = 0 is always a solution. In the limit µ → ∞ we have ρ = 1 and m = tanh(4J A βm), as expected. The numerical solution of the self consistent equations for m and ρ suggests that for large positive values of the chemical potential µ there is a line of continuous transitions between the antiferromagnetic and paramagnetic phases, which ends at a tricritical point. For locating the tricritical point we did a Landau expansion of the free energy density: The second order transition line is defined when A(µ, T ) = 0. This gives: Changing variables to T = 1/β (in units of the Boltzmann constant k B ) and defining the density of holes x = 1 − ρ, we obtain the critical line T (x): This should be the phase diagram of the antiferromagnetic phase assuming a continuous transition for any x. Nevertheless, it is possible to show that the continous transition line ends at a tricritical point (µ t , T t ). This is identified imposing the simultaneous vanishing of the prefactors A(µ, T ) = 0 and B(µ, T ) = 0. For the case J B = 2J A the tricritical point can be found numerically to be at (µ t , T t ) = (1.77, 2.55). For µ < µ t the line of transitions is of first order. This is illustrated in the T −µ phase diagram in Figure 2.

B. Super-Antiferromagnetic solution
The SAF ground state has the symmetry shown in Figure 1. The system is divided in two interpenetrated sublattices A and B (black and red-green sites in the figure). Sublattice A has a uniform low-density, low-magnetized state. Sublattice B has a striped state (red and green sites in Figure 1). Then, for the SAF phase we choose the variational parameters to be µ i,j = µ + δ and η i,j = 0 for all sites (i, j) in the sublattice A. For all sites (i, j) in the sublattice B we set µ i,j = µ and η i,j = ±η for sites belonging to alternated columns. Then, for sites belonging to the sublattice A S i,j A 0 = m A = 0 and: For sites belonging to the sublattice B: and 2e βµ cosh (βη) 1 + 2e βµ cosh (βη) = 2 cosh (βη) e −βµ + 2 cosh (βη) .
In this case f 0 = − 1 2β ln 1 + 2e β(µ+δ) + ln 1 + 2e βµ cosh (βη) , and the variational free energy for the SAF phase results: The stationarity conditions ∂φ SAF ∂η = 0 and ∂φ SAF ∂δ = 0 give: and Inserting (24) and (25) in (19), (20) and (21) we end with three coupled equations for ρ A , ρ B and m B to be solved self consistently. A complete solution can only be possible by solving numerically the set of equations. Nevertheless, a simple analysis of some limiting cases is useful to check for consistency of the equations: • The global density in the SAF state is defined to be ρ = (ρ A + ρ B )/2.
• There exists a minimum value of the chemical potential µ m such that for µ < µ m the only stable phase is the disordered one characterized by The mean field SAF phase is shown in Figure 2. From the numerical solution of the saddle point equations it turns out that the SAF order parameter m B always changes discontinuously at the transition line. In the mean field approximation the SAF-disordered transition is discontinous for any value of the chemical potential µ. This situation changes in the results from Monte Carlo simulations, as will be shown in the next section. The SAF phase has a dome like shape in the T − µ or T − x plane. In Figure 2 it can be seen that the AF and SAF phases have a large coexistence region. Both first order transition lines meet approximately at (T, µ) ≈ (2.0, 0.6). For T < 2.0 the AF and SAF phases are separated by a first order line, which is approximately a straight line at µ = 0.6 in the mean field approximation. Overall, the phase diagram has a similar shape with the ones of the iron based high temperature superconductors [24,25].

IV. MONTE CARLO SIMULATIONS
We performed Grand Canonical Monte Carlo simulations using Metropolis algorithm with Hamiltonian (1) on a square lattice with N = L × L sites and periodic boundary conditions. For each set of parameters values we let first the system to equilibrate over 2 × 10 4 Monte Carlo Steps (MCS) and then we average over M s sample points taken every 100 MCS over a single MC run. M s run between 1000 and 20000. All the calculation, except those related to hysteresis cycles, were obtained by cooling down at constant steps from high temperature keeping the chemical potential constant and taking the initial configuration at every temperature as the last one of the previous temperature value. We calculated the average antiferromagnetic staggered magnetization m s , the orientational order parameter O , with O given by Eq. (2), the associated susceptibilities (k B = 1) and the density In all cases we used J B /J A = 2. We obtained the phase diagram in the (µ, T ) plane. We summarize first the main qualitative features of the different regimes observed before presenting a more detailed analysis.
• There exists a minimum value of the chemical potential µ m such that for µ < µ m the only stable phase is the liquid (disordered) characterized by O = m s = 0. For J B /J A = 2 and L = 64, µ m ≈ −3.86025, consistent with the mean field result.
• In the interval between µ m and µ ≈ −1 a low temperature ordered SAF phase was observed, characterized by O = 0 and m s ≈ 0.
• According to the value of µ the phase transition from the SAF to the liquid phase can be of first or second order with the presence of a tricritical point. Namely, for µ m < µ < µ t the transition is first order and there is phase coexistence; for µ > µ t the transition is continuous, up to a region around µ ≈ −1. The estimated tricritical values for L = 64 are µ t ≈ −2.5, x t ≈ 0.5 and T t ≈ 1.46.
• For µ > 0 there is a stable low temperature AF phase, characterized by O ≈ 0 and m s = 0. The liquid-AF transition is continuous.
• Around µ = 0 there is a first order phase transition between the ordered phases, with strong hysteresis effects.
These results are summarized in the phase diagram shown in Fig.3. We next describe the different properties analyzed to obtain that phase diagram. In Fig.4 we illustrate the typical behavior for µ 0. In this region the staggered magnetization becomes different from zero at low enough temperatures, while the orientational order parameter remains almost zero for all temperatures. The finite size scaling behavior of the staggered susceptibility max χ s ∼ L γ/ν allows us to estimate the critical exponent γ/ν = 1.72±0.05, in agreement with the exact value for the 2D Ising model γ/ν = 7/4 = 1.75, as expected. At variance with the mean field results, we did not observe evidences of a first order AF-liquid phase transition, nor of a tricritical point. However, large fluctuations close to the region where the transition line joins the SAF-liquid one makes difficult to exclude the existence of a tricritical point in that region.
Next we analyzed the order parameters behavior for µ m < µ < µ t . In Figs.5 and 6 we show how both the density and the orientational order parameter display a discontinuity that disappears when µ → µ t . For µ > µ t the density is always continuous. The AF order parameter remains almost zero in this region. In Fig.7 we show the finite size scaling of the orientational order parameter and the associated susceptibility, whose maximum max χ O ∼ L 2.03 agrees with the expected behavior in a first order transition [26] In Fig.8 we show an example of the finite size behavior of the orientational order parameter in when µ > µ t , but not very close to µ = 0, where the transition between both ordered phases occurs. The AF order parameter remains almost zero in this region. We see that the maximum of the orientational susceptibility scales with an exponent γ/ν ≈ 1.4, well below both from d and from 7/4, consistently with a continuous phase transition in a universality class different from that of the 2D Ising model. However, this last conclusion should be checked for other values of µ in the same region, since this value could be influenced by the proximity of the tricritical point. But this is a complicated and tricky analysis, since we have to discriminate both the finite size and the crossover effects related to the tricritical point. For instance, in Fig.9 we show the finite size scaling behavior of the maximum of χ O for µ = −1 and sizes ranging from L = 16 and L = 256. If we consider all the points, the power law fitting gives and exponent γ/ν = 1.5 ± 0.1. So, it seems that when we depart from the tricritical point the exponent increases, suggesting an approach to the 2D Ising universality class. However, if we discard the smallest size point we get a better fitting (larger r 2 ) with an exponent γ/ν = 1.35 ± 0.07, which appears to be consistent with the previous value found for µ = −2.
Finally, we analyzed the transition from the AF to the SAF phase at low temperature. In general we observed that the behavior in the transition region is very noisy and relaxation times are rather large. Hence, it becomes extremely difficult to determine the nature of the transition, as well as to get a reliable estimation of the transition temperature at fixed chemical potential by means of the standard methods, such as energy and/or order parameters histograms, response functions finite size scaling, etc. Since according to the MF results the transition is a first order one, we can at least check that by looking at hysteresis effects. We then performed chemical potential cycles at constant temperatures in the corresponding region of the phase diagram according to the following protocol. We first let the system to equilibrate M e MCS at high temperatures in the liquid phase (T=1.5). Then we cool the system down to a final temperature T f at constant steps ∆T = 0.05 for fixed value of the chemical potential µ = µ min = −2, letting the system to evolve M e MCS at every temperature value. The value of T f was chosen in such a way that the system starts the cycle well inside the SAF phase. Then we increase the chemical potential at constant ∆µ = 0.05 steps keeping the temperature fixed up to a maximum value µ = µ max = 2 (well inside the AF phase). For every value of µ we first let the system to evolve M e MCS and then we calculate the average of both order parameters, taking 10 4 sampling points every 100 MCS along the same MC run. Once we reach µ max we repeat all the protocol decreasing µ down to µ min . Along the whole process we used M e = 4 × 10 5 . The results are illustrated in Fig.10 for T f = 0.4. We observe a strong hysteresis effect, consistently with a first order transition. From these data we got a rough estimation of the transition temperature as an average of the µ coordinate of the center of mass of the cycles for every order parameter. The transition points shown in Fig.3 were obtained in that way. The associated error bars were estimated as the average half-width of the cycles. Actually, that procedure overestimates the error and those bars rather give an estimation of the location of the spinodal lines.

V. DISCUSSION
We introduced an effective lattice-gas model of spin 1/2 particles doped with holes, to analyze the effect of the competition between an exchange antiferromagnetic interaction and a pairing interaction mediated by holes. Through this toy model we showed that those two ingredients are enough to reproduce the main features of the global topology of the phase diagram of many high T c superconducting compounds. Besides some differences in the order of the involved transitions, both mean field and Monte Carlo analysis provide a consistent phase diagram, with a dome in the transition line from the paired induced phase (SAF) and a decreasing AF-liquid transition line as the doping increases (decreasing chemical potential). In order to turn the model more closely related with the universal characteristics of high T c superconductors, different modifications and/or generalizations can be envisaged. One possibility is to include another antiferromagnetic interaction between second neighbours in the square lattice, like in the J 1 − J 2 model [27,28]. The J 1 − J 2 model has a SAF phase, which is interpreted as spin or electronic stripes, similar to those observed in the cuprates and pnictides. Furthermore, in presence of an external field, the model shows a kind of "electron nematic" phase, also usually present in the phase diagram of high T c compounds [29]. The introduction of holes or doping in the cited models can be an interesting route to explore.