Spin-glass model with essential short-range competing interactions

We consider a proton-glass model with an arbitrary range of competing interactions within two-particle cluster approximation for the free energy and within symmetric replica approach. We show that for the thermodynamic characteristics, the Gaussian approximation for distribution functions of cluster fields can be used. For a simpler proton-glass model with weak long-range interactions (linear approximation) we study the effect of the long-range interactions and random internal fields on the phase diagram and thermodynamics characteristics. It is shown that small fluctuations of internal fields can lead to essential smearing of proton-glass transition region. The adaptability of the theory for a description of proton-glasses of Rb1−x(NH4)xH2PO4-type is discussed.


Introduction
The systems with spin glass phase transitions have been studied for about 30 years.The hydrogen-bonded crystals of the Rb 1−x (NH 4 ) x H 2 PO 4 type, in which a spin glass state is realized at low temperatures, have been intensively studied both experimentally [1][2][3] and theoretically [2,[4][5][6][7] for more than a decade.
In early works [4,5] the energy structure of a cluster, consisting of the PO 4 tetrahedron and protons on hydrogen bonds was taken into consideration.In subsequent papers simple spin models with an infinite range of the mean value and dispersion of random interactions were used [2,6,7].An adequate consideration of all types of interactions (both short-range and long-range ones) in such systems remains an important theoretical problem.
The major features of spin-glass systems can be described based on an Ising-like Hamiltonian H with random parameters h i and k ij In the case of Bravais lattice the indices i and j denote the sites of the lattice.
In the Sherrington-Kirkpatrick model (SK) [8,9] a Gaussian distribution for k ij is used, whereas their mean value k ij and dispersion (pair cumulant) do not depend on i, j.Thus this model corresponds to a system with long-range interactions.The phase diagram, magnetization, parameter of spin glass ordering q, entropy, and specific heat have been obtained for this model within the framework of symmetric replica approach.However, the entropy and specific heat calculated within this approach become negative at low temperatures.In [10] the replica symmetry breaking solution for the SK model was suggested, which is stable at low temperatures.Later, the correctness of the replica symmetry breaking solutions was confirmed by numerical simulations.
To describe the spin glass transition in hydrogen bonded systems of the Rb 1−x (NH 4 ) x H 2 PO 4 type (proton glass) we should consider the internal random fields with a Gaussian dispersion ∼ x(1 − x).These fields are induced by a substitutional disorder due to the difference in the ionic radii of Rb and NH 4 [11].In [6] it was shown that in the presence of Gaussian random fields with zero mean value, the proton-glass transition is smeared out, i.e., the susceptibility is rounded off, and the proton-glass order parameter remains finite at temperatures above the nominal freezing temperature (stability limit of the replica-symmetry solution).
Edwards-Anderson model [12] takes into account only the nearest neighbor interactions k ij = k.For k, the Gaussian distribution functions or the P (k) = (1−c)δ(k + 1) + cδ(k − 1) distribution are usually used (k = (−1, 1)).It is convenient to explore this model on a Bethe lattice, because the two-site cluster approximation for a free energy is exact for such lattices.In [13] an integral equation was obtained for the distribution function R(σ, 1) of a single effective field ϕ 1r exerted on the site 1 by the neighboring site r.In [14] analytical solutions were obtained (considering the first harmonics of Fourier transformation) at T = 0 and for the coordination numbers z = 4, 5, 6 and c = 0.5.These solutions contain δ-functions and a symmetric continuous part and correspond to the spin-glass state.The asymmetric solutions with a continuous part at T = 0 for z=3, h=0 and an arbitrary c(δ(σ)) found in [15] have the form state.It is necessary to note that for the Bethe lattice for a pure ferroelectric phase the function R(σ, 1) contains only one δ(σ − ϕ(T )) function, whereas in a purely ferroelectric phase ϕ(T ) = 0.Much later Mezard and Parisi [16] suggested the first step of replica symmetry breaking (1RSB) for the Bethe lattice.It was shown that earlier results [13][14][15] correspond to the replica symmetry (RS) solution.Recently [17] the EA model for the Bethe lattice with Gaussian distribution of interaction constants was studied numerically based on an algorithm suggested in [16].The difference between the free energy within the RS and 1RSB solutions is less then 4% for z=4,6 (figure 1).A similar conclusion for the k = (−1, 1) distribution was made earlier.The difference between RS and 1RSB solutions increases when z increases.We suppose that in order to describe the systems with competing short-range and long-range interactions, in which the prevailing role in forming the spin-glass state is played by short-range interaction, we can use the replica symmetry approach.We expect that such a situation takes place in a hydrogen bonded system of the Rb 1−x (NH 4 ) x H 2 PO 4 type (proton glass).In this compound the PO 4 tetrahedra and their environment play a decisive role in the formation of the energy levels of the systems.Our aim is to calculate within RS the thermodynamic characteristics and phase diagrams of Ising-like systems with an arbitrary radius of competing interaction.

Theory
Within the two-particle cluster approximation, the free energy F can be written as Hereafter we use the notations 3) The cluster fields ϕ ir i , ϕ jr j are exerted on the sites i, j by the sites r i , r j of different coordination spheres.They obey the equation Let us first consider a one-sublattice system (all sites are equivalent after the configurational averaging).We use the following distribution functions Here R n (σ, z n − 1) is the distribution function of z n − 1 cluster fields for the nth coordination sphere, R(σ, z n − 1) is the distribution function of all cluster fields except for the one field from the nth coordination sphere; z n is the coordination number of the nth coordination sphere; n = 1, . . ., M. We shall neglect a correlation between the cluster fields.Then, from (2.4)-(2.6)we can easily find that ( Rn (ζ, 1) is the Fourier transform of R n (σ, 1)) If M coordination spheres are taken into account, the relations (2.7), (2.8) yield M integral equations for M distribution functions R(σ, z n − 1).The free energy can be expressed in terms of the avarages with the corresponding distribution functions R(σ, z n − 1) and R(σ, z n ) = R(σ) Using (2.4) and (2.6), the function R(σ) can be expressed via R(σ, z n − 1) only: We solve the integral equation numerically in the case when only the nearest neighbor interactions are taken into account (the reference system with M = 1, z = z 1 ) with the distribution We solve this equation using the iteration method.As a trial function (zeroth approximation), the Gaussian function R (0) (σ, z − 1) is used with the parameters ϕ (0) , q (0) determined from the system (2.12) Here R (1) (σ, z − 1) is the first iteration expression (equation (2.11), where in the right hand side we put the R (0) (σ, z − 1) function).At low temperatures the functions R (i) (σ, z − 1), starting from the first iteration, contain δ−like peaks.For correct integration of such functions we use the following relation between the function f (x) and its Fourier-transform f (ζ): In figure 2 we show the results of the 5th iteration for the function R(σ, z − 1), the results obtained from the formulas given in [15], and the results obtained using the Gaussian trial function (zeroth iteration).Already the first iteration yields a qualitatively correct structure of the R(σ, z − 1) spectrum, consisting of δ−like peaks of finite widths at ±1, ±3.That is so due to the fact that at t → 0 (t = T /z) the solutions of the system of equations for ϕ 1r are close to ± 1. Hence, the sum of z − 1 cluster fields can have only the following values: ±(z − 1), ±(z − 3), . ... As one can see in figure 2, R(σ, z − 1) is a symmetric function of σ.Thence, at these values of the parameters (z = 4, c = 0.5) the system is in a spin glass state ( σ R(σ,z−1) = 0, σ 2 cum R(σ,z−1) = 0).δ-like peaks can be observed only at very low temperatures.Already at t ∼ 0.1 the shape of the R(σ, z − 1) curve is very close to the Gaussian function (see figure 3 for z = 6).At c = 0.95 the distribution function of cluster fields is a non-symmetric function of σ (figures 4,5).Thus, at these values of the parameters (z = 4, c = 0.95) the system is in a ferroelectric state ( σ R(σ,z−1) = 0, σ 2 cum R(σ,z−1) = 0).Using the Gaussian approximation for R(σ, z 1 − 1) we can obtain the well-known analytical results for the temperatures of the transitions from pure paraelectric phase to the spin glass phase (t g ) and from pure paraelectric to mixed ferroelectric phase (t c ) For the reference system we have shown that the iteration procedure for the integral equation for the distribution functions of (z − 1) cluster fields converges rapidly if the Gaussian form for the trial function is used (the 5th and high iterations do not differ).The free energy calculated with the Gaussian distribution function, which is determined from the extremum condition, and the free energy calculated with the distribution function determined from the integral equation, are close (see figure 6).It permits to use the Gaussian approximations for all distribution functions; parameters: Therefore, when M coordination spheres are taken into account, the 2M equations for 2M parameters can be obtained from the conditions of the free energy extremum.To further simplify the problem, we shall explore a simpler spin glass model with essential nearest neighbor interactions (the reference system) and weak interactions for other coordination spheres.
We consider a system consisting of two sublattices (figure 7) (f = 1 and f = 2 after configurational averaging; sites 2 are from the first coordination sphere for sites 1, and vice versa), and use a linear approximation for the interactions with all coordination spheres but the first one.It corresponds to the two-particle cluster approximation for the first coordination sphere and the mean field approximation for other coordination spheres.The free energy per unit cell F/N (at h = 0) can be written as 0 (σ+zϕ 1 +ϕ L,1 ) . (2.15) Here From the condition of F/N extremum we find the system of equations for the parameters ϕ f , ϕ L,f , q f .Hereafter we take into account only ferroelectric and antiferroelectric ordering.For these cases we can write that (2.18) The condition of free energy extremum yields the system of three equations for the three quantities: mean cluster field ϕ and dispersion q (pair cumulant) for the first coordination sphere, mean cluster field ϕ 0 for the other (for the long-range interactions).For instance, for the antiferroelectric ordering, the system of equations reads .

