Numerical study of the enlarged O(5) symmetry of the 3-D antiferromagnetic RP2 spin model

We investigate by means of Monte Carlo simulation and Finite-Size Scaling analysis the critical properties of the three dimensional O(5) non linear sigma model and of the antiferromagnetic RP2 model, both of them regularized on a lattice. High accuracy estimates are obtained for the critical exponents, universal dimensionless quantities and critical couplings. It is concluded that both models belong to the same Universality Class, provide that rather non standard identifications are made for the momentum-space propagator of the RP2 model. We have also investigated the phase diagram of the RP2 model extended by a second-neighbor interaction. A rich phase diagram is found, where most phase transitions are first order.


Introduction
Universality is sometimes expressed in a somehow defectively simple way: some critical properties (the universal ones) of a system are given by space dimensionality and the local properties (i.e. near the identity element) of the coset space G/H, where G is the symmetry group of the high temperature phase and H is the remaining symmetry group of the broken phase (low temperature). As we shall discuss, the subtle point making the above statement not straightforward to use, is that G needs not to be the symmetry group of the microscopic Hamiltonian, but that of the coarse-grained fixed-point action.
On the spirit of the above statement, some time ago [1] a seemingly complete classification was obtained of the Universality Classes of three dimensional systems where G=O (3). In this picture, a phase transition of a vector model, with O(3) global symmetry and with an O(2) low temperature phase symmetry, in three dimensions must belong to the O(3)/O(2) scheme of symmetry breaking (classical Heisenberg model). In addition, if H=O(1)=Z 2 is the remaining symmetry, the corresponding scheme should be O(4)/O(3) which is locally isomorphic to O(3)/O (1). 1 This classification has been challenged by the chiral models [2]. However, the situation is still hotly debated: some authors believe that the chiral transitions are weakly first-order [3], while others claim [4] that the Chiral Universality Class exists, implying the relevance of the global properties of G/H.
In this letter, we shall consider the three dimensional antiferromagnetic (AFM) RP 2 model [5,6,7,8], a model displaying a second-order phase transition and escaping from the previously expressed paradigm. It is worth recalling [9,10] that one of the phase transitions found in models for colossal magnetorresistance oxides [11] belongs to the Universality class of the AFM RP 2 model. The microscopic Hamiltonian of this model has a global O(3) symmetry group, while the low temperature phase has, at least, a remaining O(2) symmetry [9]. We will show here that the model belongs to the Universality Class of the three dimensional O(5) Non Linear σ model. Some ground for this arises from a hand-waving argument, suggested to us by one of the referees of Ref. [9] (see below).
The Universality Class of the three dimensional O(5) Non Linear σ model has received less attention that O(N) models with N = 0, 1, 2, 3 and 4. In spite of that, it has been recently argued that O(5) could be relevant for the High-temperature superconducting cuprates [12]. Nevertheless, perturbative field theoretic-methods have been used to estimate the critical exponents [13,14,15,16]. From the numerical side, only a rather unconvincing Monte Carlo simulation [17] was available until very recently. Fortunately, there has been a recent, much more careful study [18]. Yet, the scope of Ref. [18] was to determine whether an interaction explicitly degrading the O(5) symmetry to an O(3) ⊕ O(2) group was relevant in the Renormalization-Group sense. To that end, those authors concentrated in producing extremely accurate data on small lattices.
Our purpose is to study in greater detail the critical properties of the three dimensional O(5) Non Linear σ model, and of the AFM RP 2 model. We improve over previous studies of both models, obtaining more accurate estimates for critical exponents, universal dimensionless quantities and non universal critical couplings. As symmetries play such a prominent role, we will also explore the possibilities of changing those of the low-temperature phase by adding a second-neighbors coupling to the Hamiltonian of the AFM RP 2 model.

