text:
# Electron Glass in Ultrathin Granular Al Films at Low Temperatures

## Abstract

Quench-condensed granular Al films, with normal-state sheet resistance close to 10 k$`\mathrm{\Omega }`$/$`\mathrm{}`$, display strong hysteresis and ultraslow, non-exponential relaxation in the resistance when temperature is varied below 300 mK. The hysteresis is nonlinear and can be suppressed by a dc bias voltage. The relaxation time does not obey the Arrhenius form, indicating the existence of a broad distribution of low energy barriers. Furthermore, large resistance fluctuations, having a $`1/f`$-type power spectrum with a low-frequency cut-off, are observed at low temperatures. With decreasing temperature, the amplitude of the fluctuation increases and the cut-off frequency decreases. These observations combine to provide a coherent picture that there exists a new glassy electron state in ultrathin granular Al films, with a growing correlation length at low temperatures.

The electron glass was first predicted to exist in disordered interacting systems nearly two decades ago . Such nonergodic behavior is very interesting because one normally expects electron systems to relax rather rapidly. Over the years, a number of studies have been reported in which electrons display glassy dynamics that are often associated with non-exponential relaxation extending over many decades in time, such as the field-effect conductance measurements in compensated GaAs , amorphous indium-oxide films , and ultrathin Bi/Ge and Pb/Ge films . The glassy behavior is believed to arise from the electron-electron (e-e) interactions and the Coulomb gap . Recently, the electron glass has received renewed interest  as the subject of e-e interactions has become a central topic in understanding the metal-insulator transition in two-dimensions . However, the precise role that the Coulomb gap plays in the observed glassy behavior is not clear since there is no simultaneous measurement of the conductance relaxation and the single particle density of states. In addition, there is no direct experimental effort to probe the correlation length in the glassy phase.

In this Letter we report glassy behavior in the normal state of quench-condensed weakly insulating granular Al films of sheet resistance, $`R_{\mathrm{}}`$, of about 10 k$`\mathrm{\Omega }`$/$`\mathrm{}`$ at 300 mK. We measure the relaxation of $`R_{\mathrm{}}`$ after the temperature was varied. We have focused on weakly insulating films because $`R_{\mathrm{}}`$ on the order of 10 k$`\mathrm{\Omega }`$/$`\mathrm{}`$ is easy to measure using sensitive ac lock-in techniques, which turn out to be crucial in measuring the resistance fluctuations described below. We observed that, below 300 mK, the resistance was strongly hysteretic and displayed ultraslow, non-exponential relaxation as the temperature was varied. We have also observed strong nonlinear behavior in the hysteretic regime. What was unique to our work was the first observation in quench-condensed metal films of large resistance fluctuations below 100 mK. We argue that these observations indicate the existence of a glassy electron state in ultrathin granular Al films with a growing correlation length at low temperatures.

Our Al films were quench-condensed onto glass substrates using 99.999$`\%`$ purity Al sources in UHV conditions inside a dilution refrigerator, with the substrates being held near 20 K during evaporation. The film thickness was near 25 Å. After warming up to room temperature, the films showed a granular morphology with a typical grain size of about 300 Å, as seen in scanning force microscopy studies. The films had a multi-lead pattern with an area of 3$`\times `$3 mm2 between the neighboring leads. An analog lock-in amplifier, operating at 27 Hz with a time constant of 3 seconds, was used to measure the four-probe ac resistance. The ac probe current was fixed at 1 nA, producing an ac voltage of about 10 $`\mu `$V across a film section. For nonlinear studies described below, an additional dc bias voltage was applied. All the data described below were measured in the normal state, with superconductivity being suppressed by a magnetic field above the spin-paramagnetic limit  of about 48 kG. To date, we have performed detailed studies on films of $`R_{\mathrm{}}`$ $``$ 10 k$`\mathrm{\Omega }`$/$`\mathrm{}`$. Such films appeared to be very uniform, with $`R_{\mathrm{}}`$ varying less than 5$`\%`$ among the various sections of the multi-lead pattern. These films are far above the percolation threshold since $`R_{\mathrm{}}`$ scales with inverse film thickness.

In the inset to Fig. 1, we show a typical resistance hysteresis measured on one Al film as temperature was cycled between 500 and 30 mK, in the high-field normal state with a field, $`H_{}`$ $`=`$ 52 kG, applied perpendicular to the film plane. Throughout all our experiments, the rates of cooling and heating were kept constant on a logarithmic temperature scale, and the cooling and the heating cycles each took 1.5 hr. We can exclude temperature lag or heating as the cause of the hysteresis, because we have not observed any hystersis in Li films of $`R_{\mathrm{}}`$ $``$ 10 k$`\mathrm{\Omega }`$/$`\mathrm{}`$ and Be film of $`R_{\mathrm{}}`$ $`=`$ 3 $``$ 500 k$`\mathrm{\Omega }`$/$`\mathrm{}`$, both quench-condensed on the same glass substrates. For relaxation studies, the films were first cooled from 500 mK to the desired temperature, $`T`$, using the above logarithmic rate, at which $`R_{\mathrm{}}`$ was measured as a function of time. In Fig. 2, we show the $`R_{\mathrm{}}`$ vs. time curves for a number of fixed $`T`$ in a field of $`H_{}`$ $`=`$ 52 kG. We found that the curves at higher $`T`$, such as 100 mK, could be fit very well to a stretched-exponential form, $`[R_{\mathrm{}}(\mathrm{})R_{\mathrm{}}(t)]/[R_{\mathrm{}}(\mathrm{})R_{\mathrm{}}(0)]=exp[(t/\tau )^\gamma ]`$, over three decades in time. With decreasing $`T`$, the time constant $`\tau `$ increased sharply. The exponent $`\gamma `$ scattered between 0.6 and 0.8 at low $`T`$ without a clear trend. We note that below 100 mK, the large resistance fluctuations seen in Fig. 2 made it difficult to accurately determine $`\tau `$ and $`\gamma `$. In the main part of Fig. 1, we plot $`\tau `$ vs. $`1/T`$ for two field values in both perpendicular and parallel field orientations. Data in Fig.1 appear to fall into two groups: one at 52 kG and the other at 90 kG, suggesting that $`\tau `$ was larger in higher fields and was insensitive to field orientation. The almost linear dependence of $`\tau `$ on $`1/T`$ in Fig. 1 indicates that, with decreasing $`T`$, $`\tau `$ increased much slower than the Arrhenius law, $`\tau (T)=\tau _0exp(E_a/T)`$, where $`E_a`$ is known as the activation energy describing the typical energy barrier to relaxation . Such non-Arrhenius behavior can be explained only if there exists a very broad distribution of low-energy barriers. Relaxation occurs over lower energy barriers with decreasing $`T`$, leading to a relaxation time that increases slower than the Arrhenius law.

The above conclusion is also supported by the nonlinear effects in the hysteresis. For nonlinear studies, we have measured $`R_{\mathrm{}}`$ in the presence of a dc bias voltage, $`V_{bias}`$, across the films, and have observed that the hysteresis can be suppressed if $`V_{bias}`$ is large enough. The experimentally relevant $`V_{bias}`$ was of the order of 0.5 mV. It leads to a bias current of 0.05 $`\mu `$A and a Joule heating of 2.5$`\times `$10-11 W. We believe that the nonlinear effects were not caused by heating. Consider the hysteresis shown in the inset to Fig. 1, we observed that, between 100 mK and 200 mK, the heating curve shifted downward and the cooling curve shifted upward with increasing $`V_{bias}`$. It is obvious that had heating by the bias current been significant the cooling curve should have also shifted downward with increasing $`V_{bias}`$, since the resistance of weakly insulating films should decrease with increasing $`T`$. In Fig. 3(a), we plot the width of the hysteresis, $`\mathrm{\Delta }R`$, versus $`T`$ for various $`V_{bias}`$ values. Such data are obtained by subtracting the cooling curve from the heating curve in a thermal cycle such as the one shown in the inset to Fig. 1. It is quite revealing in that the bias voltage does not uniformly suppress the hysteresis across the entire range of $`T`$. Instead, a small $`V_{bias}`$ suppresses the low-$`T`$ components of the hysteresis only. With increasing $`V_{bias}`$, progressively higher-$`T`$ components of the hysteresis are suppressed. We believe that relaxation processes over low-energy barriers are suppressed by the bias voltage in the same way the low-$`T`$ components of the hysteresis are suppressed. Nevertheless, we have yet to develop quantitative methods to analyze the data in Fig. 3(a) to obtain the distribution of low-energy barriers. We note that the falling off of $`\mathrm{\Delta }R`$ at low $`T`$ in the $`V_{bias}`$ $`=`$ 0 curve in Fig. 3(a) was an experimental artifact because we did not wait at 30 mK for the system to relax after cooling down from 500 mK. Had we waited long enough before heating up towards 500 mK, $`\mathrm{\Delta }R`$ would have been the largest at 30 mK for the $`V_{bias}`$ $`=`$ 0 curve.

If we make a vertical cut at a fixed $`T`$, such as 175 mK, in Fig. 3(a), we obtain $`\mathrm{\Delta }R`$ as a function of $`V_{bias}`$. We have observed that the absolute value of $`\mathrm{\Delta }R`$ varied with the strength and the orientation of the magnetic field. However, if we scale $`\mathrm{\Delta }R`$ by its value in the $`V_{bias}`$ $`=`$ 0 limit, $`\mathrm{\Delta }R(0)`$, the resulting scaled width, $`\mathrm{\Delta }R`$/$`\mathrm{\Delta }R(0)`$, appears to be independent of the strength and the orientation of the field. This is demonstrated by the good data collapse in Fig. 3(b) in which we plot four sets of such scaled data at 175 mK as a function of $`V_{bias}`$ for two field values in both perpendicular and parallel orientations. In this semi-log plot, $`\mathrm{\Delta }R`$ falls off steeply for $`V_{bias}`$ $`>`$ 0.7 mV. This suggests that there is a field-independent threshold bias, $`V_0`$ $``$ 0.7 mV, that suppresses the hysteresis. By making vertical cuts at various temperatures in Fig. 3(a), we observe that $`V_0`$ decreases with decreasing $`T`$. This is consistent with the picture that relaxation occurs over lower energy barriers at lower $`T`$; as a result, the value of $`V_0`$ for suppressing the hysteresis is also lowered. We plot $`V_0`$ vs. $`T`$ in Fig. 3(c). The flattening of the data below 100 mK is again due to the experimental artifact discussed earlier. For the high-$`T`$ part of the data, $`V_0`$ can be fit to a power law temperature dependence, $`V_0`$ $``$ $`T^{3.1\pm 0.3}`$, as shown by the solid line in the inset to Fig. 3.

The results discussed above suggest the existence of a glassy state in our quench-condensed ultrathin granular Al films. The low temperature scale of our experiments excludes the possibility of glassy structural relaxation, such as atomic diffusion between metastable configurations similar to those observed in a-Si:H . Tunneling studies have shown in such Al films , as well as in other films of similar sheet resistance , that an anomaly in the density of state exists at the Fermi energy. However, a true Coulomb gap has not been found in such weakly insulating films. Thus it is not clear what role, if any, the Coulomb gap plays in the observed glassy relaxation. In addition, we did not find any hysteresis in quench-condensed Li films of $`R_{\mathrm{}}`$ $``$ 10 k$`\mathrm{\Omega }`$/$`\mathrm{}`$ and Be films of $`R_{\mathrm{}}`$ $`=`$ 3 $``$ 500 k$`\mathrm{\Omega }`$/$`\mathrm{}`$, even though a true Coulomb gap has been observed in Be films . While we were not able to investigate the morphology of the Li films, which become unstable in air, scanning force microscopy studies of the Be films warmed up to room temperature did not find any granular structure down to 1 nm. Because of the morphological difference, we propose that granularity plays an important role in the glassy behavior of our Al films.