(2.19)
The expressions for the free energy, mean spin m 1 = −m 2 and mean square of spin

.21)
Dispersion for the long-range interactions is zero due to the linear approximation used.Therefore, these approximations can be used if the spin glass state is formed mainly by the short-range interactions.Consideration of the dispersion for the higher coordination spheres will be subject of our further studies.

Discussion
We construct the phase diagrams and perform a numerical analysis of the thermodynamic functions for the following types of the distribution function for the nearest neighbors interaction To compare the obtained results, we use the following relations between the Gaussian and the k = (−α, 1) distributions: For the proton glass model we have to take into account the random internal Gaussian fields with zero mean values and with the dispersion Q h .Now we have to replace σ i → σ i + h i in the expression (2.15) for the free energy and perform an additional averaging with the Gaussian fields h i .Solutions of the system of equations for ϕ 0 , ϕ, q are as follows: ϕ 0 = ϕ = 0, q = 0 (P-pure -pure paraelectric phase), existing only at Q h = 0, ϕ 0 = ϕ = 0, q = 0 (SG -spin glass phase), ϕ 0 , ϕ = 0, q = 0 (Fpure -pure ferroelectric phase, existing only at c = 1, or F-pure -pure ferroelectric phase, existing only at c = 0), ϕ 0 , ϕ = 0, q = 0 (F -ferroelectric or antiferroelelcltric phase, depending on the conditions (2.18)).
Figures 8 and 9 contain the phase diagrams for the cubic lattice for the symmetric k = (−1, 1) and antisymmetric k = (−0.5, 1) distributions, respectively, for the corresponding Gaussian distribution with non-random small positive long-range interaction J(J − = 2/3 • J). t g for the Gaussian distribution is lower, because in this case the values of k, which are close to zero, are more probable.An increase of non-random positive long-range interactions narrows the spin glass region.At large J the P-SG line disappears, the F-SG and AF-SG lines cross, and the SG phase exists only below the F-SG and AF-SG.Figures 10 and 11 illustrate the collapse of the SG region on the increase of the positive long-range interactions (figure 10) and the absence of the SG region for the quasi-one-dimensional model (figure 11).In figure 12 we show the temperature dependence of the mean spin value m, susceptibility χ, mean square of polarization Q = S i 2 TERM conf at c = 0.73 and c = 0.95 for different values of zJ.One can see that at c = 0.73, in contrast to the case of c = 0.95, small values of zJ essentially affect the system characteristics.The phase transitions between the paraelectric and ferroelectric phases and between the ferroelectric and spin glass phases at c = 0.73 (low temperature branch) are of the second order.The temperature dependences of thermodynamic characteristics of the system are shown in figure 13 for the symmetric Gaussian distribution and the distribution k = (−1, 1) at c = 0.5.We note that the Gaussian distribution not only lowers down the SG-transition, but also smears out the peaks of specific heat and susceptibility at the transition point.
Figures 14 and 15 contain the phase diagrams for the cubic lattice for an asymmetric distribution k = (−0.5, 1) and at different values of dispersion Q h of internal random fields at small non-random positive long-range interactions J(J − = 2/3 • J).In figures 16-19 we plot the temperature dependences of the mean square of polarization Q = S i 2 TERM conf and of susceptibility χ = ∂/∂h S i TERM conf for the cubic lattice and symmetric distribution k = (−1, 1) and for the Gaussian internal fields at zJ = 0.225 and c = 0.5 (transition from the high temperature SG-phase to the low-temperature SG-phase) and at c = 0.735 (transition from the high temperature SG-phase via the ferroelectric phase to the low-temperature SG-phase).Small internal fields essentially smear out the SG phase transition (figures 16,17).There is a transitional region between the high-temperature SG-phase (above the upper dashed line in figures 14, 15, q ≪ 1) and the low-temperature SG-phase (below the lowest dashed line in figures 14,15).The t -boundaries of the transitional region are determined as the upper and lower inflection points in the temperature curve of susceptibility (figure 17).At the same time, the internal random fields hardly affect the shape of the susceptibility peak at the transition between the SG and ferroelectric phase (figure 19).