The models
We are considering a system of N-component normalized spins { v i } placed in a three dimensional simple cubic lattice of size L with periodic boundary conditions. The actions of our lattice systems are where the sums are extended to all pairs of nearest neighbors. Our sign convention is fixed by the partition function: d v being the rotationally invariant measure over the N-dimensional unit sphere.
To construct observables, in addition to the vector field v i , we consider the (traceless) tensorial field The interesting quantities related with the order parameters can be constructed in terms of the Fourier transforms of the fields For RP N −1 models, the local gauge invariance v i → − v i implies that the relevant observables are constructed in terms of the tensor field. However, for O(N) we have found very interesting to consider as well quantities related with the tensor field.
We construct the scalars (under global O(N) transformations and spatial translations) which, in addition to the action, are the only quantities measured during the simulation. Their mean values yield the propagators: In the thermodynamic limit and at the critical point, the propagator is expected to have poles at p f 0 = (0, 0, 0) and, for the antiferromagnetic model, at where ξ δp ≪ 1, and the exponents η f and η s correspond to independent wave function renormalization at each pole. Note that close to the critical point ξ f and ξ s are expected to remain proportional to each other (this will be explicitly checked numerically).
The (non-connected) susceptibilities are simply: In a finite lattice an extremely useful definition of the correlation length can be obtained from the (discrete) derivative of G(p). Using δp = (2π/L, 0, 0) one obtains [19,20] We also compute the cumulants 2 In the remaining part of this section, if a subindex V or T does not explicitly appears, it will imply that the equation is valid both for vector or tensor quantities. The same convention will apply for superscript (f) (ferromagnetic) and (s) (staggered). Staggered quantities are useful only for antiferromagnetic models.
Finally, the energy per link is We have computed β-derivatives of observables through their connected correlation with the action. Furthermore, we have extrapolated mean-values from the simulation coupling, to a neighboring value of β using the standard reweighting techniques, that cover all the relevant part of the critical region [20].
The relationship between the O(5) model and the RP 2 model arises from the Landau-Wilson-Fisher Hamiltonian for the RP 2 system [9]. Indeed, at the Mean Field level [9,6], the ferromagnetic quantities are simple functions of the staggered ones. This suggests to construct the Landau-Wilson-Fisher Hamiltonian from the staggered magnetization, which is a traceless, real, symmetric 3 × 3 matrix: Note that M s has 5 independent quantities. It is therefore a simple matter to obtain a five-components real vector v such that v 2 = tr M 2 s . The less trivial part regards the fourth-order interaction terms. In principle, the O(3) symmetry of the microscopic Hamiltonian would allow for a trM 4 s term and a [trM 2 s ] 2 . Surprisingly enough, both terms are proportional to ( v 2 ) 2 . Thus, assuming that sixth-order terms are irrelevant, the Landau-Wilson-Fisher Hamiltonian is expected to have a O(5) symmetry group and both models belong to the same Universality Class. This does not only implies that both models have the same critical exponents but also that the L → ∞ limit of U O(5) 4,V and of U s,RP 2 4 (evaluated at their respective critical couplings) coincide.

