Dynamical generation of a gauge symmetry in the Double-Exchange model

It is shown that a bosonic formulation of the double-exchange model, one of the classical models for magnetism, generates dynamically a gauge-invariant phase in a finite region of the phase diagram. We use analytical methods, Monte Carlo simulations and Finite-Size Scaling analysis. We study the transition line between that region and the paramagnetic phase. The numerical results show that this transition line belongs to the Universality Class of the Antiferromagnetic RP(2) model. The fact that one can define a Universality Class for the Antiferromagnetic RP(2) model, different from the one of the O(N) models, is puzzling and somehow contradicts naive expectations about Universality.


Introduction
The problem of the dynamical restoration of a gauge symmetry (see e.g. [1,2] and references therein) has received considerable attention in the recent 10 years, because of the problem of introducing a chiral gauge theory in the lattice. Although the Ginsparg-Wilson [3] method has somehow superseded this approach, the question remains an interesting one. Naively, the problem could seem a trivial one [2]: the non gauge-invariant terms in the action generate a high-temperature-like gauge-invariant expansion with a finite radius of convergence. The subtle question is whether the radius of convergence of this expansion will remain finite in the continuum limit or not. In this letter, we want to address a related, but different question, namely, the generation of a local invariance in the low-temperature (broken-symmetry) phase of a system without any explicitly gauge-invariant term in the action. We have found this intriguing phenomenon in a numerical study of a simplified version of the Double-Exchange Model [4,5](DEM), one of the most general models for magnetism in condensed matter physics, still under active investigation [6]. The local invariance does not follow from a high-temperature-like expansion, but from the infinite degeneracy of the ground-state (see next section), which occurs at a unique value of the control parameter at zero temperature, then extending to a finite region of the phase diagram at finite temperature. This phenomenon reminds one of the so-called Quantum-Critical Point phenomenology [7].
We have studied the model using Monte Carlo simulations and Finite-Size Scaling techniques [8,9,10]. We have found that the critical exponents are fully compatible with the ones [11] of the antiferromagnetic (AFM) RP 2 model in three dimensions [11,12,13], which has an explicit local Z 2 invariance. This might not be surprising, given the strong similarities in the ground-state of both models (see next section). However, the fact that one can explicitly show that there is a Universality Class associated to the AFM-RP 2 model is puzzling. Indeed, the most ambitious formulation of the Universality Hypothesis states that the critical properties of a system are fully determined by the space dimensionality and its symmetry groups at high-temperature (G) and low-temperature (H). Moreover, systems with a locally isomorphic G/H are expected to have the same critical behavior. In our case, G=O(3), and from our numerical study H seems to be O(2) [14], although the O(2) residual symmetry could be also broken (to O(1)=Z 2 ). In the former case, the Universality Class should be the one of the O(3) non-linear σ model, while in the latter case one expects O(4) non-linear σ model-like behavior [15]. In recent years, the Universality Hypothesis (as stated above) has been challenged in a number of frustrated, chiral models [16]. Yet, detailed numerical simulations have shown that typical transitions are weakly first-order [17], which is hardly surprising, because the typical critical exponents proposed for chiral systems [16] are fairly similar to the effective exponents one expects in weak first-order transitions [18]. On the other hand, we have no doubt that the transition here studied is continuous, but we have no alternative explanation for our results.