Consider a simple conduction model for granular metals in which an electron moves from one neutral grain to another nearby one to create a charge-anticharge pair . The energy cost for creating such a pair is 2$`E_c`$, where $`E_c=e^2/2C`$ is the grain charging energy, with C being the capacitance of a grain. $`E_c`$ can be significant if the grains are small, and can vary strongly with location due to the randomness in grain sizes. This results in a rough potential background for the charge carriers. Transport is mediated via the ionization of such pairs , with the energy associated with ionization being the Coulomb attraction between the pair, which, in two dimensions, has a logarithmic form with a cutoff length . This model has motivated experiments searching for a finite-$`T`$ Kosterlitz-Thouless-Berezinskii charge unbinding transition in arrays of Josephson junctions  and granular films . However, in a recent theoretical study of capacitively coupled grains, Granato and Kosterlitz  have found that finite-$`T`$ transitions are suppressed by disorder. They found instead a charge glass with a correlation length that diverges as temperature decreases to zero. They have also predicted that nonlinear behavior sets in at a characteristic voltage $`V_cT^{1+n}`$, with a thermal critical exponent $`n`$ $``$ 1.7 characterizing the $`T`$ $`=`$ 0 charge glass transition. It is interesting to note the closeness of this predicted power-law dependence with the power-law fit in Fig. 3(c). However, we caution that our fit has a limited range. Within this picture of granular transport, we can qualitatively understand the behavior of the films in an applied magnetic field as shown in the main part of Fig. 1. The field splits the singlet states, effectively reducing the density of states. As a result, the relaxation time increases with increasing magnetic field, however, it is insensitive to the orientation of the field.

We note that the data in Fig. 2 show drastic increases in the size and decreases in the characteristic frequency of the resistance fluctuations with decreasing $`T`$, as shown by the 60 mK and 30 mK curves. This is a unique feature in our data that has not been reported in similar quench-condensed films. We have analyzed such noise data using a fast Fourier transform (FFT). We first obtain the noise time traces by subtracting the stretched-exponential fitting curves from the $`R_{\mathrm{}}`$ vs. time traces. The noise time traces are then analyzed using FFT. Figure 4 shows a set of such spectra obtained from data at a number of temperatures with $`H_{}`$ $`=`$ 52 kG, including the three traces shown in Fig. 2. There are two distinguishing features in the spectra. The first is that, the $`T`$ $`=`$ 30 mK spectrum, which has the broadest frequency range, clearly has a $`1/f`$-type form  over a wide range of frequency, indicating that fluctuations occur at all frequencies. The second is that all the spectra in Fig. 4 level off at low frequencies, indicating a low-frequency cut-off in the spectra which clearly shifts to lower frequencies as $`T`$ is lowered. The increase in the magnitude of the low-frequency fluctuations and the decrease in the cut-off frequency with decreasing $`T`$ in Fig. 4 indicate that the number of low-energy states increases. We suggest that these two features arise from the collective hopping of many correlated electrons. The cut-off frequency should be a measure of the characteristic frequency of fluctuations occurring on the length scale of the correlation length. A decreasing cut-off frequency with decreasing $`T`$ indicates an increasing correlation length. We also point out that the power spectra for all the temperatures shown in Fig. 4 nearly follow a single curve for $`f>`$ 10-2 Hz, with a frequency dependence of $`1/f^{1.72\pm 0.08}`$. Such behavior deviates significantly from the $`1/f`$ behavior seen in the 30 mK curve at low frequencies. Although the origin of such behavior is unknown, it is possible that this nearly temperature independent high-frequency feature is not related to collective electron hopping. Rather, it is the emerging $`1/f`$ part of the spectra at low temperatures that represents the appearance of collective electron hopping. In general, we expect the following relation between the cut-off frequency, $`f_c`$, the correlation length, $`\xi `$, and $`T`$: $`f_c\xi ^zT^{\nu z}`$, where $`\nu `$ is the correlation length exponent and $`z`$ is the dynamical critical exponent. Ideally, one measures the noise using a sensitive bridge setup  in which the large and non-fluctuating resistance background is eliminated, rather than the simple four-terminal methods used in this study. In addition, the fluctuations should be measured when the system is in a steady-state, instead of in the initial relaxation process as shown in Fig. 2. Future experiments taking the above mentioned concerns into account should provide an accurate determination of $`f_c`$ and the critical exponents.

In conclusion, we have found strong evidence for a new glassy electron state in quench-condensed ultrathin granular Al films, with a growing correlation length at low temperatures. We gratefully acknowledge numerous invaluable discussions with S. Teitel, Y. Shapir, Y. Gao, and P. Adams. We thank S. Zorba and Y. Gao who performed scanning force microscopy studies of the quench-condensed Al and Be films.


format: text

text:
# Exact Expressions for Minor Hysteresis Loops in the Random Field Ising Model on a Bethe Lattice at Zero Temperature

## I Introduction

Hysteresis is a common phenomenon but uncommonly difficult to treat analytically. Theoretically, hysteresis is expected to vanish as the frequency of the driving field goes to zero, or its period goes to infinity. However, in many cases hysteresis persists over the longest experimental time scales. For example, there are indications that hysteresis would be observed in a permanent magnet even if the applied field were to be cycled so slowly as to take the entire life of an experimentalist to complete one loop. It is of practical importance to make a theory for this kind of hysteresis which persists over the longest practical time scales, and where the effect of temperature on the area and the shape of the hysteresis loop is small. A reasonable starting point for the theory may be the Glauber dynamics of the Ising model at temperature $`T`$, driven by a field of frequency $`\omega `$. However, this is difficult to treat analytically. A simpler version of this model which appears to be adequate for our purpose is obtained by taking the limit $`T=0`$ first, and then the limit $`\omega =0`$. The zero-temperature, zero-frequency version produces realistic looking hysteresis loops if one incorporates a Gaussian distribution of on-site quenched random field in the model. Thus the zero-temperature single-spin-flip dynamics of the ferromagnetic random field Ising model (RFIM) was proposed as a model of hysteresis and Barkhausen noise by Sethna et. al. sethna . Anti-ferromagnetic RFIM is also interesting in the context of relaxation without Barkhausen noise shukla1 . The difficulty in the analytical treatment of the models comes from the presence of the quenched random field. The zero-temperature dynamics of RFIM can not be solved exactly in two or three dimensions (so far). We have obtained the major hysteresis loop for the ferromagnetic RFIM in one dimension shukla2  and on a Bethe lattice shukla3 . The Barkhausen noise on the major loop has been analyzed shukla4 . Minor hysteresis loops in the ferromagnetic RFIM have been obtained in one dimension shukla5 . Anti-ferromagnetic RFIM is apparently more difficult, and its analytic solution is limited so far to the major hysteresis loop in one dimension for a rectangular distribution of the quenched field of width $`2\mathrm{\Delta }`$ where $`\mathrm{\Delta }|J|`$, $`J`$ being the strength of the nearest neighbor interaction.

In the following, we present the solution of minor loops in the ferromagnetic RFIM on a Bethe lattice. This is an extension of the work presented in reference shukla5 , and the completion of a problem which remained open in reference shukla3 ; therefore the reader is assumed to be familiar with these references. Several important aspects of the analysis and simulations of the model which are discussed in references shukla3  and shukla5  are not repeated here in order to save space, or mentioned only briefly for the sake of completeness. It was shown in one dimension shukla5  that a reversal of the applied field by an amount $`2J`$ from anywhere on the major loop brings the system on the opposite half of the major loop. The reversed trajectory merges with the opposite half of the major loop for the portion of the reversed field exceeding $`2J`$. Thus, the width (along the applied field axis) of a minor loop which touches both halves of the major loop (but does not merge with either of them) is constant and equal to $`2J`$ irrespective of the position of the minor loop. The shape of the minor loop does depend upon its position inside the major loop. This physically interesting result is now shown to hold on all Bethe lattices irrespective of the coordination number $`z`$ of the lattice.

## II Major hysteresis loop

In this section, we recall the model briefly, and the solution of the major hysteresis loop obtained earlier shukla3 . The RFIM in an external field $`h_{ext}`$ is characterized by the Hamiltonian:

$$H=J\underset{i,j}{}s_is_j\underset{i}{}h_is_ih_{ext}\underset{i}{}s_i$$
(1)
The sum in the first term is restricted to pairs of nearest neighbors on a Bethe lattice of coordination number $`z`$. The external field $`h_{ext}`$ is cycled from $`\mathrm{}`$ to $`+\mathrm{}`$ and back to $`\mathrm{}`$. This takes the system around its major hysteresis loop. Spins turn up on the lower half of the loop, and turn down again on the upper half. The applied field changes very slowly. Equivalently, at each value of the external field, the system is allowed adequate time to attain a relaxed state with each spin pointing along the net field at its site. In the relaxed state at $`h_{ext}=h`$, the probability that an arbitrary site $`i`$ is up is given by,

$$p(h)=\underset{n=0}{\overset{z}{}}\left(\begin{array}{c}z\\ n\end{array}\right)[P^{}(h)]^n[1P^{}(h)]^{zn}p_n(h)$$
(2)
Here $`P^{}(h)`$ is the conditional probability that a nearest neighbor of site $`i`$ is up before site $`i`$ is relaxed , and $`p_n(h)`$ is the probability that the site $`i`$ with quenched field $`h_i`$ can turn up in applied field $`h`$ if $`n`$ of its nearest neighbors are up. The starting state on the lower half of the major loop has all sites down, and the starting state on the upper half has all sites up. Thus, on the lower half of the major loop, $`P^{}(h)`$ denotes the conditional probability that a nearest neighbor of site $`i`$ turns up before site $`i`$. On the upper half, $`P^{}(h)`$ denotes the conditional probability that a nearest neighbor of site $`i`$ turns down after site $`i`$. We distinguish the two situations by putting a subscript $`l`$ for the lower half and $`u`$ for the upper half. This gives

$$P_l^{}(h)=\underset{n=0}{\overset{z1}{}}\left(\begin{array}{c}z1\\ n\end{array}\right)[P_l^{}(h)]^n[1P_l^{}(h)]^{z1n}p_n(h),$$
(3)
and

$$P_u^{}(h)=\underset{n=0}{\overset{z1}{}}\left(\begin{array}{c}z1\\ n\end{array}\right)[P_u^{}(h)]^n[1P_u^{}(h)]^{z1n}p_{n+1}(h)$$
(4)
We note that $`p_{n+1}(h)=p_n(h+2J)`$, and therefore $`P_u^{}(h)=P_l^{}(h+2J)`$. Here, $`P_u^{}(h)`$ is the conditional probability that given a site which is up on the upper half of the major loop at field $`h`$, its nearest neighbor is also up. Similarly $`P_l^{}(h+2J)`$ is the conditional probability that given a site which is down on the lower half of the major loop at field $`h+2J`$, its nearest neighbor is up.

## III Starting a Minor Loop

Suppose we are on the lower part of the major loop when the applied field is reversed from $`h`$ to $`h^{}`$ ($`h^{}h`$) to generate the upper half of a minor hysteresis loop. We ask the question, what is the probability that an arbitrary site $`i`$ which was up at $`h`$ turns down at $`h^{}`$. In order to compute this probability correctly, we must take into account the irreversibility of the zero-temperature dynamics. Consider a site $`i`$ which is down on the lower half of the major loop at an applied field $`h\delta h`$ but turns up at $`h`$, where $`\delta h`$ is an arbitrarily small field. Once site $`i`$ has turned up, it may not turn down immediately if the field is rolled back to the value $`h\delta h`$. The reason is as follows. When site $`i`$ turns up, it increases the net field on each of its nearest neighbors by an amount $`2J`$. The increased field may cause one or more nearest neighbors of site $`i`$ to turn up. Suppose $`n_a`$ nearest neighbors of site $`i`$ were already up before site $`i`$ turned up, and $`n_b`$ nearest neighbors turn up after site $`i`$ turns up. Clearly, $`n_b`$ must lie in the range $`0n_bzn_a`$. The $`n_b`$ neighbors increase the local field at site $`i`$ by a finite amount $`2n_bJ`$. Therefore an infinitesimal decease in the applied field will not cause site $`i`$ to turn down unless $`n_b=0`$. A site $`i`$ with $`n_b>0`$ will turn down in decreasing applied field only after all of the $`n_b`$ nearest neighbors have turned down. When the field is reversed to $`h^{}<h`$, none of the $`n_a`$ neighbors (which turned up before site $`i`$ turned up) could possibly turn down before site $`i`$ turns down. This leaves the other $`n_b`$ neighbors which turned up after site $`i`$. The $`n_b`$ neighbors turned up because the field on them increased by an amount $`2J`$ after site $`i`$ turned up. In decreasing field $`h^{}`$, the $`n_b`$ neighbors will turn down before site $`i`$ turns down. At $`h^{}=h2J`$, all of the $`n_b`$ neighbors would have turned down with the result that the nearest neighbors of a site $`i`$ which are up at $`h^{}=h2J`$ are precisely those which were up before site $`i`$ flipped up. These neighbors will turn down on the reverse trajectory only after site $`i`$ turns down. Thus, given an up site $`i`$ at $`h^{}=h2J`$ on the upper half of the minor loop, the conditional probability $`P_u^{}(h2J)`$ that a nearest neighbor of $`i`$ is up is equal to $`P_l^{}(h)`$, where $`P_l^{}(h)`$ is the conditional probability that a nearest neighbor of site $`i`$ is up at field h given that site $`i`$ is down on the lower half of the major loop at field $`h`$. The probability that the site $`i`$ is up at $`h2J`$ is given by the equation,