Conclusions
We propose a spin glass model with an arbitrary range of competing interactions, considered within the two-particle cluster approximation for the free energy and within the symmetric replica approach.If M coordination spheres are taken into account, the free energy contains the distribution functions R(σ, z n − 1) (the distribution function of all cluster fields but one field from the nth coordination sphere), which order is not higher than M. For them, the system of M integral equations is derived within the non-correlated cluster fields approximation.In the case when only the first coordination sphere with the coordination number z=4,6 is taken into account, and the following distribution function of the short-range interactions P (k) = (1 − c)δ(k + 1) + cδ(k − 1) is used, the solutions for the R(σ, z − 1) function are studied numerically.At low temperatures it contains, in addition to the continuous part, the δ−like peaks of finite widths at ±(z − 1), ±(z − 3), . ... On the increase of temperature, only a single peak persists, with the envelope curve close to the Gaussian distribution.The free energy, calculated with the Gaussian distribution function determined from the extremum conditions, and the free energy calculated with the distribution functions determined from the integral equations, are close.This fact permits us to use the Gaussian approximations for the distribution functions, entering the expression for the free energy.
A proton glass model with essential short-range interactions (first coordination sphere with the distribution (−αk, k) or with the Gaussian distribution of interactions) and weak long-range interactions (linear approximation) is considered.We calculate the phase diagram, entropy, order parameter, specific heat, and susceptibility.Asymmetry of the phase diagram arises due to the long-range interactions taken into account and due to the asymmetry of the distribution function of short-range interactions.For this model, in the presence of non-random long-range interactions, we explore the effect of Gaussian fluctuations of internal fields with the dispersion ∼ x(1 − x) existing in hydrogen bonded systems of the Rb x (NH 4 ) 1−x H 2 PO 4 type.These fields are created by the structural disorder due to the difference between the ionic radii of Rb and NH 4 .Due to the internal fields, in the high-temperature range, instead of the paraelectric phase, the high-temperature proton glass exists (small but different from zero value of the spin glass parameter Q = S i 2 TERM conf .At the given value of x, on lowering the temperature, instead of the transition point from the paraelectric phase to the spin glass phase, a region of the transition from the high-temperature spin glass phase (small Q) to the low-temperature spin glass phase (large Q) emerges.It is shown that small fluctuations of internal fields can essentially smear out the region of such a phase transition.At the same time, the internal random fields hardly affect the shape of the susceptibility peak at the transition between the high temperature proton glass phase and the ferroelectric phase.
It should be noted that in proton glasses of the Rb x (NH 4 ) 1−x H 2 PO 4 type, the lower boundary of spin glass phase is experimentally determined from the peak in the temperature curve of the imaginary part of the low-frequency dielectric permittivity.Our next paper will be devoted to the calculation of frequency dependence of susceptibility of proton glasses with different types of competing interactions.

Figure 1 .
Figure 1.Energy and magnetization in the ground state as functions of µ = k [17].Simulations based on the 1 step of RSB are compared with solutions of the integral equation within the RS approximation (solid line).

Figure 7 .
Figure 7.Effective fields exerted on the sites i, j by the other sites from different coordination spheres (a lattice fragment before the configurational averaging is shown).We consider a model, which separates into two sublattices after averaging (× and • -sites on different sublattices).

Figure 12 .
Figure 12.Temperature dependence of the mean spin value m, susceptibility χ, mean square of polarization Q = S i 2 TERM conf at z = 6 for the distribution k = (−1, 1) at c = 0.73 (left column, SG region is to the left of the left peak, F-phase is in the central region) and c = 0.95 (right column, F-phase is to the left of the peak) for zJ = 0.0; 0.075; 0.15; 0.225 (dotted line).

Figure 13 .Figure 14 .
Figure 13.Temperature dependences of entropy S, specific heat C, mean square of polarization Q, and susceptibility χ for z = 6, c = 0.5 (SG phase) with the distribution k = (−1, 1) and with the Gaussian distribution for zJ = 0.0; 0.225 (in the SG phase only χ depends on J, t g for Gaussian distribution is lower.)

Figure 15 .Figure 16 . 2 TERM
Figure 15.The phase diagram at z = 6 for the distribution k = (−0.5, 1) and for the Gaussian distribution of internal fields with the dispersion Q h at zJ = 0.225.

Figure 17 .
Figure 17.Temperature dependence of susceptibility χ = ∂/∂h S i TERM conf at z = 6 for the distribution k = (−1, 1) and for the Gaussian distribution of internal fields with dispersion Q h at zJ = 0.225 and c = 0.5 (transition from the high-temperature SG-phase to the low-temperature SG-phase).

Figure 18 . 2 TERM
Figure 18.Temperature dependence of the mean square of polarization Q = S i 2 TERM conf at z = 6 for the distribution k = (−1, 1) and for the Gaussian distribution of internal fields with dispersion Q h at zJ = 0.225 and c = 0.735 (transition from the high-temperature SG-phase to the low-temperature SG-phase).

Figure 19 .
Figure 19.Temperature dependence of susceptibility χ = ∂/∂h S i TERM conf at z = 6 for the distribution k = (−1, 1) and for the Gaussian distribution of internal fields with dispersion Q h at zJ = 0.225 and c = 0.735 (transition from the high-temperature SG-phase to the low-temperature SG-phase).