The Model
Although some powerful techniques have been developed [19] for the Double-Exchange model [4] (involving dynamical fermions), lattice sizes beyond L = 16 are extremely hard to study with the present generations of computers. Thus, one may turn to the simplified version proposed by Anderson [4], where a purely bosonic Hamiltonian is considered. This simplified model has been recently studied [5] by extensive Monte Carlo simulation. Yet, previous studies missed several phases in the phase diagram (see Fig. 1 and below).
Specifically, we consider a three dimensional cubic lattice of side L, with periodic boundary conditions. The dynamical variables, φ i , live on the sites of the lattice and are three-component vectors of unit modulus. The Hamiltonian contains the Anderson version of the Double-Exchange model plus an antiferromagnetic first-neighbor Heisenberg interaction [20]: where < i, j > means first-neighbor sites on the lattice. The partition function reads the integration measure being the standard rotationally-invariant measure on the sphere. In the following, expectation values will be indicated as . . . .
The zero temperature limit is dominated by the spin configurations that minimize the energy. Exploratory Monte Carlo simulations showed that these configurations have a bipartite structure. Indeed, let us call a lattice site even or odd, according to the parity of the sum of its coordinates ( r i = (x i , y i , z i ), x i + y i + z i even or odd). Then the spins on the (say) even lattice are all parallel, while the spins on the odd sublattice are randomly placed on a cone of angle Ω ( φ i · φ j = cos Ω) around the direction of the even lattice. The energy when T goes to zero is simply Now, Ω(J) is obtained by minimizing H 0 (Ω). For J > − √ 2/4, Ω = 0, mean-ing a ferromagnetic vacuum. For − √ 2/4 > J > −0.5, one has 0 < Ω < π/2 (ferrimagnetic vacuum), while for all J < −0.5, it is π/2 < Ω < π (antiferrimagnetic vacuum). The full antiferromagnetic configuration (Ω = π) is never stable at zero temperature. The J = −0.5 point is very peculiar: much like for the AFM-RP 2 model [13,12,11], spins in the even sublattice are randomly aligned or anti-aligned with the (say) Z axis, while the spins in the odd sublattice are randomly placed in the X, Y plane. Since spins in the two sublattices are perpendicular to each other, one can arbitrarily reverse every spin. A local Z 2 symmetry is dynamically generated and, as we will see, it extends to finite temperatures. An operational definition of dynamical generation of a gauge symmetry is the following. One must calculate the correlation-length for non gauge-invariant operators (ξ NGI ) and compare it to the correlation-length corresponding to gauge-invariant quantities (ξ GI ). In the continuum-limit (ξ GI → ∞), one should have ξ NGI /ξ GI → 0. We have checked that the correlation-length associated to the spin-spin correlation function (non Z 2 gauge-invariant) is smaller than 0.3 for all temperatures at J = −0.5. The alert reader will notice that the symmetry group at the point (T, J) = (0, −0.5) is rather a local O(2), besides the local Z 2 previously discussed. However, the associated correlation-length at finite temperature grows enormously when approaching the critical temperature (tenths of lattice-spacings already at T = 0.9T c ), and probably diverges. More details on this will be given in Ref. [14].
A further analytical evidence for this fact can be obtained by performing a Taylor expansion of the action at J = −1/2 assuming that φ i · φ j , which just vanishes at J = −1/2 and zero temperature, is small: We can assume that this expansion has a finite radius of convergence and so we can extend this series to the non zero temperature region. Notice that the first term in the expansion is just the AFM-RP 2 Hamiltonian modified by terms that are no longer gauge invariant (those with odd powers in the scalar product). We can argue that those terms, which break explicitly the gauge invariance, are irrelevant operators, in the Renormalization Group sense, at the PM to AFM-RP 2 critical point and so, our model at finite temperature should belong to the same Universality Class as the AFM-RP 2 model. Obviously, were the transition of the first order, the argument would not apply.
In the Z 2 gauge-invariant phase, the vectorial magnetizations defined as are zero. Thus, we define proper order parameters, invariant under that gauge symmetry, in terms of the spin field, φ i and the related spin-2 tensor field (which is Z 2 invariant): Then we define: The different phases we find (see [14] for  [7]. A detailed analysis of this phase diagram will appear soon [14].
Let us end this section by defining the observables actually used in the simulation. They are obtained in terms of the Fourier transform of the tensor fields: and their propagators where p = 2π L n , n i = 0, . . . , L − 1 .
Then we have the usual (χ u ) and staggered (χ s ) susceptibilities, Having those two order parameters, we must expect the following behavior for the propagators in the thermodynamic limit, in the scaling region and for T > T c : where ξ u and ξ s are correlation-lengths, t = (T −T c )/T c is the reduced temperature and Z u,s , A u,s and B u,s are constants. On the other hand, the anomalous dimensions η u and η s need not be equal: we can relate them to the dimensions of the composite operators following the standard way: d − 2 + η u = 2dim(τ u ) and d − 2 + η s = 2dim(τ s ), and in general dim(τ s ) = dim(τ u ).
With those quantities in hand, one can define a finite-lattice correlationlength [21] for the staggered, and non staggered sectors: Other quantities of interest are the dimensionless cumulants Besides the above quantities, we also measure the energy (1), which is used in a reweighting method [22], and temperature derivatives of operators through their crossed correlation with the energy.

Critical behavior
For an operator O that diverges as |t| −x O , its mean value at temperature T in a size L lattice can be written, in the critical region, assuming the finite-size scaling ansatz as [8] O where F O is a smooth scaling function and ω is the universal leading correctionto-scaling exponent.
In order to eliminate the unknown F O function, we use the method of quotients [9,10,11,23]. One studies the behavior of the operator of interest in two lattice sizes, L and rL (typically r = 2): Then one chooses a value of the reduced temperature t, such that the correlationlength in units of the lattice size is the same in both lattices [11]. One readily obtains Notice that the matching condition Q ξ = r can be easily tuned with a reweighting method. The usual procedure consists on fixing r = 2, and obtaining the  Table 1 Results of the infinite volume extrapolation with Eq. (19) to obtain T c and ω. Q is the quality-of-fit parameter. Our final values are the bold values.
above quotients for several L values in order to perform an infinite volume extrapolation.
In order to obtain the critical exponents, we use as operators χ u,s (x χ = γ = 2 − η), ∂ T ξ s (x ∂ T ξs = ν + 1). Notice that several quantities can play the role of the correlation-length in Eq. (18): ξ u , ξ s , Lκ u and Lκ s . This simply changes the amplitude of the scaling corrections, which will turn out to be quite useful.
Another interesting quantity is the shift of the apparent critical point (i.e. rξ(L, T L,r c ) = ξ(rL, T L,r c )), with respect to the real critical point:

The Simulation
We have studied the model (1)  We have carried out 20 million full-lattice sweeps (measuring every 2 sweeps) at each lattice size at T = 0.056. For the L = 64 lattice we have also performed 20 million sweeps at T = 0.0558. The largest autocorrelation time measured is about 1400 sweeps (corresponding to χ s ). To ensure the thermalization we have discarded a minimum of 150 times the autocorrelation time.
The computation was made on the RTN3 cluster of 28 1.9GHz PentiumIV processors at the University of Zaragoza and the total simulation time was equivalent to 11 months of a single processor.

Critical exponents
The first step, as usual, is to estimate ω from Eq. (19). For this, one needs a rough-estimate of ν. Since our data for the quotient of ∂ T ξ s using ξ s as correlation-length show very small scaling corrections (see Fig. 3), one can temporally choose ν = 0.78 and proceed with the determination of ω. Having four possible dimensionless quantities, ξ u /L, ξ s /L, κ u and κ s , we can perform a joint fit to Eq. (19) constrained to yield the same T c . This largely improves the accuracy of the final estimate. The full covariance matrix is used in the fit, and errors are determined by the increase of χ 2 by one-unit. Our results are summarized in Fig. 2 and table 1. Although scaling corrections are clearly visible, good fits are obtained from L min = 12. Thus, we conclude that ω = 0.82 (5) , T c = 0.055895 (5) . We are now ready for the infinite volume extrapolation of the critical exponents. As usual, one needs to worry about higher-order scaling corrections.
Here we shall follow a conservative criterion: we shall perform the fit to Eq. (18) only for L ≥ L min , and observe what happens varying L min . Once we found a L min for which the fit is acceptable and the infinite volume extrapolation for L ≥ L min and L > L min are compatible within errors, we take the extrapolated value from the L ≥ L min fit, and the error from the L > L min fit. Our results can be found in Fig. 3 and in table 2. As well as for the critical temperature, we used all four dimensionless quantities in a single fit constrained to yield a common infinite volume extrapolation. Our final estimates are ν = 0.781 (18)(1), η s = 0.032(5)(2), η u = 1.337(6)(7) , were the second error is due to the uncertainty in ω.
One can compare these results with other models, once it is decided what is going to play the role of our η u in the O(N) models. Our candidate is the  Table 2 Infinite volume extrapolation for the critical exponents, using Eq. (18)

Conclusions
We have studied numerically a bosonic version of the DEM, Eq. (1), by Monte Carlo simulations, obtaining its full phase-diagram, missed in previous stud- 1 In O(N ) models it is possible to compute η φ⊗ φ using field theoretical methods.
ies [5]. We have studied its critical behavior with Finite-Size Scaling techniques. As Eq. (21) and Eq. (22) show, our results for the critical exponents are fully compatible with the results for the AFM-RP 2 model, barely compatible with the O(4) model, and fully incompatible with the results for the O(3) model. Our results in the low temperature phase [14] seem to indicate that the scheme of symmetry-breaking is O(3)/O(2), which contradicts Universality. Most puzzling is the excellent agreement between the present results and the estimates for the AFM-RP 2 model. This seems to indicate that the AFM-RP 2 model really represents a new Universality Class in three dimensions, in plain contradiction with the Universality-Hypothesis (at least in its more general form). This seems to imply that the local isomorphism of G/H is not enough to guarantee a common Universality Class. Of course, it could happen that we were seeing only effective exponents and that in the L → ∞ limit a more standard picture arises. Yet, we do not find any obvious reason for two very different models to have such a similar effective exponents.
Another intriguing effect is that the augmented local Z 2 symmetry of the point (T = 0, J = −0.5) extends to the finite temperature plane, which is recalling a Quantum Critical Point [7]. Indeed, we find that the Universality Class is the one of an explicitly gauge-invariant model. To our knowledge this is a new effect in (classical) Statistical Mechanics, and deserves to be called dynamical generation of a gauge symmetry.