$$p(h2J)=\underset{n=0}{\overset{z}{}}\left(\begin{array}{c}z\\ n\end{array}\right)[P_l^{}(h)]^n[1P_l^{}(h)]^{zn}p_n(h2J)$$
(5)
or, using the identity $`P_l^{}(h)=P_u^{}(h2J)`$,

$$p(h2J)=\underset{n=0}{\overset{z}{}}\left(\begin{array}{c}z\\ n\end{array}\right)[P_u^{}(h2J)]^n[1P_u^{}(h2J)]^{zn}p_n(h2J)$$
(6)
It follows from the above equation that the reverse trajectory will meet the upper half of the major loop at $`h^{}=h2J`$ and merge with it for $`h^{}<h2J`$. Thus, the task of computing the minor hysteresis loop is reduced to range $`h2Jh^{}h`$. We return to the question asked at the beginning of this section. What is the probability that a site $`i`$ which is up at $`h`$ turns down at $`h^{}`$ ? This is given by,

$$q^{}(h^{})=\underset{n=0}{\overset{z}{}}\left(\begin{array}{c}z\\ n\end{array}\right)[P_l^{}(h)]^n[D^{}(h^{})]^{zn}\left[p_n(h)p_n(h^{})\right]$$
(7)
Here $`D^{}(h^{})`$ is the probability that a nearest neighbor of site $`i`$ turns down on the reverse trajectory before site $`i`$. $`D^{}(h^{})`$ is determined by the equation,

$`D^{}(h^{})={\displaystyle \underset{n=0}{\overset{z1}{}}}\left(\begin{array}{c}z1\\ n\end{array}\right)[P_l^{}(h)]^n[1P_l^{}(h)]^{z1n}\left[1p_{n+1}(h)\right]`$ (10)
$`+{\displaystyle \underset{n=0}{\overset{z1}{}}}\left(\begin{array}{c}z1\\ n\end{array}\right)[P_l^{}(h)]^n[D^{}(h^{})]^{z1n}\left[p_{n+1}(h)p_{n+1}(h^{})\right]`$ (13)
Given a site $`i`$ which is up at $`h`$, the first sum above gives the conditional probability that a nearest neighbor of site $`i`$ is down at $`h^{}=h`$, i.e. at the very start of the reverse trajectory (and hence remains down for $`h^{}<h`$). The second sum takes into account the situation that the nearest neighbor in question is up at $`h`$, but turns down before site $`i`$ turns down on the return loop. Note that all the nearest neighbors of a site i which turned up after it turned up on the lower major loop must turn down before site $`i`$ turns down on the upper minor loop.

The magnetization on the reverse trajectory is given by,

$$m^{}(h^{})=2[p(h)q^{}(h^{})]1$$
(14)
## IV Completing the Minor Loop

We reverse the field $`h^{}`$ to $`h^{\prime \prime }`$ ($`h^{\prime \prime }>h^{}`$) to trace the lower half of the return loop. The magnetization on the lower half of the return loop may be written as,

$$m^{\prime \prime }(h^{\prime \prime })=2[p(h)q^{}(h^{})+p^{\prime \prime }(h^{\prime \prime })]1$$
(15)
where $`p^{\prime \prime }(h^{\prime \prime })`$ is the probability that an arbitrary site $`i`$ which turned up at $`h`$ and turned down at $`h^{}`$, turns up again at $`h^{\prime \prime }`$.

$$p^{\prime \prime }(h^{\prime \prime })=\underset{n=0}{\overset{z}{}}\left(\begin{array}{c}z\\ n\end{array}\right)[U^{}(h^{\prime \prime })]^n[D^{}(h^{})]^{zn}\left[p_n(h^{\prime \prime })p_n(h^{})\right]$$
(16)
Here $`U^{}(h^{\prime \prime })`$ is the conditional probability that a nearest neighbor of a site $`i`$ turns up before site $`i`$ turns up on the lower return loop. It is determined by the equation,

$`U^{}(h^{\prime \prime })=P^{}(h){\displaystyle \underset{n=0}{\overset{z1}{}}}\left(\begin{array}{c}z1\\ n\end{array}\right)[P_l^{}(h)]^n[D^{}(h^{})]^{z1n}[p_n(h)p_n(h^{})]`$ (19)
$`+{\displaystyle \underset{n=0}{\overset{z1}{}}}\left(\begin{array}{c}z1\\ n\end{array}\right)[U^{}(h^{\prime \prime })]^n[D^{}(h^{})]^{z1n}\left[p_n(h^{\prime \prime })p_n(h^{})\right]`$ (22)
The rationale behind equation (12) is similar to the one behind equation (8). Given that a site $`i`$ is down at $`h^{}`$, the first two terms account for the probability that a nearest neighbor of site $`i`$ is up at $`h^{\prime \prime }h^{}`$. Note that the neighbor in question must have been up at $`h`$ in order to be up at $`h^{}`$, and if it it is already up at $`h^{}`$ then it will remain up on the entire lower half of the return loop, i.e. at $`h^{\prime \prime }h^{}`$. The third term gives the probability that the neighboring site was down at $`h^{}`$, but turned up on the lower return loop before site $`i`$ turned up. It can be verified that the lower return loop meets the lower major loop at $`h^{\prime \prime }=h`$ and merges with it for $`h^{\prime \prime }>h`$, as may be expected on account of the return point memory. The exact expressions given above have been checked against numerical simulations of the model in selected cases resulting in excellent agreement in all cases which were tested.

## V Concluding Remarks.

The method of calculating the minor loop described above may be extended to obtain a series of nested minor loops. The key point is that whenever the applied field is reversed, a site $`i`$ may flip only after all neighbors of site $`i`$ which flipped in the wake of site $`i`$ (on the immediately preceding sector) have flipped back. The neighbors of site $`i`$ which remained firm after site $`i`$ flipped previously do not yield before site $`i`$ has flipped again on the return loop. We have obtained above expression for the return loop when the applied field is reversed from $`h_{ext}=h`$ on the lower major loop to $`h_{ext}=h^{}`$ ($`h2Jh^{}h`$), and reversed again from $`h_{ext}=h^{}`$ to $`h_{ext}=h^{\prime \prime }`$ ($`h^{\prime \prime }h`$). When the applied field is reversed a third time from $`h^{\prime \prime }`$ to $`h^{\prime \prime \prime }`$ ($`h^{\prime \prime \prime }<h^{\prime \prime }`$), expressions for the magnetization on the nested return loop follow the same structure as the one on the trajectory from $`h`$ to $`h^{}`$. Qualitatively, the role of $`P^{}`$ on the first leg ($`h`$ to $`h^{}`$) is taken up by $`U^{}`$ on the third leg ($`h^{\prime \prime }`$ to $`h^{\prime \prime \prime }`$) of the nested return loop.

We conclude by comparing the results obtained above with numerical simulations on Bethe lattices of coordination number $`z=3`$ and $`z=4`$, and also contrast these results with those obtained in the one dimensional case shukla5 . Let us specifically choose a Gaussian distribution of the random field with mean value zero and variance $`\sigma ^2`$. One generally expects the solution of an Ising model with nearest neighbor interactions on a Bethe lattice ($`z3`$) to be qualitatively different from its solution in one dimension ($`z=2`$), and to be similar to the mean-field solution of an infinite-range model. However, these expectations are not born out in the case of hysteresis in RFIM. The mean-field solution sethna  does not show any hysteresis for $`\sigma \sigma _c(\mathrm{})=\sqrt{\frac{2}{\pi }}`$. In contrast to this, there is hysteresis on Bethe lattices for all values of $`\sigma `$. Moreover, the behavior on lattices with $`z=2`$ and $`z=3`$ turns out to be qualitatively similar. For lattices with $`z4`$, the magnetization on each half of the hysteresis loop has a jump discontinuity for $`\sigma \sigma _c(z)`$; the jump discontinuity is absent on lattices with $`z3`$ for any finite value of $`\sigma `$. Figure 1 shows the major as well as a minor hysteresis loop on a Bethe lattice with $`z=3`$, and $`\sigma =2`$. The minor loop starts on the lower half of the major loop at $`h=1.5`$ and meets the upper half at $`h=.5`$ as may be expected from the theoretical prediction. As the magnetization on the lower half of the major loop is a single valued function of the applied field in Figure 1, the point where the applied field is reversed determines the minor loop uniquely. The one dimensional case ($`z=2`$) is similar shukla5 . However, the situation is different for $`z4`$. For $`z=4`$, we have $`\sigma _c(4)=1.78`$ approximately. Figure 2 shows the major loop for $`z=4`$ and $`\sigma =1.70`$ with a jump discontinuity at a critical field $`h_c`$ which is slightly higher than unity for $`\sigma =1.70`$. There are two values of the magnetization at $`h_c`$ both lying on the lower half of the major loop. If we reverse the applied field from the value $`h_c`$, we must specify the state of the magnetization from where the return is made giving us two possible return loops originating at $`h_c`$. Figure 2 shows two minor loops starting at $`h=.95`$ (slightly less than $`h_c`$), and $`h=1.05`$ (slightly greater than $`h_c`$) on the lower half of the major loop. Both return trajectories touch the upper half of the major loop when the field has been reversed by an amount $`2J`$ as expected from the theoretical analysis. The overall agreement between the simulations and the theory is also quite good indicating that the model considered here is self-averaging shukla3 ; shukla4 ; shukla5 .


text:
# Variational Mean Field approach to the Double Exchange Model

## I Introduction

Doped manganites show many unusual features, the most striking being the colossal magnetoresistance (CMR) in the ferromagnetic (FM) phase . In addition, the manganites have a rich phase diagram as function of band filling, temperature and chemical composition. The broad features of these phase diagrams can be understood in terms of the double exchange model (DEM) , although Jahn-Teller deformations  and orbital degeneracy may also play a role. A remarkable property of these compounds is the existence of inhomogeneities in the spin and charge distributions in a large range of dopings, compositions and temperatures. In fact, for materials displaying the largest CMR effects, the size of the phase-separated domains is so large ($`0.5\mu \mathrm{m}`$ ), that the electrostatic stability of the material should be addressed by theorists. At band fillings where CMR effects are present, $`x0.20.5`$, these compounds can be broadly classified into those with a high Curie temperature and a metallic paramagnetic (PM) phase, and those with lower Curie temperatures and an insulating magnetic phase.

The double exchange mechanism was introduced by Zenner  through the following Kondo lattice type model

$$H_{\mathrm{K}M}=\underset{i,j,\alpha }{}t_{ij}c_{i\alpha }^{}c_{j\alpha }+J_\mathrm{H}\underset{i,\alpha ,\alpha ^{}}{}𝑺_ic_{i\alpha }^{}𝝈_{\alpha \alpha ^{}}c_{i\alpha ^{}},$$
(1)
where $`t`$ and $`J_\mathrm{H}`$ are, respectively, $`e_\mathrm{g}`$ electron’s hopping and Hund’s coupling between the $`e_\mathrm{g}`$ and the localized $`t_{2g}`$ electrons responsible for the core spin $`𝑺`$.

When $`J_\mathrm{H}`$ is larger than the width of the conduction band, the model can be reduced to the double exchange model with weak inter-atomic antiferromagnetic (AFM) interactions. Early investigations showed a rich phase diagram, with AFM, canted and FM phases, depending on doping and the strength of the AFM couplings. More recent studies have shown that the competition between the double exchange and the AFM couplings leads to phase separation into AFM and FM regions, suppressing the existence of canted phases. In addition, the double exchange mechanism alone induces a change in the order of the FM transition, which becomes of first order, and leads to phase separation, at low doping. Note, however, that a systematic study of the nature of the transition at finite temperature was not addressed until recently , despite its obvious relevance to the experiments. In fact, in Ref. it was shown that a small AFM uniform superexchange interaction between the localized $`t_{2g}`$ spins is crucial to understand some of the more relevant features of the phase diagram of the manganites. In particular a first-order phase transition is found between the PM and FM phases in the range $`x0.20.5`$ is found. This transition does not involve a significant change in electronic density, so that domain formation is not suppressed by electrostatic effects. Therefore, we find a phase-separation of a rather different type of the previously discussed, not driven by a charge instability, but by a magnetic instability. In addition to this phase transition, we recover those previously discussed.

