# Effect of Dilution on First Order Transitions: The Three Dimensional Three States Potts Model

## I Introduction

The study of the impurities effect on the critical behavior of a pure material is an important issue since frequently real systems cannot be considered as pure. Nowadays the effect of dilution (disorder coupled to the energy density) on second order phase transitions is well understood. The Harris’ criterion  states that if the specific-heat of the pure system presents a power-like divergence (i.e. $`\alpha _{\mathrm{pure}}>0`$) the disorder induces a new Universality Class. Otherwise ($`\alpha _{\mathrm{pure}}<0`$) the critical behavior of the model remains unchanged. The criterium does not decide in the marginal case $`\alpha _{\mathrm{pure}}=0`$. Moreover, it is possible to show rigorously that for all the continuous phase transitions in presence of disorder, the correlation-length critical exponent $`\nu `$ verifies $`\nu 2/d`$, $`d`$ being the dimensionality of the space .

When the pure model shows a first order phase transition the situation is more complicated. However, there are a set of important results, both numerical and analytical. For instance, Aizenman and Wehr  showed rigorously that in two dimensions when introducing disorder, its conjugated density becomes a continuous function of the thermodynamic parameters. For instance, in a site diluted model the conjugate density is the energy, while in the Random Field Ising Model (RFIM), it is the magnetization. If the phase transition becomes continuous, one may ask about its Universality Class. A widely studied model in this context has been the $`q`$-states Potts model, whose pure version in $`d=2`$ undergoes a first order phase transition for $`q5`$. In recent numerical simulations using Monte Carlo or transfer matrix methods, the $`\nu `$ exponent has been found compatible within errors with the pure Ising value ($`\nu =1`$) independently of $`q`$, but for the magnetic exponent, $`\beta `$, the numerical results for $`q>2`$ are significantly different from the pure Ising value .

Unfortunately, for $`d>2`$ the situation is not so clear. Cardy and Jacobsen  (see also Cardy ) have put forward a picture of the general behavior by means of a mapping from the diluted $`q`$-states Potts model to the RFIM. Their mapping being asymptotically exact for large $`q`$, their results are also expected to hold for phase transitions with large latent heat. According to them, when a system that undergoes a first-order phase transition gets weakly diluted, the latent heat decreases. For larger dilutions, it will eventually vanish at a so-called tricritical point. Cardy and Jacobsen relate the latent heat with the order parameter of the RFIM, whose behavior is governed by the zero-temperature fixed point. In this way they are able to predict the critical exponents for the tricritical point $`\beta =\beta _{\mathrm{RFIM}}`$ and $`\nu =\nu _{\mathrm{RFIM}}/(2\alpha _{\mathrm{RFIM}}\beta _{\mathrm{RFIM}})`$. The critical-behavior of the system for larger dilutions remains unaddressed in their work. A physical realization of this scenario is provided by some magnetic semiconductors like Zn<sub>1-x</sub>Mn<sub>x</sub>Te . The magnetic atoms of these materials behave as Heisenberg spins living in a fcc lattice, with antiferromagnetic interactions. In these highly frustrated systems a first order phase transition is found in pure samples that gets second order upon dilution.

We finish this overview describing the results obtained by Elderfield and Sherrington for the diluted Potts Model in the Mean Field approximation . They found that the phase transition is first order for $`q>2`$, for all values of spin-concentration. Another interesting and related model is the Potts glass for $`q>4`$ that, according to Mean Field theory, undergoes a first order transition with no latent heat while for $`q=3`$ and $`4`$ the transition is continuous .

In this work, we shall consider the effects of site-dilution in the three dimensional three state Potts model, whose pure version presents a weak first order transition (small latent heat or large correlation length at the critical point). Our choice has been motivated by the ubiquity of weak first order transitions in nature. In particular, the pure $`q=3`$ Potts model shows very different experimental realizations appearing in very distant fields. For instance, we can cite the de-confining phase transition in quenched Quantum Chromodynamics. It can also characterize different situations in solid state physics. For instance a cubic ferromagnet with three easy axes of magnetization when a magnetic field in the diagonal of the cubic lattice is turned on (e.g. DyAl<sub>2</sub>), structural phase transitions (e.g. SrTiO<sub>3</sub>) and some fluid mixtures of five (suitable chosen) components .