Numerical methods
In the O(5) model we have studied lattice sizes L = 6, 8, 12, 16, 24, 32, 48, 64, 96 and L = 128, at β = 1.1812. We have combined a Wolff's single cluster update with Metropolis. Our elementary Monte Carlo step (EMCS) consists of (10L+1) Wolff's cluster update and then a full-lattice Metropolis sweep. We take measurements after every EMCS. Since the average size of clusters grows as L 2−η ≈ L 2 , 80% of simulation time we are tracing clusters for all L, while Metropolis accounts for 10% of time and measurements for the remaining 10%. The total simulation time time has been the equivalent to 600 days of Pentium IV at 3.2 GHz. The number of EMCS ranges from 10 8 for L = 6 to 1.4 × 10 6 at L = 128. The integrated autocorrelation times for the susceptibility and for the energy are smaller than 1 EMCS for all the simulated lattice sizes.
Since we are interested in high accuracy estimates, we have used double precision arithmetics. One also needs to worry about the pseudo random-number generator. We have therefore implemented a Schwinger-Dyson test. It turned out that the 32-bits Parisi-Rapuano pseudo random-number generator [21] produces biased results. Either the Parisi-Rapuano plus congruential generator [22] or the 64 bits Parisi-Rapuano generator cured this bias. The 64 bits Parisi-Rapuano generator is faster and it has been our final choice.
For the antiferromagnetic RP 2 model, no efficient cluster method is available. We have simulated in lattice sizes from L = 8, 12, 16, 24, 32, 48 and L = 64 at β = −2.41. We used a multi-hit Metropolis sequential algorithm. Making a new spin proposal completely independent from the previous spin value, we achieve an acceptance of about 30%. We have used 2 hits what ensures a 50% acceptance. The observables have been measured every two Metropolis full lattice-sweep (our EMCS).
The number of EMCS ranges from 10 8 for L = 8 to 7 × 10 8 for L = 64. In units of the integrated autocorrelation time τ (for the order parameter) we have more than 10 6 τ for L = 64. The data up to L = 48 were obtained in Pentium IV clusters (simulation time was roughly equivalent to a 1000 days of a single processor). For the largest lattice, data were obtained in the Mare Nostrum computer of the Barcelona Supercomputing Center (simulation time was roughly equivalent to 3000 days of a single processor).
We perform a Finite-Size Scaling analysis, using the quotients method [5,6,20]. In this approach, one compares the mean value of an observable, O, in two systems of sizes L 1 and L 2 , at the value of β where the correlation-length in units of the lattice sizes coincides for both systems. If, for the infinite volume system, O (β) ∝ |β − β c | −x O , the basic equation of the quotient method is where the dots stand for higher order scaling corrections, ν is the correlationlength critical exponent, −ω is the (universal) first irrelevant critical exponent, while A O is a non universal amplitude. In a typical application, one fixes the ratio s = L 2 /L 1 to 2, and consider pairs of lattices L and 2L. A linear extrapolation in L −ω is used to extract the infinite volume limit. One just needs to make sure that the minimum lattice size included in the extrapolation is large enough to safely neglect the higher-order corrections. Of course, any quantity scaling like ξ at the critical point, such as LU 4 , may play the same role in Eq. (13). However, usually ξ yields smaller scaling corrections than U 4 .
The extrapolation method based on Eq. (13) is feasible for the antiferromagnetic RP 2 model. Unfortunately, for the O(5) model the amplitude A O is surprisingly small. In fact, resummation of the ǫ-expansion yields ω = 0.79(2) [13], while blind use of Eq.(13) on our numerical data would predict ω ≈ 2. We have then considered an additional correction term,Ã O L −σ . The exponent σ is an effective way of taking into account a variety of higher-order scaling corrections of similar magnitude (an L −2ω contribution, subleading universal irrelevant critical corrections, analytic corrections, effects of the non-linearity of the scaling fields, etc. [20]). Its utility will be in that it allow us to give sensible error estimates for the infinite-volume extrapolations, instead of bluntly taking A O = 0.
The most precise way of extracting the critical exponent, ω, and the critical point β c is to consider the crossing point of dimensionless quantities such as ξ/L and U 4 . Indeed, comparing their values in lattices L 1 and L 2 , one finds that they take a common value at The non universal amplitude B depends on the considered dimensionless quantity. Again, one usually take pairs of lattices L and 2L, and extrapolates to infinite volume using Eq. (14), maybe performing a joint fit for the crossing points of several dimensionless quantities. Again, for the O(5) model the amplitudes B are exceedingly small, and we need to add to Eq.(14) an analogous higher-order term, where σ plays the role of ω, and with amplitudeB.