In this work we give a detailed exposition of the new mean-field technique  and emphasis is made on the importance of the Berry phase for the existence of first order phase transitions near half filling. We have been able of achieving a more complete description of the phase diagram than in previous work, because we have taken full profit of a very particular feature of the DEM, namely, fermions are bilinearly coupled to classical degrees of freedom (the Mn spins). This allows to trace-out the fermions, thus obtaining a non local effective Hamiltonian for the spins, that can be explicitly written down in terms of the density of states of the fermionic hopping matrix. What we propose to do is to calculate exactly the effective spin-Hamiltonian, for a given (disordered in general) configuration of the spins, using the so called moments-method  complemented with an standard truncation procedure . This technique can be directly used for other models, like for instance models of classical spins and lattice vibrations coupled to fermions without direct interactions , or also in contexts different from the manganites physics like the pyrochlores or double perovskytes. Once electrons are traced-out, we study the spin thermodynamics using the variational version of the Weiss Mean-Field method .

The main difference of our approach with previous work is in that we use the exact spin-Hamiltonian, while an approximated form was used up to now. For instance, the effective crystal approximation, which amounts to consider that electrons move on a perfect crystal, with a magnetically reduced hopping, was employed in  (see section VIII for a comparison with our method). A more accurate estimate of the density of states, well-known from the physics of disordered-systems and which becomes exact on the infinite dimensions limit , is the CPA. Notice that for non self-interacting electrons the Dynamical Mean-Field Approximation(DMFA)  is also equivalent to a CPA calculation on a given mean-field, which is determined via self-consistency equations (both the mean-field and the CPA density of states are obtained self-consistently). Although those approaches are reasonable, and provide useful information, they are approximated in two different ways, which is undesiderable because two different effects are entangled. First, even for only-spins models of magnetism (like the Ising Model) where the exact spin Hamiltonian can be evaluated easily, the Mean-Field approximation neglects the spatial correlations of the statistical fluctuations of the order parameter (see e.g. ). But, in addition, for an electron-mediated magnetic interaction, the evaluation of the effective spin Hamiltonian is only accurate in the limit of infinite dimensions. With our approach, the correlations on the magnetic fluctuations are to be blamed for all the differences between our results and the real behavior of the model. On the other hand, in section VIII we show how the failure of the effective-crystal approximation on finding the first-order phase transition at half-filling is due to the inaccuracy on the calculation of the density of states. Moreover, we are able to study directly in three dimensions some rather subtle details, like the non-negligible effects of keeping the Berry phase on the DEM Hamiltonian.

The structure of the paper is as follows. In Section II we introduce the DEM and our notations. In Section III we present our Mean-Field approximation. The very non-trivial part of the work, the computation of the Density of States, is explained in Section IV. It requires numerical simulations that can be performed on large lattices with a high accuracy. They are described in Section V. The effects of the Berry phase are analyzed in Section VI. Section VII is devoted to the study of the influence of the Berry phase in the phase diagram of the DEM. The comparison of the Mean-Field approach studied in this work with the de Gennes’s  and with the DMFA is carried out in Section VIII. The conclusions are summarized in Section IX.

## II Model

We study a cubic lattice with one orbital per site. At each site there is also a classical spin. The coupling between the conduction electron and this spin is assumed to be infinite, so that the electronic state with spin antiparallel to the core spin can be neglected. Finally, we include an AFM coupling between nearest neighbor core spins. We neglect the degeneracy of the conduction band. Thus, we cannot analyze effects related to orbital ordering, which can be important in the highly doped regime, $`x>0.5`$ (see however Ref. ). We also neglect the coupling to the lattice. We focus on the role of the magnetic interactions only. As mentioned below, magnetic couplings suffice to describe a number of discontinuous transitions in the regime where CMR effects are observed. These transitions modify substantially the coupling between the conduction electrons and the magnetic excitations. Thus, they offer a simple explanation for the anomalous transport properties of these compounds. Couplings to additional modes, like optical or acoustical phonons, will enhance further the tendency towards first order phase transitions. We consider that a detailed understanding of the role of the magnetic interactions is required before adding more complexity to the model. Note that there is indeed evidence that, in some compounds, the coupling to acoustical phonons or to Jahn-Teller distortions is large.

The Hamiltonian of the DEM is

$$=\underset{ij}{}𝒯(𝑺_i,𝑺_j)c_i^{}c_j+\underset{ij}{}\stackrel{~}{J}_{\mathrm{A}F}S^2𝑺_i𝑺_j$$
(2)
where $`S=3/2`$ is the value of the spin of the core, Mn3+, and $`𝑺`$ stands for a unit vector oriented parallel to the core spin, which we assume to be classical. In the following, we will use $`J_{\mathrm{A}F}=\stackrel{~}{J}_{\mathrm{A}F}S^2`$. Calculations show that the quantum nature of the core spins does not induce significant effects. In one of the earliest studies of this model , the superexchange coupling was chosen FM between spins lying on the same $`z=constant`$ plane, and AFM between spins located on neighboring planes. This is a reasonable starting point for the study of $`\mathrm{L}a_{1x}\mathrm{C}a_x\mathrm{M}n\mathrm{O}_3`$ if $`x<0.16`$, where A-type antiferromagnetism is found. For larger doping, $`0.16<x<0.5`$, which is our main focus, the magnetism is uniform and there is no a priori reason for favoring a particular direction.

The function

$$𝒯(𝑺_i,𝑺_j)=t\left[\mathrm{cos}\frac{\theta _i}{2}\mathrm{cos}\frac{\theta _j}{2}+\mathrm{sin}\frac{\theta _i}{2}\mathrm{sin}\frac{\theta _j}{2}\mathrm{e}^{\mathrm{i}(\phi _i\phi _j)}\right]$$
(3)
stands for the overlap of two spin 1/2 spinors oriented along the directions defined by $`𝑺_i`$ and $`𝑺_j`$, whose polar and azimuthal angles are denoted by $`\theta `$ and $`\phi `$, respectively. It defines a hopping matrix $`𝒯`$, whose matrix elements are $`𝒯_{ij}=𝒯(𝑺_i,𝑺_j)`$. The hopping function can be written as

$$𝒯(𝑺_i,𝑺_j)=t\mathrm{cos}\frac{\vartheta _{ij}}{2}\mathrm{exp}(\mathrm{i}\varphi _{ij}),$$
(4)
where $`\vartheta _{ij}`$ is the relative angle between $`𝑺_i`$ and $`𝑺_j`$, and $`\varphi _{ij}`$ is the so called Berry phase. It is sometimes assumed that the Berry phase can be set to zero without essential loss. It is therefore interesting to study the model that ignores the Berry phase, the hopping matrix being

$$𝒯_{ij}^{\mathrm{m}od}=\left|𝒯_{ij}\right|=\mathrm{cos}\frac{\vartheta _{ij}}{2}=\sqrt{\frac{1+𝑺_i𝑺_j}{2}}.$$
(5)
In the following sections we will analyze both models, with Berry phase (hopping matrix $`𝒯`$) and without Berry phase (hopping matrix $`𝒯^{\mathrm{m}od}`$).

## III Mean-field approximation

Our approach to the problem follows the variational formulation of the Mean-Field approximation, described for instance in . We start by writing the Grand Canonical partition function for the DEM:

$$𝒵_{\mathrm{G}C}=[dS]\mathrm{T}r^{(Fock)}\mathrm{exp}\left[(\mu 𝒩)/T\right],$$
(6)
where $`\mu `$ is the electronic chemical potential, $`𝒩=_ic_i^{}c_i`$ is the electron number operator, $`T`$ is the temperature and we use units in which the Boltzmann constant is one. The trace, taken over the electron Fock space, defines an effective Hamiltonian for the spins,

$$\mathrm{exp}\left[^{\mathrm{e}ff}(𝑺)/T\right]=\mathrm{T}r^{(Fock)}\mathrm{exp}\left[(\mu 𝒩)/T\right],$$
(7)
that can be computed in terms of the eigenvalues, $`E_n`$, of the hopping matrix, $`𝒯`$:

$`^{\mathrm{e}ff}(𝑺)`$ $`=`$ $`J_{\mathrm{A}F}{\displaystyle \underset{ij}{}}𝑺_i𝑺_j`$ (9)
$`T{\displaystyle \underset{n}{}}\mathrm{log}\left\{1+\mathrm{exp}\left[(E_n(𝑺)\mu )/T\right]\right\}.`$
Introducing the Density of States (DOS) of $`𝒯`$:

$$g(E;𝑺)=\frac{1}{V}\underset{n=1}{\overset{V}{}}\delta \left[EE_n(𝑺)\right],$$
(10)
where $`V`$ is the volume of the lattice, the effective Hamiltonian can be written as

$`^{\mathrm{e}ff}(𝑺)`$ $`=`$ $`{\displaystyle \underset{ij}{}}J_{\mathrm{A}F}𝑺_i𝑺_j`$ (12)
$`TV{\displaystyle 𝑑Eg(E;𝑺)\mathrm{log}\left[1+e^{(E\mu )/T}\right]}.`$
The Grand Canonical partition function becomes an integral in spin configuration space:

$$𝒵_{\mathrm{G}C}=[dS]\mathrm{exp}\left[_{\mathrm{e}ff}(𝑺)/T\right].$$
(13)
Thermodynamics follows from Eq. (13) as usual. The free energy, $``$, and the electron density, $`x`$, are given by:

$$=\frac{T}{V}\mathrm{log}𝒵_{\mathrm{G}C}$$
(14)
$$x=\frac{}{\mu }=𝑑Eg(E;𝑺)\frac{1}{1+\mathrm{exp}\left[(E\mu )/T\right]}$$
(15)
where $`\mathrm{}`$ stands for spectation value over equilibrium spin configurations.

The variational Mean-Field approach consist on comparing the actual system with a set of simpler reference models, whose Hamiltonians, $`_𝒉`$, depend on external parameters, $`𝒉_i`$. For simplicity, we choose the model:

$$_𝒉=\underset{i}{}𝒉_i𝑺_i.$$
(16)
The variational method is based on the inequality

$$_𝒉+^{\mathrm{e}ff}_𝒉_𝒉,$$
(17)
where $`_𝒉`$ is the free energy of the system with Hamiltonian (16), and the expectation values $`\mathrm{}_𝒉`$ are calculated with the Hamiltonian $`_𝒉`$. The inequality (17) follows easily from the concavity of the exponential function . The best approximation to the actual free energy with the ansatz of Eq. (16) is

$$=\underset{𝒉}{\mathrm{min}}\left\{_𝒉+^{\mathrm{e}ff}_𝒉_𝒉_𝒉\right\}.$$
(18)
Since, for technical reasons that will become clear in the following, it is not possible to work with one field $`𝒉_i`$ per site, we must select some subsets that contain only a few independent parameters (see Section VII). The choice is of paramount importance since it is an ansatz that will artificially restrict the behavior of the system. We have chosen the following four families  of fields, depending on a parameter, $`𝒉`$:

$`𝒉_i`$ $`=`$ $`𝒉,`$ (19)
$`𝒉_i`$ $`=`$ $`(1)^{z_i}𝒉,`$ (20)
$`𝒉_i`$ $`=`$ $`(1)^{x_i+y_i}𝒉,`$ (21)
$`𝒉_i`$ $`=`$ $`(1)^{x_i+y_i+z_i}𝒉,`$ (22)
which correspond, respectively, to FM, A-AFM, C-AFM, and G-AFM orderings. There is an order parameter (magnetization) associated to each of these orderings. We will denote them by $`M_\mathrm{F}`$, $`M_\mathrm{A}`$, $`M_\mathrm{C}`$, and $`M_\mathrm{G}`$, respectively. As a shorthand, they will be denoted generically by $``$. The order parameter is related to the corresponding $`𝒉`$ by

$$=\frac{1}{\mathrm{tanh}h}\frac{1}{h},$$
(23)
where $`h=|𝒉|/T`$. Thus, the free energy can be written in terms of $``$ instead of $`h`$ and Eq. (18) implies that it must be minimized with respect to $``$. The free energy has three contributions: the fermion free energy (FFE), the superexchange energy and the entropy of the spins:

$$()=_{\mathrm{F}er}()+NJ_{\mathrm{AF}}^2T𝒮_𝒉(),$$
(24)
where $`N`$ is , respectively, 3,$`3`$, 1 and $`1`$ for FM, G-AFM, A-AFM and C-AFM orderings.