We now describe briefly the phase diagram of the three dimensional three states Potts model in the temperature–concentration plane $`(T,p)`$ (see Fig. 1). The pure model undergoes a (weak) first order phase transition, at a critical temperature $`T_\mathrm{c}(p=1)`$ separating the paramagnetic high-temperature phase from the low temperature ordered one. This first order transition can be, in principle, continued inside the $`(T,p)`$ plane, where the critical temperature $`T_\mathrm{c}(p)`$ will lower for smaller $`p`$. The latent-heat for the first-order phase transition will decrease until the tricritical point. At this point the model suffers a second order phase transition that continues (belonging to another Universality Class) until the $`T_\mathrm{c}(p_\mathrm{c})=0`$ percolation limit. We remark that, in the most economic picture, this phase diagram presents three different Universality Classes: site percolation in three dimensions (which has been studied in the literature, e.g. in Ref. ), the Universality Class of the tricritical point (conjectured in Ref. ) and the Universality Class that controls the critical behavior in the region between the tricritical point and the percolation point.

In this paper we will restrict ourselves to the study of the second-order line. Although an experimental realization of the site diluted Potts model is not yet known (disorder tends to couple with the order parameter rather than with the energy), whenever it will appear it will be interesting to have clear theoretical predictions at hand. The techniques used in this paper are well suited for second order transitions, but they should be modified in the concentration range for which the phase transition is first order. Work is in progress to study this region .

## II The Model and Observables

We have studied the three dimensional Site Diluted three state Potts Model, whose Hamiltonian defined on a cubic lattice with volume $`V=L^3`$ is

$$=\mathrm{R}e\left[\underset{<i,j>}{}ϵ_iϵ_jz_iz_j^{}\right],$$
(1)   
and periodic boundary conditions are applied. In Eq. (1) $`z_i`$’s are complex roots of $`z^3=1`$, and $`ϵ_i`$’s are uncorrelated quenched random variables which are $`1`$ with probability $`p`$, and $`0`$ with probability $`1p`$. The Boltzmann weight is proportional to $`\mathrm{exp}(\beta )`$.

We have used clusters algorithms in order to update the system. In a diluted system, the set of occupied sites can present regions that are lightly connected to the percolating cluster. These regions are very difficult to equilibrate just with a single cluster algorithm . We have found that a single cluster algorithm combined with a Heat Bath sweep per measure is efficient for large concentrations. However, for small concentrations ($`p0.5`$) the previous method is not efficient enough due to the presence of intermediate-sized clusters, and we have used the Swendsen-Wang algorithm .

We have simulated at $`p=1.0,0.9,0.8,0.7,0.6,0.5`$ and $`0.4`$ in lattices $`L=8,16,32`$ and $`64`$. For $`p=0.8,0.7,0.4`$ we have also run in $`L=128`$ lattices. We have performed $`N_I=200`$ nearly independent measures in every single disorder realization. For $`p0.8`$ the number of these realizations has been $`N_S=10000`$, except for $`p=0.8`$, $`L=128`$, where we have fixed $`N_S=1000`$. In the $`p=0.9`$ case we have measured in 2000 different disorder realizations. The total amount of CPU time has been the equivalent of 16 years of 200 Mhz Pentium-Pro processor. For small dilutions we have performed the usual $`\beta `$ extrapolation  while for $`p0.5`$ we used a $`p`$ extrapolation method . Let us recall that when planning a disordered model simulation, one should balance two competing effects. First, to minimize statistical errors, it is better to work in a $`N_IN_S`$ regime. On the other hand, if $`N_I`$ is too small, the usual calculation of $`\beta `$ derivatives and extrapolations is biased. We follow the same procedure of Ref.  to eliminate the bias. With our simulation strategy ($`N_IN_S`$), it is crucial to check that the system is sufficiently thermalized while taking measures. In order to ensure this, we have systematically compared the results coming from hot and cold starts: half of our statistics for the largest lattices have been obtained with hot starts while the other half comes from cold starts.