Results for the O(5) model
The first step is the location of the critical point and the scaling corrections exponent. In table 1 we show the crossing points β L,2L c for the dimensionless quantities ξ V /L and U 4,V . To study the finite size corrections, we need to fit them to FixingB = 0 yields ω larger than 2 which is unacceptable given the field theory estimate ω = 0.79(2) [13]. Our interpretation is that B is too small to be observed even with our 6-digit accuracy. We have therefore fixed ω = 0.79(4) and we have taken as fit parameters β c , B,B and σ. We have doubled the field theory error in ω for safety. To further constrain β c (and σ), we have performed a joint fit of the crossing points for both ξ V /L and U 4,V , with the L β L,2L Notice that, for using equation (15), an estimate of ν is needed. Fortunately a rough estimate ν ≈ 0.78 (see below) is enough, given the uncertainty in ω and σ. As for the amplitudes of the leading correction term, we find while the amplitudesB are of order one. We then see that the L −σ term is crucial in order to obtain a sensible error estimate in the L → ∞ extrapolation.
At this point we may obtain two universal quantities, U * 4,V and ξ * V /L, namely the L → ∞ limit of U 4,V and ξ V /L evaluated exactly at the critical coupling. Again, due to the smallness of the leading scaling-corrections, we extrapolated to L → ∞ using the following functional forms: and Our numerical estimates for U 4,V (β L,2L c , L) and ξ V (β L,2L c , L)/L are displayed in the third and fourth columns of table 1. Although scaling-corrections are tiny, they can be clearly observed. Using Eqs. (18,19) we obtain while the amplitudes of the leading scaling-corrections are To obtain the critical exponents, we consider the operators ∂ β ξ V and χ V,T whose associated exponents are x ∂ β ξ V = ν + 1 and x χ V,T = γ V,T = ν(2 − η V,T ) . Taking the base 2 logarithm of the quotients, see Eq.(13), we obtain the effective size dependent exponents shown in table 2. In order to obtain their infinite volume value, we use (13), including an explicit L −σ term in the fit: (note that we have absorbed a constant factor 2 x O /ν into the amplitudes A for scaling-corrections). We obtain ν = 0.780(2), σ = 2.15 (19), A ∂ β ξ V = −0.04 (7) Less than a 5% of the total error is due to the error in ω = 0.79(4) for both ν and η V . For η T it is about a 30%. There are two points to be made about the extrapolation: • The O(5) model is not an improved action [23] (in the sense of exactly vanishing leading scaling-corrections), since A χ T is clearly non-zero (this is also illustrated in Fig. 1). Had we not considered the tensorial operators, this would have been completely missed. • It is somehow disappointing to compare the accuracy of the effective exponent ν in table 2 with the error in the extrapolation (23). Indeed, could we safely set A ∂ β ξ V = 0, the final result would have been ν = 0.7813(4).

Results for the antiferromagnetic RP 2 model
As we said before, qualitative arguments suggest that the antiferromagnetic RP 2 model belongs to the O(5) Universality Class. Our aim is to make the most astringent possible test of this hypothesis, thus we perform here an update of a previous study [6] of the RP 2 critical quantities. We report here largely improved estimates for critical coupling and exponents. Furthermore, we give estimates for the dimensionless quantity U s, * 4 , that can be directly compared with the U * 4,V obtained for the O(5) model. The ω used in the X axis correspond to the optimal value. Note that the minimum pair plotted is (12,24).
In this case, the extrapolation to the infinite volume limit is more standard (see for instance [24]) than for the O(5) model, because the amplitude of the leading scaling-corrections are much larger in most cases. To estimate ω and β c , we consider pair of lattices L and 2L, performing a joint fit to Eq.(14) of the crossing points of all four dimensionless quantities, imposing a common value of β c and ω (see Fig. 2). To control for systematic errors due to higher-order corrections, we follow the following procedure. We perform the fit using data for L ≥ L min , seeking a value of L min where a reasonable value of χ 2 /dof is found. Furthermore, we require that the fit performed for L > L min yield compatible results. In that case, we report the central value from the L ≥ L min fit, but taking the enlarged errors from the L > L min fit. We found that L min = 12 is enough for the extrapolation of β c . We obtain: (13) , ω = 0.78(4) , χ 2 /dof = 8.5/10 .
Once we have determined ω, we proceed to extrapolate U s, * 4 , and the critical exponents, using the analog of Eqs. (18,19,22) without the effective L −σ term. Although one could consider all four types of crossing points, β L,2L c , the resulting quotients would be highly correlated, making join fits scarcely useful. We concentrate on the crossing point of ξ s /L, which seems the most natural quantity, as we are dealing with an antiferromagnetic model and it can be obtained using only the two points correlation function. We have checked that other choices for β L,2L c yield compatible results, with slightly larger errors. Our extrapolations are shown together with the effective L dependent estimates in table 3. Error estimates in the extrapolation include the effect of the uncertainty in ω. For exponent ν, scaling corrections are completely buried in the statistical errors. We extrapolated with a simple linear fit, using L min = 8. The situation is rather different for η s . For that exponent, enlarging L min systematically increases the asymptotic estimate. On the other hand, a fit quadratic in L −ω yields a linear term compatible with zero. The linear extrapolation with L min = 16 is identical to the quadratic extrapolation from L min = 8. This is the result indicated in table 3. As for η f , we have rather strong leading scaling corrections. Indeed, a fit linear in L −ω yields basically identical results for L min = 8 and L min = 12 (this is the result reported in table 3). Furthermore, a fit quadratic in L −ω including all points, yielded η f = 1.331 (5). The extrapolation for U s, * 4 is equally simple.  Table 3 Size