The entropy of the spins can be easily computed in terms of the Mean-Field:

$`𝒮_𝒉()`$ $`=`$ $`(_𝒉_𝒉_𝒉)/T`$ (25)
$`=`$ $`\mathrm{log}\left[\mathrm{sinh}h()/h()\right]h(),`$ (26)
but the FFE,

$$_{\mathrm{F}er}()=T𝑑Eg(E;𝑺)_𝒉\mathrm{log}\left[1+\mathrm{e}^{(E\mu )/T}\right],$$
(27)
must be estimated numerically.

The non trivial part of the computation is the average of the DOS $`g(E;𝑺)_𝒉`$. The key is that it can be computed by numerical simulations on large lattices with high accuracy, ought to two basic facts: 1) the Mean Field Hamiltonian (16) describes uncorrelated spins and therefore equilibrium spin configurations can be easily generated on very large lattices, and 2) the DOS is a self-averaging quantity. This last point means that the mean value of the DOS can be obtained on a large lattice by averaging it over a small set of equilibrium configurations. Once the DOS is computed, the integral of Eq. (27) can be performed numerically to get the FFE as a function of $``$.

## IV Computation of the averaged DOS

The DOS can be accurately computed for any given spin configuration with the technique that we describe in the following . From its definition, Eq. (10), the DOS is a probability distribution in the variable $`E`$, whose moments are

$$\mu _k(𝑺)=𝑑Eg(E;𝑺)E^k=\frac{1}{V}\mathrm{T}r𝒯^k$$
(28)
Now, it is easy to show that the eigenvalues of the hopping matrix verify $`6tE_n<6t`$. A probability distribution of compact support can be reconstructed from its moments using the techniques of Stieljes . In practice, we only know the first $`p`$ moments, but the method of Stieljes allows us to find a good approximation to the distribution if $`p`$ is large enough.

To compute the averaged DOS we follow four steps:

* Generate spin configurations, $`\left\{𝑺\right\}`$, according to the Mean-Field Boltzmann weight, $`\mathrm{exp}(_𝒉/T)`$. This can be achieved very efficiently with a heat bath algorithm, since all the spins are decoupled in the Mean-Field Hamiltonian. In this way, one obtains spin configurations in perfect thermal equilibrium with the Boltzmann-weight given by the Mean-Field Hamiltonian.
* For each $`\left\{𝑺\right\}`$, one would calculate the moments of the DOS, using Eq. (28), and then apply the techniques of reconstruction of Stieljes . However, this is impractical since, although the matrix $`𝒯`$ is sparse, the trace in Eq. (28) would require to repeat the process $`V`$ times, and we would end-up with an algorithm of order $`V^2`$. We use instead an stochastic estimator. First, we extract a normalized random vector $`|v`$, with components

$$v_i=\frac{\alpha _i}{_{j=0}^V\alpha _j^2}$$
(29)
where the $`\alpha _i`$ are random numbers extracted with uniform probability between $`1`$ and $`1`$. Let us now call $`|n`$ to the eigenvector of eigenvalue $`E_n`$ of the matrix $`𝒯`$. It is easy to check that (the overline stands for the average on the random numbers $`\alpha _i`$)

$`\overline{v_iv_j}`$ $`=`$ $`{\displaystyle \frac{\delta _{i,j}}{V}},`$ (30)
$`\overline{n|vv|m}`$ $`=`$ $`{\displaystyle \frac{\delta _{n,m}}{V}}.`$ (31)
Then we introduce a $`v`$-dependent density of states:

$$g(E;v;𝑺)\underset{n=1}{\overset{V}{}}|n|v|^2\delta (EE_n).$$
(32)
From Eq. (31), it follows immediately that

$$\overline{g(E;v;𝑺)}=g(E;𝑺).$$
(33)
From $`v|v=1`$, and from Eq. (32), we see that $`g(E;v;𝑺)`$ is a perfectly reasonable distribution function, whose moments are

$$𝑑EE^kg(E;v;𝑺)=v|𝒯^k|v$$
(34)
Numerically, the algorithm is of order $`k\times V`$, since, as mentioned, $`𝒯`$ is sparse and only $`O(V)`$ operations are required to multiply $`v`$ by $`𝒯`$. This method allows to compute a large number of moments on large lattices. However, notice that the actual calculation is not performed this way (round-off errors would grow enormously with the power of $`𝒯`$), but as explained in the appendix.
* Reconstruct $`g(E;v;𝑺)`$ from the moments, by the method of Stieljes. The DOS is obtained in a discrete but very large number of energies, $`E`$. The cost of refining this set of energies is negligible. Hence, the integral over $`E`$ that gives the FFE, Eq. (27), can be approximated numerically with high accuracy.
* Average $`g(E;v;𝑺)`$ over the spin configurations $`𝑺`$ and over the random vectors $`|v`$. In practice we only use a random vector per spin configuration: it is useless to obtain an enormous accuracy on the density of states for a particular spin configuration, that should be spin-averaged, anyway. Since the errors due to the fluctuations of the spins and the fluctuations of the $`|v`$ are statistically uncorrelated, both of them average out simultaneously . It is crucial that the DOS is a self-averaging quantity, what means that its fluctuations are suppressed as $`1/\sqrt{V}`$. Hence, its average over a few equilibrium configurations is enough to estimate it with high accuracy.

This program can be carried out successfully on lattices as large as $`64^3`$ and even $`96^3`$, computing a large enough number of moments, $`50`$ or $`75`$. As we shall see, this suffices to achieve an excellent accuracy in the averaged DOS and in the FFE.

The whole process is repeated for several values of the Mean Fields $`𝒉`$ within each family. The FFE computed in a discrete set of magnetizations is extended to a continuous function of $``$ in the interval $`[0,1]`$ through a polynomial fit, as we will show in next section. The magnetic phase diagram of the model will come out easily then.

## V Extracting the Fermionic Free Energy

Let us discuss in this section the numerical results for the averaged DOS and for the FFE. We carried out the program designed in the previous section for 20 values of $`𝒉`$ within each family of Mean Fields, chosen in such a way that the corresponding magnetizations, $``$, are uniformly distributed in $`[0,1]`$. The expectation value $`g(E;𝑺)_𝒉`$ is estimated by averaging over 50 equilibrium (with respect to the Mean Field distribution) spin configurations. This is enough to have the statistical errors under control on a $`64^3`$ lattice.

Figure 1 displays the averaged DOS, computed on a $`64^3`$ lattice for four values of $`𝒉`$, corresponding to FM, ($`M_\mathrm{F}=0.5`$ and $`M_\mathrm{F}=1`$), paramagnetic, and G-AFM ($`M_\mathrm{G}=0.5`$) phases. They were reconstructed with its 50 first moments . Note that the DOS is even in $`E`$, as required by the the particle-hole symmetry, and that it becomes narrower in going from the FM to the G-AFM phase, as expected. The width of the density of states in the PM phase is $`8t`$, two thirds of the width corresponding to the perfect ferromagnetic, in agreement with the results of full diagonalization .

From the averaged DOS it is straightforward to compute the FFE by performing the integral entering Eq. (27) numerically. In this way, we computed $`_{\mathrm{F}er}`$ for the chosen values of $``$. Given that the free energy can be shifted by a term independent of $``$, we use $`_{\mathrm{F}er}()_{\mathrm{F}er}(0)`$ instead of $`_{\mathrm{F}er}()`$ .

To have an analytic expression for $`_{\mathrm{F}er}`$, we fit the data with a polynomial of order sixth in $``$ , with coefficients that depend on $`T`$ and $`\mu `$:

$$_{\mathrm{F}er}^{(I)}(M_I)=A_2^{(I)}M_I^2+A_4^{(I)}M_I^4+A_6^{(I)}M_I^6.$$
(35)
The index $`I`$ denotes the type of ordering: $`I=\mathrm{F},A,C`$, or $`\mathrm{G}`$.

Figure 2 shows the FFE at $`T=0`$ and $`\mu =0`$, which correspond to half filling, $`x=1/2`$, as a function of $`M_\mathrm{F}`$ or $`M_\mathrm{G}`$. The points are the result of the numerical computation and the lines are the best fits of the form (35). The high quality of the fits is remarkable. Note that the fermions favor FM order. The coefficients $`A_j^{(I)}(T,\mu )`$ of Eq. (35), which will play a major role in the exploration of the phase diagram, are displayed in Figure 3 in the cases $`I=\mathrm{F},G`$, for $`T=0`$, as a function of $`\mu `$. We always found $`A_2^{(F)}(T,\mu )<0`$ and $`A_2^{(G)}(T,\mu )>0`$, in agreement with the FM nature of the spin interaction induced by the double exchange mechanism.

Let us estimate the errors of our numerical approach. We have three sources of errors:

* Finite size of the lattice.
* Statistics, arising from the numerical simulation.
* Truncation of the infinite sequence of moments.

Finite size errors have been estimated comparing the results on a $`64^3`$ and a $`96^3`$ lattice. They turn out to be negligible, as expected given the sizes of the lattices. To estimate the statistical errors we performed three different simulations using the 50 lowest order moments. One more simulation, this time with the 75 first moments, was done in order to study the systematic error associated to the truncation of the sequence of moments. As an example, Figure 4 displays the averaged DOS in the FM phase $`(M_\mathrm{F}=0.5)`$ extracted from two different simulations, one with 50 and the other with 75 moments. We see only tiny differences, which can hardly be appreciated on the scale of the figure. This small error propagates to the FFE, which is shown in Figure 5 for $`T=0`$ and $`\mu =0`$. There are four sets of points plotted, corresponding to the four mentioned simulations. Again, the differences cannot be appreciated. The errors in the computed FFE give rise to uncertainties in the coefficients $`A`$’s of Eq. (35), which can be appreciated in Fig. 3. The largest errors appear at $`\mu =0`$ (half-filling). We have also checked that fits to a polynomial of eight-order do not change the values of $`A_j^{(I)}`$ significantly.    FIG. 4.: Averaged DOS vs. energy for $`M_\mathrm{F}=0.5`$ on a $`64^3`$ lattice from a simulation using 50 moments (solid line) and another one with 75 moments (dashed). The curves can hardly be distinguished on this scale.

To summarize, we have checked that the numerical uncertainties inherent to our numerical approach are well under control and therefore all the conclusions are robust in this sense. The discrepancies between our analysis and the true behavior of the system, if they are important, must be attributed only to the Mean-Field ansatz.

## VI Effect of the Berry phase

Sometimes it is stated that one can ignore the Berry phase in the DEM without loosing anything important. To investigate this point, let us repeat our analysis of the DEM setting the Berry phase to zero. All we have to do is to compute the DOS of the hopping matrix $`𝒯_{ij}^{\mathrm{m}od}`$ of Eq. (5). The rest is identical to what we have discussed in the previous sections.

Figure 6 displays the averaged DOS at $`M_\mathrm{F}=0`$ (PM phase) and $`M_\mathrm{F}=0.5`$ for hopping with and without Berry phase. At first sight, the differences, although noticeable (they are much bigger than the errors, cf. Fig. 4), do not seem very important. However, it happens that the results are very sensitive to small modifications of the DOS. We shall see indeed that the presence or absence of the Berry phase is crucial for some features of the phase diagram. Thus, the analysis of errors of the previous section turns out to be extremely important to give a meaning to our results. Notice that the effect of the Berry phase is stronger in the disordered phase, as expected.

Figure 7 shows the effects of the Berry phase on the FFE at $`T=0`$ and $`\mu =0`$. These effects modify the coefficients $`A_k^{(I)}`$ entering $`_{\mathrm{F}er}`$, which can be seen in Fig. 8. These coefficients are very sensitive to the modifications of the DOS induced by the Berry phase. Of especial relevance is $`A_4^{(F)}`$, which, as we shall see in the next section, governs the possibility of having first order PM-FM transitions. In particular, in the vicinity of $`\mu =0`$ this coefficient is negative. Notice however that without the Berry phase $`A_4^{(I)}`$ is negative around $`\mu =0`$ in a smaller region than with Berry phase, and it is closer to zero. This fact induces important differences in the nature of the phase transitions of the model, as we shall see in the next section.

## VII Phase diagram

The equilibrium states are determined by the absolute minima of the free energy, Eq. (24), respect to the order parameters, $``$. The minima determine the phases and the phase boundaries. Given that we know $``$ as a function of $``$, the problem of determining the equilibrium states is reduced to numerical minimization of a function of a single variable. This is indeed the way we proceed. It is however illuminating to get some insight by a semi-analytic treatment of the problem. As we have seen, to a very good approximation, the FFE is a polynomial of sixth degree in $``$. The entropy (26) can also be expanded in powers of $``$ around $`=0`$:

$$𝒮_𝒉()=\left(\frac{3}{2}^2+\frac{9}{20}^4+\frac{99}{350}^6+\mathrm{}\right).$$
(36)
Hence, we find the Landau expansion of the free energy in powers of the order parameter:

$$()=c_2^2+c_4^4+c_6^6+\mathrm{}.$$
(37)
The coefficients of the expansion are

$`c_2`$ $`=`$ $`{\displaystyle \frac{3}{2}}T+NJ_{\mathrm{A}F}+A_2^{(I)}(T,\mu ),`$ (38)
$`c_4`$ $`=`$ $`{\displaystyle \frac{9}{20}}T+A_4^{(I)}(T,\mu ),`$ (39)
$`c_6`$ $`=`$ $`{\displaystyle \frac{99}{350}}T+A_6^{(I)}(T,\mu ),`$ (40)
where $`N`$ was defined right after Eq. (24).

The free energy has the symmetry $``$. At high $`T`$, the entropic term dominates and the minimum of $``$ is at $`=0`$. As the temperature decreases, the internal energy becomes more important and the absolute minimum of $``$ can be located at $`0`$. The phase transition will be continuous when the absolute minimum at the origin changes to a maximum, i.e.:

$$c_2=\mathrm{\hspace{0.33em}0}\frac{3}{2}T_\mathrm{c}+NJ_{\mathrm{A}F}+A_2^{(I)}(T_\mathrm{c},\mu )=\mathrm{\hspace{0.33em}0}$$
(41)
At a first order transition three minima, $`=0`$ and $`=\pm _00`$, are degenerate. Three conditions must hold:

* Minimum at $`=0`$:

$$c_2>0.$$
(42)
* Minimum at $`_0`$:

$$2c_2+4c_4_0^2+6c_6_0^4=0.$$
(43)
* Degeneracy ($`(0)=(_0)`$):

$$c_2+c_4_0^2+c_6_0^4=0.$$
(44)

The solution of these three equations is:

$`_0^2=2c_2/c_4`$ (45)
$`c_6=c_4^2/(4c_2)`$ (46)
$`c_2>0;c_4<0;c_6>0`$ (47)
Eq. (45) gives the spontaneous magnetization; Eq. (46) determines the critical temperature $`T_\mathrm{c}`$ of the first order transition as a function of $`\mu `$ and $`J_{\mathrm{A}F}`$; Eq. (47) sets necessary conditions (real $`_0`$) for a first order transition to happen. We see that to have a first order PM-FM transition we must have $`A_4^{(I)}(T_\mathrm{c},\mu )<0`$. Fig. 3 shows that in particular this is possible around half filling.

The boundaries between first and second order lines are tricritical points. They are determined by the conditions:

$`c_2=0`$ $``$ $`J_{\mathrm{A}F}^\mathrm{t}={\displaystyle \frac{3}{2N}}T_\mathrm{t}{\displaystyle \frac{A_2^{(I)}(T_\mathrm{t},\mu )}{N}}`$ (48)
$`c_4=0`$ $``$ $`T_\mathrm{t}={\displaystyle \frac{20}{9}}A_4^{(I)}(T_\mathrm{t},\mu )`$ (49)
With these ingredients, we are able to discuss the phase diagram of the DEM. It has been shown in  that the phase diagram of double exchange systems is richer than previously anticipated and differs substantially from that of more conventional itinerant ferromagnets. Moreover, it is consistent with the magnetic properties of manganites  (see section IX).
FIG. 10.: Phase diagram of the DEM without Berry phase in the plane $`(x,J_{\mathrm{A}F})`$ at $`T=0`$.  We shall not repeat here the analysis of the phase diagram of the DEM carried out in . Let us concentrate on the effects of the Berry phase. Figure 9 displays the phase diagram in the plane $`(x,T/t)`$, for several values of $`J_{\mathrm{A}F}`$. The left part corresponds to the model that includes the Berry phase and on the right the Berry phase is neglected. For $`J_{\mathrm{A}F}<0.06`$ both phase diagrams are very similar. The transition temperature is slightly higher (13%) at half filling in the model that ignores the Berry phase. At low filling, both phase diagrams are almost identical. However, for $`J_{\mathrm{A}F}>0.06`$ important differences arise. A first order PM-FM transition develops around half filling if the Berry phase is taken into account. The onset for such behavior is $`J_{\mathrm{A}F}=0.06`$. In contrast, the PM-FM transition around half filling remains continuous for any value of $`J_{\mathrm{A}F}`$ if the Berry phase is neglected. The explanation is simple: the coefficient $`A_4^{(F)}`$, although negative, has a too small absolute value to drive a first order transition. This negative coefficient would become effective if the critical temperature were much lower, what can be achieved by increasing the AF superexchange coupling. But in this case the competition between FM and AF is so strong that the transition at half filling takes place between PM and A-AFM phases, and it is second order. First order PM-FM transitions only appear if the Berry phase is properly taken into account.

At $`T=0`$, the phase diagrams are similar in both cases, and we only display that of the model without Berry phase, in Fig. 10. The discussion of  applies to this case without any modification.

Let us end this section with the analysis of the phase transitions at finite applied magnetic field, $`B`$. The first order PM-FM transition around half filling survives under an applied magnetic field. In this case, the order parameter, $`M_\mathrm{F}`$, is non-zero in both phases, but suffers a jump on a line in the plane $`(B,T)`$. The line ends at a critical point, $`(B^{},T^{})`$, which has a certain magnetization $`M_\mathrm{F}^{}`$. The critical field can be measured and is of interest . Let us compute it. The free energy in the presence of a magnetic field, $`B`$, is:

$$(M_\mathrm{F})=c_2M_\mathrm{F}^2+c_4M_\mathrm{F}^4+c_6M_\mathrm{F}^6BM_\mathrm{F}.$$
(50)
The magnetic field shifts the three degenerate minima of the zero-field PM-FM first order transition and lifts the degeneracy. By tuning (increasing) the temperature it is possible to get two degenerate minima again, and a first order transition takes place. In this way, we get a transition line in the $`(B,T)`$ plane. Increasing $`B`$, the two degenerate minima become closer. At the critical field, $`B^{}`$, both minima coalesce at some point $`M_\mathrm{F}^{}`$, and the transition disappears. When this happens, the three first derivatives of $``$ respect $`M_\mathrm{F}`$ vanish, and the fourth is positive. These conditions read:

$`^{}(M_\mathrm{F}^{})=0`$ $``$ $`B^{}=2c_2M_\mathrm{F}^{}+4c_4M_\mathrm{F}^{\mathrm{\hspace{0.17em}3}}+6c_6M_\mathrm{F}^{\mathrm{\hspace{0.17em}5}}`$ (51)
$`^{\prime \prime }(M_\mathrm{F}^{})=0`$ $``$ $`2c_2+12c_4M_\mathrm{F}^{\mathrm{\hspace{0.17em}2}}+30c_6M_\mathrm{F}^{\mathrm{\hspace{0.17em}4}}=0`$ (52)
$`^{\prime \prime \prime }(M_\mathrm{F}^{})=0`$ $``$ $`24c_4M_\mathrm{F}^{}+120c_6M_\mathrm{F}^{\mathrm{\hspace{0.17em}3}}=0`$ (53)
$`^{(iv)}(M_\mathrm{F}^{})>0`$ $``$ $`24c_4+360c_6M_\mathrm{F}^{\mathrm{\hspace{0.17em}2}}>0.`$ (54)
These equations determine $`B^{}`$, $`T_\mathrm{c}`$ and $`M_\mathrm{F}^{}`$ as a function of $`\mu `$ (or $`x`$) and $`J_{\mathrm{A}F}`$. The critical temperature $`T_\mathrm{c}`$ varies very little from its value at $`B=0`$. Fig. 11 displays $`B^{}`$, in units of $`10^4t`$, versus $`x`$, for $`J_{\mathrm{A}F}=0.08t`$. In physical units, using $`t0.166\mathrm{e}V`$, $`B^{}`$ varies from 0.6 Tesla at $`x=0.33`$ to 2.2 Tesla at $`x=0.5`$. Recent measurements in $`\mathrm{L}a_{0.6}Y_{0.07}Ca_{0.33}MnO_3`$ gave a critical field of 1.5 Tesla .

## VIII Comparison with other calculations

### A Rigid band mean field approximation.

The main conclusion of the variational Mean-Field technique applied to the DEM is the prediction of a first order PM-FM transition at half-filling and its vicinity for $`J_{\mathrm{A}F}/t[0.06,0.1]`$. This is in sharp contrast with the widely used Mean-Field approach devised by de Gennes in 1960 , which predicts a second order PM-FM at half-filling for any value of $`J_{\mathrm{A}F}`$. Let us see briefly what are the differences between these two approaches that yield different qualitative behavior.

The difficulty in the Mean-Field approach to the DEM lies in calculating the contribution of the fermions to the Mean-Field free energy. In the variational method discussed here, we compute it exactly through a numerical simulation. As we have already mentioned, the only approximation is the Mean-Field ansatz for the Boltzmann weights of the spin configurations. On the other hand, de Gennes suggested that the fermion free energy might be well approximated by the free energy of an assembly of fermions propagating on a crystal with an homogeneous hopping parameter given by the average of the spin-dependent hopping parameter over the Mean-Field spin configurations. The de Gennes’ method neglects the influence of the Berry phase. In this approach. the electronic DOS depends on the spin configuration only through the hopping parameter (in this case, without Berry phase): $`g(E;𝑺)=g(E;𝒯^{\mathrm{m}od}(𝑺_i𝑺_j))`$. In mathematical terms, de Gennes’s approximation is carried out through the following substitution:

$`g(E;\left|𝒯(𝑺_i𝑺_j)\right|)_𝒉`$ $``$ $`g(E;\left|𝒯(𝑺_i𝑺_j)\right|_𝒉)`$ (56)
$`=g_0(E;𝒯_0)`$
where $`g_0(E;𝒯_0)`$ is the DOS of free fermions with hopping

$`𝒯_0(h)`$ $`=`$ $`\left|𝒯(𝑺_i𝑺_j)\right|_𝒉`$ (57)
$`=`$ $`{\displaystyle \frac{2}{J_{1/2}^2(\mathrm{i}h)}}{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{J_{l+1/2}^2(\mathrm{i}h)}{(2l1)(2l+3)}},`$ (58)
and $`J_\nu (z)`$ is the Bessel function.

Since the hopping is homogeneous, the fermion free energy is known analytically. At $`T=0`$ and half-filling ($`\mu =0`$) it is:

$`_{\mathrm{F}er}`$ $`=`$ $`{\displaystyle _{𝒯_0W_0}^0}𝑑Eg(E;𝒯_0)E`$ (59)
$`=`$ $`𝒯_0(h){\displaystyle _0^{W_0}}𝑑Eg_0(E;1)E`$ (60)
All the dependence in the magnetization is contained in $`𝒯_0(𝒉)`$. The expansion in powers of the magnetization follows straightforwardly from Eqs. (58) and (23). It yields:

$$\frac{1}{t}𝒯_0(h(M_\mathrm{F}))=\frac{2}{3}+\frac{2}{5}M_\mathrm{F}^2\frac{6}{175}M_\mathrm{F}^4\frac{18}{875}M_\mathrm{F}^6+\mathrm{}$$
(61)
The coefficient of $`M_\mathrm{F}^4`$ in $`_{\mathrm{F}er}`$ is positive. Hence, the PM-FM phase transition at half-filling can only be continuous. We have also checked that this remains true when we keep the contribution of all powers of $`M_\mathrm{F}`$ to $`_{\mathrm{F}er}`$.

The fermions in de Gennes’s approach propagate only on perfect crystals. In the truly variational Mean-Field presented in this work, the fermions propagate on the disordered spin background generated by the Mean-Field $`𝒉`$. This appears as an important ingredient that leads the predictions closer to the phenomenology, as we have shown in Ref. .

### B Dynamical Mean Field Approximation.

This method allows for an improvement on the treatment of the electronic contribution to the self energy. In the PM phase, the density of states is proportional to that in the fully ferromagnetic phase, like in de Gennes treatment. The only difference is that the constant of proportionality is $`1/\sqrt{2}`$ and not $`2/3`$. Below $`T_\mathrm{c}`$, the density of states is calculated self consistently, through a self energy which can be written as:

$$\mathrm{\Sigma }(E;𝑺_i)=|𝒯(𝑺_i𝑺_j)|^2g(E;𝑺_j)_{𝑺_j}$$
(62)
and:

$$g(E;𝑺_i)=\frac{1}{E\mathrm{\Sigma }(E;𝑺_i)}$$
(63)
Finally, the average $`\mathrm{}_{𝑺_j}`$ is carried out defining a probability distribution, $`𝒫(𝑺_j)`$, which depends self consistently on the free energy associated with a site with magnetization $`𝑺_j`$ immersed in the lattice described by $`𝒫(𝑺_j)`$.

Our approach is similar to the dynamical mean field approximation, but differs from it in two aspects:

i) The electronic density of states is calculated in a cubic lattice, instead of using the semielliptical DOS valid in the Bethe lattice with infinite coordination.

ii) We use a variational ansatz for $`𝒫(𝑺_j)`$, instead of determining it fully self consistently.

Point i) allows us to consider effects of the lattice geometry, and the influence of the Berry’s phase, as discussed above. At zero temperature, where both approaches become exact for their respective lattices, we find phases which can only be defined in a 3D cubic lattice.

If the transition is continuous, the distribution $`𝒫(𝑺_j)`$ can be expanded on the deviation from the isotropic one, $`𝒫(𝑺_j)=\mathrm{c}onst.`$ , in the PM phase. The ansatz that we use has the correct behavior sufficiently close to $`T_\mathrm{c}`$, so that both approaches will predict the same value of $`T_\mathrm{c}`$, for a given lattice. One must be more careful in the study of discontinuous transitions. Our ansatz introduces an approximation in the ordered phase (which disappears at $`T=0`$). However, near a first order transition we do not expect divergent critical fluctuations, so that our approach should give qualitative, and probably semiquantitative correct results as compared to the DMFA, in lattices where the latter is exact.

### C Hierarchy of approximations.

We are tempted to design a hierarchy of approximations, ordered according to the coefficient $`A_4^{(F)}`$ of the $`M_\mathrm{F}^4`$ term in the Landau expansion of the fermion free energy, as follows:

* de Gennes’s approximation: $`A_4^{(F)}>0`$ and the PM-FM transition is second order.
* Exact variational computation without Berry phase: $`A_4^{(F)}<0`$ but $`|A_4^{(F)}|`$ too small to produce first order PM-FM transitions, see Eq. (39).
* Exact variational computation with Berry phase: $`A_4^{(F)}<0`$ and $`|A_4^{(F)}|`$ large enough to produce first order PM-FM transitions.

## IX Conclusions

We have presented a detailed analysis of the variational Mean Field technique. This method can be useful in any situation where non self-interacting fermions are coupled to classical continuous degrees of freedom. Within this method, the fermionic contribution to the free energy is calculated exactly, and, later on, the variational Mean Field method is applied to the classical degrees of freedom. As an example, we have chosen the Double Exchange Model, both with and without Berry phase. The phase diagram has been obtained in both situations.

We have shown that the Berry phase is crucial in order to get first order PM-FM phase transitions around half filling. Such transitions are second order if the topological effects associated to the Berry phase are neglected. Thus, the dimensionality of the lattice plays a very important role in the structure of the phase diagram.

Some earlier Mean-Field computations  approximate the fermion free energy by that of an assembly of fermions propagating on a perfect crystal with an homogeneous hopping parameter averaged over the spin configurations. They yield second order FM-PM transitions in the vicinity of half-filling. The propagation of the fermions in the disordered spin background generated by the Mean-Field is another crucial ingredient to get discontinuous PM-FM transitions at half-filling. More modern approaches, as the DMFA  cannot deal with three dimensional effects as the Berry phase either.

As shown in , the variational Mean-Field described in the present work leads to results that are consistent with the phenomenology of the magnetic properties of the manganites $`\mathrm{L}a_{1x}(Sr,Ca)_x\mathrm{M}nO_3`$, in the range $`0.3x0.5`$, in particular with the fact that for materials with a high transition temperature, the PM-FM transition is continuous while for those with low $`T_\mathrm{c}`$ is not. Moreover, the order of magnitude of our estimate of the critical field for which histeretic effects disappear agrees with the experimental findings in $`\mathrm{L}a_{0.60}Y_{0.07}Ca_{0.33}MnO_3`$ . Also the phase diagram obtained by substitution of a trivalent rare earth for another one with smaller ionic radius (i.e. compositional changes that do not modify the doping level) is in remarkable agreement with our results.

Of course, the DEM itself can also be highly improved. For instance, one should include the orbital degeneracy, which is known to play an important role, and other elements like phonons and Jahn-Teller distortions. The variational Mean-Field approach can be applied with the same techniques presented in this work whenever the bosonic fields that interact with the electrons can be treated as classical.

Acknowledgements

We are thankful for helpful conversations to L. Brey, J. Fontcuberta, G. Gómez-Santos, C. Simon, J.M. de Teresa, and especially to R. Ibarra and V. S. Amaral. V. M.-M. is a MEC fellow. We acknowledge financial support from grants PB96-0875, AEN97-1680, AEN97-1693, AEN97-1708, AEN99-0990 (MEC, Spain) and (07N/0045/98) (C. Madrid).

## A The method of moments

In this appendix we include, for completeness, some details on the method of moments . For a complete mathematical background we refer to Ref. . The method of moments allows to obtain some statistical properties of large matrices (as the density of states or the dynamical structure factors, in general quantities depending on two-legged Green functions), without actually diagonalizing the matrices. Regarding the density of states, once one recognizes that it is a probability function whose moments can be obtained by iteratively multiplying by $`𝒯`$ the initial random vector $`|v`$, it is clear that the classical Stieljes techniques  can be used. Ought to the fact that the matrix $`𝒯`$ is sparse, and using the random vector trick, Eq. (34), the moments can be calculated with order $`V`$ operations. Since the spectrum of the matrix $`𝒯`$ lies between $`6t`$ and $`6t`$ for any spin configuration, the Stieljes method is guaranteed to converge. The procedure is as follows: one first introduce the resolvent

$$R(z)=_{6t}^{6t}𝑑E^{}\frac{g(E^{})}{zE^{}},$$
(A1)
that has a cut along the spectrum of $`𝒯`$, with discontinuity

$$2\pi g(E)=m\underset{ϵ0}{lim}\left[R(E\mathrm{i}ϵ)R(E+\mathrm{i}ϵ)\right].$$
(A2)
The resolvent can be obtained from the orthogonal polynomials of the $`g(E)`$, with the monic normalization:

$$P_n(E)=E^n+C_{n1}E^{n1}+\mathrm{}$$
(A3)
$$\delta _{n,m}_{6t}^{6t}𝑑Eg(E)P_n(E)P_m(E),$$
(A4)
with $`P_0=1`$, $`P_1=0`$, and $`n,m=0,1,\mathrm{}`$ The polynomials verify the following recursion relation:

$$P_{n+1}(E)=(Ea_n)P_n(E)b_nP_{n1}(E),$$
(A5)
with the coefficients $`a_n`$ and $`b_n`$ given by:

$`a_n`$ $`=`$ $`{\displaystyle \frac{_{6t}^{6t}𝑑Eg(E)EP_n^2(E)}{_{6t}^{6t}𝑑Eg(E)P_n^2(E)}},`$ (A6)
$`b_n`$ $`=`$ $`{\displaystyle \frac{_{6t}^{6t}𝑑Eg(E)P_n^2(E)}{_{6t}^{6t}𝑑Eg(E)P_{n1}^2(E)}}.`$ (A7)
The coefficient $`b_0`$ is arbitrary and is conventionally settled to one.

The resolvent has a representation in terms of a continued fraction as follows:

$$R(z)=\frac{1}{za_0\frac{b_1}{za_1\frac{b_2}{za_2\mathrm{}}}}.$$
(A8)
If one truncates the continued fraction, the resolvent would be approximated by a rationale function, which does not have a cut and use Eq. (A2) is impossible. Fortunately, when, as in this case, the density of states does not have gaps, the coefficients $`a_n`$ and $`b_n`$ tends fastly to their asymptotic values $`a`$ and $`b`$ . Thus, one can end the continued fraction  with a truncation factor $`T(z)`$, that verifies

$$T(z)=\frac{b}{zaT(z)}.$$
(A9)
Since the previous equation is quadratic in $`T(z)`$, we find that $`T(z)`$ has a branch-cut between $`a2\sqrt{b}`$ and $`a+2\sqrt{b}`$, which are the limits of the spectrum.

One should not use the moments of the $`g(E)`$ calculated with Eq. (34) to obtain the orthogonal polynomials (and hence the $`\{a_n,b_n\}`$), since this is an extremely unstable numerical procedure. It is better to use the recurrence relation:

$$P_{n+1}(𝒯)|v=(𝒯a_n)P_n(𝒯)|vb_nP_{n1}(𝒯)|v,$$
(A10)
starting with

$$P_1(𝒯)|v=\mathrm{\hspace{0.33em}0},P_0(𝒯)|v=|v.$$
(A11)
From this, one immediately gets

$`a_n`$ $`=`$ $`{\displaystyle \frac{P_n(𝒯)v|𝒯P_n(𝒯)v}{P_n(𝒯)v|P_n(𝒯)v}},`$ (A12)
$`b_n`$ $`=`$ $`{\displaystyle \frac{P_n(𝒯)v|P_n(𝒯)v}{P_{n1}(𝒯)v|P_{n1}(𝒯)v}}.`$ (A13)
In this way, one generates the $`N^{\mathrm{t}h}`$ orthogonal polynomial of the matrix (times $`v`$) recursively, at the price of $`N`$ multiplications per $`𝒯`$. The cost of this procedure is always of order $`V`$ operations. For each random-vector, one first extracts the density of states through Eq. (A2), which is subsequently averaged over the different $`|v`$ and spin realizations.

Let us finally point out that the above recursion relation is virtually identical to the Lanzcos method (the only difference lies on the normalizations). It should therefore not be pursued for a large number of orthogonal polynomials, without reorthogonalization. For the relative modest number of coefficients calculated in this work, this has not been needed.


text:
# Order by disorder from non-magnetic impurities in a two-dimensional quantum spin liquid

\[

## Abstract

We consider doping of non-magnetic impurities in the spin-1/2, 1/5-depleted square lattice. This structure, whose undoped phase diagram offers both magnetically ordered and spin-liquid ground states, is realized physically in CaV4O9. Doping into the ordered phase results in a progressive loss of order, which becomes complete at the percolation threshold. By contrast, non-magnetic impurities introduced in the spin liquids create a phase of weak but long-ranged antiferromagnetic order coexisting with the gapped state. The latter may be viewed as a true order-by-disorder phenomenon. We study the phase diagram of the doped system by computing the static susceptibility and staggered magnetization using a stochastic series-expansion quantum Monte Carlo technique.

\]

Doping of low-dimensional quantum spin liquids leads to a variety of phenomena, including unconventional superconductivity when localized spins are replaced by mobile holes , and effective random spin systems, whose low-temperature properties are dominated by large effective spin degrees of freedom, in the presence of static impurities. Such low-temperature phases have been observed recently in quasi-one-dimensional (1d) spin systems. In the two-leg ladder compound Ca14-xSrxCu24O41, a superconducting transition has been measured at high pressure . Zn-doping of the quasi-1d ladder material $`\mathrm{SrCu}_2\mathrm{O}_3`$ , and of the dimerized chain CuGeO3 , yields low-temperature phases with weak antiferromagnetic (AF) order, arising from effective couplings of random strength and sign between unpaired spins.

The question naturally arises of whether similar phenomena can be expected in quasi-2d systems. There are several prominent examples of 2d spin liquids, including (i) structurally frustrated materials, such as $`\mathrm{SrCu}_2(\mathrm{BO}_3)_2`$ , containing AF coupled spin-1/2 triangles, (ii) interaction-frustrated systems, such as the $`J_1`$-$`J_2`$ Heisenberg model on a square lattice  in the gapped regime around $`J_2/J_1=0.5`$, and (iii) structurally discretized materials, such as the 1/5-depleted AF square lattice formed by the spin-1/2 vanadium ions of CaV4O9 (Fig. 1) . These systems have resonant-valence-bond (RVB) ground states with a finite spin gap. Their AF correlation length is finite, decreasing with increasing gap size, and the low-temperature uniform susceptibility is exponentially activated. We find that non-magnetic impurities in such systems act to induce a long-ranged AF order, which exists on the subsystem of uncombined spins, and coexists with the dominant, gapped phase of the spin liquid. This order arising from disorder is the central result of our analysis. As in the quasi-1d case , its origin lies in the unfrustrated nature of the exponentially weak interactions between unpaired spins introduced when impurities break the RVB singlet units in a bipartite lattice: these are ferromagnetic when the free spins lie on the same sublattice, and AF for opposite sublattices .