Regarding the observables, in addition to the energy we have measured the complex magnetization

$$M=\underset{i}{}ϵ_iz_i,$$
(2)   
from which we obtain the real susceptibility as

$$\chi =\frac{1}{V}\overline{|M|^2}.$$
(3)   
We have denoted with $`(\mathrm{})`$ the thermodynamical average with fixed disorder and with $`\overline{(\mathrm{})}`$ the average over the disorder.

The formulæ for the cumulants read

$`g_2`$ $`=`$ $`{\displaystyle \frac{\overline{|M|^2^2}\overline{|M|^2}^2}{\overline{|M|^2}^2}},`$ (4)   
$`g_3`$ $`=`$ $`{\displaystyle \frac{\overline{M^3}}{\overline{|M|^2^{3/2}}}},`$ (5)   
$`g_4`$ $`=`$ $`2{\displaystyle \frac{\overline{|M|^4}}{\overline{|M|^2^2}}},`$ (6)   
$`g_4`$ being the standard Binder cumulant, $`g_2`$ measures whether the susceptibility is or not a self-averaging quantity and $`g_3`$ has been introduced since the three states Potts model is invariant under a global transformation of the $`Z_3`$ group. The other cumulants, $`g_2`$ and $`g_4`$, are also trivially invariant since we have used the modulus of the complex magnetization in their construction.

We have used a quotient method , in order to compute the critical exponents. We recall briefly the basis of this method. Let $`O`$ be a quantity diverging as $`t^{x_O}`$ ($`t`$ being the reduced temperature) in the thermodynamical limit. We can write the dependence of $`O`$ on $`L`$ and $`t`$ in the following way

$$O(L,t)=L^{x_O/\nu }\left[F_O\left(\frac{L}{\xi (\mathrm{},t)}\right)+𝒪(L^\omega ,\xi ^\omega )\right],$$
(7)   
where $`F_O`$ is a (smooth) scaling function and $`(\omega )`$ is the biggest non positive eigenvalue of the Renormalization Group transformation (the corrections-to-scaling exponent). This expression contains the not directly measurable term $`\xi (\mathrm{},t)`$, but if we have a good definition of the correlation length in a finite box $`\xi (L,t)`$, Eq. (7) can be written

$$O(L,t)=L^{x_O/\nu }\left[G_O\left(\frac{\xi (L,t)}{L}\right)+𝒪(L^\omega )\right],$$
(8)   
where $`G_O`$ is a smooth function related with $`F_O`$ and $`F_\xi `$ and we have neglected the term $`\xi _{\mathrm{}}^\omega `$ because we are working deep in the scaling region. The definition of the correlation length on a finite box that we use is :

$$\xi =\frac{\sqrt{\chi /F1}}{2\mathrm{sin}(\pi /L)},$$
(9)   
where $`\chi `$ was defined in Eq. (3) and $`F`$ is given by

$$F=\frac{V}{3}\underset{𝒌=\frac{2\pi }{L}}{}\overline{|\widehat{z}(𝒌)|^2},$$
(10)   
$`\widehat{z}(k_1,k_2,k_3)`$ being the discrete Fourier transform of $`ϵ_iz_i`$. We remark that the definition in Eq. (9) makes sense as a correlation length only in the pure paramagnetic phase of the model.

The main formula of the quotient method is