Next nearest neighbors coupling
A rather subtle question regards the symmetry of the low-temperature RP 2 antiferromagnetic phase [6,9]. A way of investigating this problem is to study the enlarged action where an additional second-neighbors coupling is considered.
The phase diagram for β 1 < 0 (Fig. 3) contains the following regions (spins are classified as even or odd, according to the parity of x i + y i + z i ): • PM: the usual (paramagnetic) disordered state, where the O(3) symmetry of the action (26) is preserved. • O(2): (say) even spins fluctuate almost parallel to (say) the Z axis, with random sense (local Z 2 symmetry), while odd spins fluctuate in the perpendicular plane (global O(2) symmetry) • O(1): two sublattices with ferromagnetic ordering in perpendicular directions, with random sense (the local Z 2 symmetry, v i → − v i is always preserved). • Skyrmion/Flux: the spins are parallel to the diagonals of the unit cube, so that they point out from/to the center (i.e. the propagator show three peaks at p 0 = (π, π, 0) and permutations). It is interesting to note the vectorial version of this phase appear in models for colossal magnetorresistance oxides [25].
The most relevant results can be summarized as follows: • We have obtained critical exponents for several points along the PM-O(2) critical lines, with significantly less accuracy than for the β 2 = 0 model. No variation was observed within errors. • The O(2)-O(1) critical line is repelled from the β 2 = 0 axis by the secondneighbors coupling. We interprete this as a competition with the order-from-disorder mechanism [6,9] behind the PM-O(2) transition. • A naive analysis suggests that the O(2)-O(1) transition line should belong to the XY universality class (consider first the limit β 1 = −∞, then the identity cos 2 θ = (1 + cos 2θ)/2 for the less ordered face-centered cubic sublattice). However, we have found that this transition line is first order, as revealed by the double-peaked histogram of the second neighbors energy. We are able to estimate a non-zero latent-heat up to β 1 = −6. At β 1 = −4 the double-peak structure is still easy to observe on small lattices. We presume that the whole line is first-order, although it could become very weak. • The skyrmion-PM transition lines turned out to be first-order at all the checked points. • Note that at β 1 = 0, we have two decoupled ferromagnetic RP 2 models on the face-centered cubic lattice (a model showing first-order transition, well known in the liquids-crystal context [26]). We should remark that a precise location of the triple point O(2)-O(1)-PM is very difficult to achieve.

Conclusions
We have obtained high accuracy estimates of critical exponents and other universal quantities for the three-dimensional O(5) and the antiferromagnetic RP 2 models, by means of Monte Carlo simulation, Finite-Size Scaling analysis and careful infinite-volume extrapolation.
In the case of the O(5) model the coupling to the leading irrelevant operator is rather weak, but non-vanishing, and one needs to consider higher-order scaling corrections to obtain sensible error estimates. In spite of that, our estimates for the critical coupling and the anomalous dimensions for the vector and tensor representations improve significantly over previous work (see table 4).
For the RP 2 model, the leading scaling corrections are sizeable. It is amusing that we are able to obtain numerically (for the first time, we believe) an estimate for the (universal) scaling-correction exponent ω = 0.78(4), of accuracy comparable to the perturbative field-theoretical estimate [13] ω = 0.79 (2). As table 4 shows, within the achieved accuracy, both models seem to belong to the same Universality Class (for comparison, we also show results from the O(4) Universality Class). To conclude this, one needs to accept that in the RP 2 model, the wavefunction renormalization for the propagator pole at p 0 = (π, π, π) is as for the the O(5) fundamental field, while at p 0 = (0, 0, 0) is as for the O(5) tensor field.
We have also obtained the phase-diagram of the RP 2 model extended with a second nearest-neighbors interaction. We have found a rich phase diagram.  Table 4 Summary of infinite-volume estimates for the 3D antiferromagnetic RP 2 , O(5) and O(4) models. We call η ′ to η f for RP 2 and to the η T for O(N ) models. FT stands for Field-Theory.