In this work we consider the effects of non-magnetic impurities in the 1/5-depleted square lattice geometry (Fig. 1), motivated by the structure of CaV4O9, and by the possibility of doping this particular material with Ti. The set of superexchange parameters required to model CaV4O9 has been found  to be dominated by the frustrating, next-neighbor AF interaction, which to a reasonable approximation allows the real system to be considered as an ensemble of weakly coupled metaplaquettes (a square of side $`\sqrt{2}`$ larger than the nearest-neighbor spacing) of 4 V4+ spins. We present first a theoretical analysis for the unfrustrated system, which contains all the qualitative features we wish to address, and return thereafter to the experimental situation.

The effective magnetic Hamiltonian,

$`H=J{\displaystyle \underset{\mathrm{}}{}}𝐒_i𝐒_j+J^{}{\displaystyle \underset{\mathrm{}\mathrm{}}{}}𝐒_i𝐒_j,`$ (1)
is a spin-1/2 AF Heisenberg model with $`J`$ is the intra-plaquette and $`J^{}`$ the inter-plaquette exchange coupling constant, in which the sums marked by $`\mathrm{}`$ and $`\mathrm{}\mathrm{}`$ run over nearest neighbors within and between the plaquettes. We will explore the full parameter range $`J^{}/J0`$, which contains a spin-liquid, plaquette RVB (PRVB) phase for $`J^{}/J0.93`$ , an AF-ordered phase in the approximate interval $`J^{}/J[0.93,1.68]`$, and a further spin liquid, the dimer RVB (DRVB) phase, at larger $`J^{}/J1.68`$ . Expansions about the plaquette and dimer limits will be compared with simulations using the stochastic series-expansion quantum Monte Carlo (SSE QMC) algorithm , on 1/5-depleted, square lattices of up to 30$`\times `$30 sites, and at temperatures down to T = 0.02$`J`$. This algorithm involves expansion of the partition function in inverse temperature, employs local and global system updates, and is significantly more efficient than conventional QMC schemes. In the presence of non-magnetic impurities, ensemble averages over typically 500 random realizations were taken to ensure numerical errors smaller than the symbol sizes in Figs. 2 and 3.

Considering first an ensemble of completely decoupled 4-spin plaquettes, a small number of randomly placed, non-magnetic dopants creates free spin-1/2 moments, corresponding to the doubly degenerate $`S^z=\pm 1/2`$ ground state of 3-spin Heisenberg clusters . Each of these “pruned” plaquettes gives a paramagnetic contribution to the susceptibility $`\chi (T)`$ at low temperature. Upon further doping, clusters containing $`n<3`$ spins are also created, with distributions $`P(n)`$ , but only unpaired spins contribute to the low-$`T`$ Curie tail. Fig. 2(a) shows $`T\chi (T)`$ for various impurity concentrations. Without impurities, this quantity extrapolates to zero as the temperature is lowered, indicative of the spin gap in the clean system. With impurities, the zero-temperature Curie constant is given by

$`C=\underset{T0}{lim}[T\chi (T)]={\displaystyle \frac{1}{4}}\left({\displaystyle \frac{P(1)+1/3P(2)+P(3)}{1\rho }}\right),`$ (2)
and depends on concentration as shown in the inset of Fig. 2(a). For all $`\rho `$, three distinct temperature regimes are observed. (i) For $`T\mathrm{\Delta }_{\mathrm{max}}`$ (not shown) the spins are essentially independent, each contributing $`\frac{\mu _B^2}{4k_B}`$ to $`C`$. $`\mathrm{\Delta }_{\mathrm{max}}=1.5J`$ is the largest energy level separation in any cluster. (ii) For $`\mathrm{\Delta }_{\mathrm{min}}<T<\mathrm{\Delta }_{\mathrm{max}}`$, some spins freeze into singlets, gradually reducing the overall magnetic response by a factor of approximately $`\rho `$ as $`T`$ decreases. $`\mathrm{\Delta }_{\mathrm{min}}=0.25J`$ is the smallest energy spacing in any cluster. (iii) For $`T<\mathrm{\Delta }_{\mathrm{min}}`$, plateaus appear in $`T\chi (T)`$ when the contributing clusters reach their ground states.

With finite inter-plaquette coupling, spin-1/2 moments on adjacent plaquettes may combine to form new RVB clusters with smaller spin gaps, involving larger numbers of spins. In the susceptibility, this effective recombination reduces $`C`$ from its value for decoupled, pruned plaquettes, further suppressing $`T\chi (T)`$ at very low $`T`$. This cluster formation may be quantified from the distribution functions for neighboring plaquettes, evaluated perturbatively in $`J^{}/J`$. These functions have a non-trivial dependence on the relative positions of the dopant vacancies in adjacent clusters. Close to the PRVB limit, the inset of Fig. 2(a) illustrates this reduction of $`C`$, which is strongest at smaller $`\rho `$ where more spins are available to form larger RVB clusters. Fig. 2(b) shows analogous results for the DRVB limit.

In Fig. 3, the $`T`$-dependence of the uniform susceptibility is shown at various doping levels in the three regimes of $`J^{}/J`$. In Figs. 3(a) and (c), results from SSE QMC simulations (symbols) are plotted along with those from perturbative expansions about the limits of completely decoupled plaquettes and dimers (solid lines). At $`\rho =0`$, $`T\chi (T)`$ clearly displays activated behavior in the two gapped regions, but a linear form in the long-range-ordered phase.

PRVB: At small $`\rho `$, a plateau develops in $`T\chi (T)`$ at low temperatures, $`T\mathrm{\Delta }_{\mathrm{min}}`$, indicating a freezing of the plaquettes into their ground state configurations. For a finite inter-plaquette coupling, the additional reduction of $`C`$ observed at ultra-low temperatures indicates that the plaquette spins become correlated on an ultra-low energy scale $`J_{\mathrm{eff}}`$. At sufficiently small $`\rho `$, this energy can be estimated  by assuming RKKY-like interactions between the impurity spins, yielding

$`J_{\mathrm{eff}}J\mathrm{exp}\left({\displaystyle \frac{(1\rho )}{\xi _{AF}(\rho ,T)}}\right),`$ (3)
where $`\xi _{AF}`$ is the short-range AF correlation length in the spin-liquid regimes. Approaching $`T=0`$, the Curie constant extrapolates to $`C=\rho /12`$, a reduction of roughly 1/3 compared to the plateau value . At larger impurity concentrations, the onset of this ultra-low-$`T`$ regime moves up, consistent with $`TJ_{\mathrm{eff}}`$ (Eq. 3). At still higher doping, the picture of isolated impurity spins, and RKKY interactions mediated via the spin-liquid RVB background, breaks down. The regimes of moment formation and correlation overlap, and the plateau feature in $`T\chi (T)`$ becomes less pronounced. On doping towards the percolation threshold $`\rho _c`$, the magnetic response is dominated by the largest percolation clusters, typically containing several pruned plaquettes. From simulations in the classical limit , we find $`\rho _c=0.26(9)`$ for the 1/5-depleted lattice. Studies of percolation thresholds in regular square lattices  suggest that the value of $`\rho _c`$ for the quantum spin system differs by less than $`1\%`$ from the classical value ($`\rho _c=0.407`$). In this regime the plaquette plateaus disappear, a trend accentuated as the long-range-ordered regime is approached ($`J^{}/J0.93`$), where $`\xi _{AF}`$ increases and (Eq. 3) the onset of the ultra-low-$`T`$ regime moves to lower energies.

DRVB: As in the PRVB regime, non-magnetic impurities introduce local moments in the gapped, spin-liquid state. A dimer “plateau” is reached at $`TJ_{\mathrm{eff}}`$, and further moment recombination occurs at ultra-low temperatures, leading to an additional reduction of $`T\chi (T)`$ as $`T0`$. The plateaus disappear on approaching the transition to the AF phase, $`J^{}/J1.68^+`$.

AF: In contrast to gapped spin liquids, isolated, non-magnetic impurities do not introduce quasi-free local moments, because all spins remain aligned with the staggered AF background . However, the low-temperature susceptibility does have a divergent contribution, albeit with $`C`$ rather smaller than in the spin liquids at given $`\rho `$, arising from free moments isolated from the AF system by dopants.

Returning to the issue of order by disorder, our simulations close to the AF regime of the pure system show that introduced moments do not simply recombine into larger RVB clusters, but form a long-range-ordered AF network at ultra-low temperatures. This doping-induced long-range order is best characterized by the staggered magnetization $`m_s`$, which is shown in Fig. 4 as a function of non-magnetic impurity doping within each of the distinct phases of the system. As also in Fig. 5, these data were obtained by finite-size extrapolation at ultra-low temperatures, and thus represent the thermodynamic limit. In the long-range-ordered AF regime, we see that random static vacancies simply reduce the average staggered magnetization, driving it to zero beyond the percolation threshold $`\rho _c=0.26(9)`$. By contrast, in the spin-liquid regimes the doping-induced moments have effective interactions mediated by the RVB background, which are unfrustrated . A long-ranged order may then be expected to emerge from the large-spin clusters at zero temperature. Indeed, the resulting AF network has an extensive staggered magnetization, which peaks around $`\rho =0.08`$. This ordered AF phase coexists  with the gapped, spin-liquid state of the majority, undoped plaquettes or dimers (Fig. 3). The mechanism leading to this moment formation may be viewed as an order-by-disorder phenomenon. The inset of Fig. 4 shows the staggered magnetization as a function of the coupling ratio at the maximal doping. The disorder-induced, ordered moments in the spin-liquid regions are significant over a wide range of $`J^{}/J`$, and become large near the phase boundaries of the pure system, where they cross continuously to the full moment of the intrinsically ordered regime.

Our results are summarized by the phase diagram in Fig. 5, where the shading illustrates the strength of the staggered magnetic order. The solid lines depicting phase boundaries may not be taken as true transitions at any finite doping: impurity-induced AF order within a spin liquid crosses smoothly to impurity-suppressed order within an AF, with no vanishing order parameter. Impurities may be considered as damaging both to the AF ordered phase, and the quantum spin liquid phase.

Experimentally, the situation in CaV4O9 may be more complex than that simulated here. Inelastic neutron scattering measurements , and subsequent refinements by comparison with simulations , show that the magnetic system is composed of metaplaquettes with $`J_2=J`$ = 14meV, which interact with their neighbors via mutually frustrating couplings $`J_1=J_1^{}`$ = 0.49$`J`$, and with next neighbors by a further weak coupling $`J_2^{}=0.25J`$. Because of the frustrated nature of the subdominant interactions $`J_1`$, it would appear reasonable that the effective coupling $`J^{}`$ between metaplaquettes is small, and thus that the physical system is well in the PRVB regime; this deduction is consistent with the robust spin gap, $`\mathrm{\Delta }=9.4`$meV, of CaV4O9. Doping with non-magnetic impurities could be effected by random replacement of some spin-1/2 V4+ ions by non-magnetic Ti4+. Our simulations of the uniform susceptibility and staggered magnetization suggest that the weak AF order induced on doping should be detectable in 5-10% Ti-doped CaV4O9, by measurements of $`\chi (T)`$, Knight shift or magnetic Bragg scattering.

In summary, we have analyzed the behavior on doping by non-magnetic impurities of the spin-liquid and ordered AF phases of the Heisenberg model on a 1/5-depleted, square lattice of spins S = 1/2. In the spin liquids, introduction of impurities creates effectively free spins which combine through exponentially weak interactions to create a long-ranged AF order at $`T=0`$, coexisting with the gapped, spin-liquid phase. In the intrinsically AF regime, impurities cause a progressive destruction of the long-range order until it vanishes at the percolation threshold. The SSE QMC algorithm which we employ proves to be a powerful and efficient means of simulating the properties of a spin system; these attributes are particularly important when many realizations of random impurity distributions are required. Finally, the unique magnetic plaquette structure of CaV4O9 provides a well-defined starting point for studying the effect of non-magnetic impurities in spin liquids, and provides valuable insight for other 2d RVB compounds, such as frustrated Heisenberg systems, where a recombination of pruned bonds into new RVB clusters is also expected to occur.

We acknowledge the financial support of the Zumberge Foundation (SW, SH), SFB 484 of the Deutsche Forschungsgemeinschaft (BN), and the Japanese Ministry of Education and Science (MS).