$$Q_O|_{Q_\xi =s}=\frac{O(sL,t)}{O(L,t)}=s^{x_O/\nu }+𝒪(L^\omega ),$$
(11)   
e.g. we compute the quotient between $`O(sL,t)`$ and $`O(L,t)`$ at the reduced temperature, $`t`$, in which $`\xi (sL,t)/\xi (L,t)=s`$. As particular cases of interest we cite the susceptibility, $`\chi `$, and the $`\beta `$-derivative of the correlation length, $`_\beta \xi `$, whose associated exponents are:

$`x_{_\beta \xi }`$ $`=`$ $`1+\nu ,`$ (12)   
$`x_\chi `$ $`=`$ $`(2\eta )\nu .`$ (13)   
A clean measure of scale invariance is provided by $`(\xi /L)|_{Q_\xi =s}`$. Let us recall that $`\xi /L`$ is a monotonically growing function of the inverse temperature. In the ordered phase it grows as $`L^{d/2}`$, while in the disordered phase decreases with growing lattice size. Therefore, for any pair of lattice sizes, there is a crossing temperature where $`Q_\xi =2`$. In a second order transition, $`\xi /L`$ at the crossing point should tend to a non-vanishing universal value. For a first order transition, the crossing temperatures tend to the transition point but $`\xi /L`$ at the crossing diverges due to the coexistence of ordered and disordered phases.

We finally analyze the quotient of the cumulants $`g_2`$, $`g_3`$ or $`g_4`$ at two different lattices, $`L`$ and $`sL`$, computed at the temperature where $`Q_\xi =s`$. Notice that for a second order phase transition the asymptotic limit ($`L\mathrm{}`$) of these quotients is $`1`$ corrected by terms like $`L^\omega `$ (see Eq. (11)).

The quotient method, Eq. (11), has several interesting features. First, we profit of the large statistical correlation between $`Q_O`$ and $`Q_\xi `$. Next, one does not need a previous estimate of the infinite volume critical point. Finally, it allows a simple control of the scaling corrections. All of this makes the method specially efficient for the measures of anomalous dimensions.

## III Numerical Results

The phase diagram that we have obtained numerically is shown in Fig. 1 and it agrees with the standard picture. The dashed line corresponds to the part of the transition line where we do not find a clear second order asymptotic behavior. A clear first order signature is very difficult to see even for dilutions as small as $`p=0.95`$. In this region the individual samples usually present a double peak structure. Consequently the thermalization process is very slow and a different algorithm for updating must be used. Work is in progress to analyze this region .

Our scope is now to compute the critical exponents in the region in which the transition is clearly second order, i.e. the study of the Universality Class between the tricritical and percolation limits. The first stage is to determine where an asymptotic second order behavior has been reached with lattice sizes up to $`L=128`$. In Fig. 2 we show the value of $`\xi /L`$, at the points for which $`Q_\xi =2`$ for the different $`(L,2L)`$ lattice pairs and as functions of $`L^\omega `$.

We have used for $`\omega `$ the corresponding value of the site diluted Ising model . For $`p0.7`$ we find that $`\xi /L`$ seems to tend to a dilution-independent value. Notice the clear divergence for $`p=1`$, where the transition is known to be first order. For $`p=0.9`$ we find a similar trend that for the pure case while for $`p=0.8`$ we find a transient behavior: for small lattices $`\xi /L`$ grows while in the largest lattices it seems to approach the universal value.

Another interesting quantity in order to clarify the second order behavior is the cumulant $`g_3`$. The $`g_3(L)`$ values at the points where $`\xi (2L,t)=2\xi (L,t)`$ are shown in Fig. 3. In this case we see a different scaling behavior for $`p=0.9`$ and $`p=0.8`$ up to the studied lattice sizes. We also guess from this figure that the $`\omega `$ value cannot be much larger than 0.4.

We have next considered the quotients of the different cumulants $`g_i`$ at the points where $`Q_\xi =2`$. We recall that these quantities should go to 1 as $`L`$ tends to infinity in a second order phase transition. We present our results in Fig. 4. At concentrations $`p=0.9`$ and $`p=0.8`$ we do not find an asymptotic behavior. For $`p=0.7`$, the behavior is yet not monotonous. Only for $`p=0.4,0.5`$ and $`0.6`$ it seems that the asymptotic behavior is reached. Unfortunately, a reliable estimate of $`\omega `$ cannot be obtained but our results point to a value near $`0.4`$. Moreover, the higher order scaling corrections are rather strong for these quantities. Finally, let us remark that the corrections to scaling and statistical errors are much larger for $`g_2`$ and $`g_4`$ than for $`g_3`$.

Therefore, for the study the second order region, we conclude that only for $`p0.6`$ an asymptotic scaling behavior for the considered lattice sizes has been found.

We report the results for the critical exponents as functions of $`p`$ and $`L`$ in tables I and II. We have applied Eq. (11) with $`s=2`$ to $`_\beta \xi `$ for computing $`\nu `$ and to $`\chi `$ for extracting $`\eta `$. We can observe that the asymptotic behavior of these estimates for $`p0.7`$ is not clear.

We have plotted the $`\nu `$ and $`\eta `$ apparent critical exponents for $`p=0.6,0.5`$ and $`0.4`$ as functions of $`1/L`$ in Figs. 5 and 6 respectively. We recall that we have found a $`\omega 0.4`$ value for the cumulants. From the Figs. 5 and 6 we see that $`\omega =1`$ for the leading scaling-corrections term could be a reasonable choice in this case. A possible explanation of this contradiction could be that for the observables used for computing the critical exponents the leading term ($`\omega 0.4`$) vanishes. In any case we should remark that we have not a precise control over the scaling corrections unlike, for example, in the investigation of the three dimensional site diluted Ising model . Fortunately, the scaling corrections for the critical exponents are rather small. Thus, it is not essential in this model to perform an infinite-volume extrapolation of our estimates. This is in marked contrast with the Ising case, where the extrapolation procedure was crucial to correctly compute the critical exponents.

From Figs. 5 and 6 we estimate the critical exponents in the second order region as the displayed in the last row of the table III. We remark that the $`\nu `$ value is indistinguishable from the 3D site diluted Ising model one. However, the $`\eta `$ value is very different from the values found for the rest of the models reported in the table III.

## IV Conclusions

We have numerically studied the three dimensional site diluted three states Potts model. The phase diagram in the temperature-concentration plane consists of a ferro-magnetically ordered phase separated from a paramagnetic, high temperature one. Between both regions there is a critical line, which is (weakly) first order in the limit of pure samples. For small concentrations, a clear second-order behavior is found, while the region with $`p0.9`$ shows a different behavior, probably corresponding to a cross-over, more difficult to analyze.

The critical exponents of the second order region have been computed using Finite Size Scaling Techniques. We have found that these exponents are dilution independent, and that they show a very mild evolution with the lattice size. That is why a sound estimate of the critical exponents can be given, in spite of the fact that we have been unable of measuring the scaling-corrections exponent $`\omega `$. This is in marked contrast with the situation in the site-diluted Ising model, where the scaling-corrections are severe but $`\omega `$ can be obtained with a $`15\%`$ accuracy.

Regarding the variation of the critical exponents with $`q`$ we have compared the results for the Potts model and the Ising case. We have found that the $`\nu `$ exponent varies slowly (or perhaps remains unchanged) with the $`q`$ value whereas the $`\beta `$ or $`\eta `$ exponents show a clear variation on this parameter. This picture strongly reminds the obtained in numerical simulations in two dimensions.

The study of the first-order region and of the critical behavior in the neighborhood of the tricritical point requires rather different numerical techniques and will be the matter of future work.

## V Acknowledgments

We gratefully acknowledge discussions with D. Belanger, J. Cardy and H. Rieger. We thank partial financial support from CICyT (AEN97-1708 and AEN97-1693). The computations have been carried out using the RTNN machines (Universidad de Zaragoza and Universidad Complutense de Madrid) and the ORIGIN2000 at the Centro de Supercomputación Complutense (CSC).