format: text

text:
# Absence of static phase separation in the high Tc cuprate YBa2Cu3O6+y

\[

## Abstract

We use 89Y NMR in YBa2Cu3O6+y in order to evaluate with high sensitivity the distribution of hole content $`p`$ in the CuO2 planes. For $`y=1`$ and $`y=0.6`$, this hole doping distribution is found narrow with a full width at half maximum smaller than $`\mathrm{\Delta }p=0.025`$. This rules out any large static phase separation between underdoped and optimally doped regions in contrast with the one observed by STM in Bi2212 and by NQR in LaSrCuO. This establishes that static electronic phase separation is not a generic feature of the cuprates.

\]

In the cuprates, the hole doping of the CuO2 planes induces the unusual features of the phase diagram: high-$`T_c`$superconductivity, strange metal behavior, or pseudogap. A large body of theoretical work argues that this hole doping could be intrinsically strongly inhomogeneous in the planes, forming segregated hole-rich and hole-poor regions on a nanoscale . This phase separation has been proposed to be essential to explain the unusual properties of the cuprates, appearing for example as stripes as argued from neutron scattering experiments . Such proposals have been recently highlighted by STM studies which reveal strong inhomogeneities of the superconducting properties at the surface of Bi2212 . The Berkeley STM group have imaged a spatial distribution of the superconducting gap which they associated with a distribution of concentration of holes ranging from $`p=0.1`$ to $`p=0.2`$ holes/planar unit cell .

The archetype of the cuprate families is YBa2Cu3O6+y (YBaCuO). If these electronic inhomogeneities do not exist in this YBaCuO family, then they result of specific disorder in other compounds. However, STM measurements are not as extensive in YBaCuO due to surface cleaving problems and oxygen loss in vacuum. The existing STM studies display an inhomogeneous or relatively homogeneous surface depending on the surface preparation procedure . These limitations do not occur if one uses other local probes such as nuclear magnetic (NMR) or quadrupolar resonance (NQR). The huge body of NMR/NQR studies done so far concentrated on the $`p`$variation of the average values of specific quantities (i.e. the static NMR shift $`K`$ or the relaxation time $`T_1`$). However, as will be emphasized hereafter, the NMR/NQR spectroscopy also allows to determine the distribution of these quantities in the bulk samples. NMR is then sensitive to the local distribution of electronic properties like STM, but not to the specificities of the surface. As an appealing example, Singer et al   evidenced recently a distribution of $`T_1`$ over the Cu NQR spectrum in bulk LaSrCuO, which can be attributed to a distribution of $`p`$ as large as the one observed on the Bi2212 surface.

We present an NMR static study of the YBaCuO family using 89Y NMR spectra. This allows us to estimate the shape of the hole distribution in the bulk of the compound. We study two doping compositions (underdoped O6.6 and slightly overdoped O7). At O7 the oxygen chains reservoirs are full, thus well ordered. At O6.6, the various properties are weakly oxygen content dependent. So both compositions should produce a minimal disorder of doping in the planes. We show indeed that our results reveal a very narrow hole distribution. These results will be compared with other experiments in the various cuprate families.

We synthesized a large sample batch of single crystal grains of YBaCuO. It was oxidized at $`T=350{}_{}{}^{}C`$, and then slowly cooled under oxygen atmosphere from 350$`{}_{}{}^{}C`$ to 270$`{}_{}{}^{}C`$ at a rate of $`0.4K`$/hr. After this procedure, the sample has a nominal maximum oxygen content $`y=1`$with a reduced $`T_c=89.1`$ K with respect to optimal doping ($`T_c=92.5`$K). A part of this batch was reduced to $`y=0.60\pm 0.02`$ ($`T_c=55.8`$K) under primary vacuum ($`P=0.13`$ mbar) at $`T=385{}_{}{}^{}C`$ . The oxygen content was measured by thermogravimetry, i.e. measurement of weight loss during the deoxidisation process until equilibrium is reached in the above conditions. For both samples, the crystallites were then aligned in stycast epoxy in the field $`H_{ext}(7.5`$ Tesla) of the NMR spectrometer . This field corresponds to a reference NMR resonance $`{}_{}{}^{89}\nu =15634.67`$ kHz for a liquid YCl3 solution. Spectra were obtained for Hext perpendicular to the $`c`$ cristallographic axis using a standard $`\pi /2\pi `$ pulse sequence and Fourier transform of half of the spin echo  The repetition time of the pulse sequence was long enough at each temperature to recover the saturation signal for any value of $`p`$. For a given T, the maximum of T$`{}_{1}{}^{}(p)`$ varies from $`17`$ to $`110`$ sec when decreasing $`T`$ from $`300K`$ down to $`120K`$ .

89Y NMR probes sensitively any doping distribution in the CuO2 planes because the Y nucleus is coupled to Cu of its adjacent CuO2 planes through the oxygen orbitals . This results in a shift of the NMR line related to the Cu spin susceptibility $`\chi `$ given by

$${}_{}{}^{89}K(T)=\frac{8A_{hf}\chi (T)}{\mu _B}+^{89}\delta $$
(1)
Here, $`{}_{}{}^{89}K(T)`$ is strongly dependent on hole content because of the strong hole doping dependence of $`\chi (T)`$. In contrast, the $`T`$independent chemical shift $`{}_{}{}^{89}\delta `$ and hyperfine coupling $`A_{hf}=1.95`$ $`kOe/\mu _B`$ between the $`{}_{}{}^{89}Y`$ and one Cu do not change much with hole doping ($`\mu _B`$is the Bohr magneton) . This relation is illustrated in Fig.1 which displays the shift of the peak value of the 89Y NMR line. It has been taken for a series of samples prepared from a single batch for different $`y`$. The shift $`{}_{}{}^{89}K(T)`$ displays the well known behaviour for $`\chi (T)`$. It is nearly constant above optimal doping. It displays a marked decrease at low T which is the signature of the pseudogap for the underdoped cases. The oxygen content $`y`$ can be converted into the hole content per plane using the parabolic dependence of $`T_c`$  given by $`p=0.16\pm \sqrt{(1T_c(y)/T_c^{\mathrm{max}})/82.6}`$. The corresponding $`p`$ and experimental $`T_c`$ are plotted in the inset of fig.1. Note that near optimal doping, the large $`p`$ dependence of $`{}_{}{}^{89}K`$ makes it a sensitive probe to any doping variation: for example, at $`T=120`$ K, a change in $`p`$ as small as $`0.01`$ corresponds to a change in $`{}_{}{}^{89}K`$ of about 20 ppm, easily detectable. In the underdoped regime, the sensitivity to the doping is smaller but still sizeable.

The NMR spectrum of a given sample consists of a histogram of the NMR shifts throughout the sample. Note that, as STM, the NMR local probe is sensitive to variations nearly on the atomic scale, as each Y nucleus is coupled only to its 8 nearest neighbour Cu sites. Then, from the correspondence shift-doping of fig.1, any local distribution of doping will lead to a similar distribution in the shift, hence to a broadening of the spectrum.

Let us now present the experimental 89Y NMR spectra for $`y=0.6`$ and $`1`$ from which we will extract quantitatively the actual doping distribution. The spectra are plotted for $`T=120`$ K and $`300`$ K for the optimally and underdoped samples in fig.2. As can be observed immediately there is no significant overlap between the two spectra at $`T=120`$ K . This already proves that the O7 ($`p=0.18`$) sample contains no appreciable amount of underdoped nanoscale regions (equivalent to $`p=0.1`$) and vice versa. In fig. 2 we also display the hole doping scales deduced from the relation between $`p`$ and $`{}_{}{}^{89}K`$ obtained from fig.1. From these hole doping scales it can be seen roughly that the full width at half maximum FWHM of the doping distribution is necessarily smaller than $`\mathrm{\Delta }p0.04`$ at optimal doping and $`0.02`$ for the underdoped case. This clearly proves without any further analysis that the observed doping distribution in YBaCuO is much smaller than in Bi2212 and LaSrCu3O where $`\mathrm{\Delta }p0.1`$.

For a more quantitative analysis, we need to consider besides the distribution of $`\chi (T)`$ any distribution of the chemical shift $`\delta `$ and hyperfine coupling $`A_{hf}`$which enter Eq.1. Such distributions will come from any local structural disorder. These distributions are fully responsible for the FWHM $`\mathrm{\Delta }K=32\pm 5`$ $`ppm`$ of the undoped compound for $`T>T_N`$ . Indeed, for oxygen content $`y<0.15`$, no hole distribution broadening is expected . This results from the fact that xoxygens introduced in the Cu reservoir layer mainly convert the two adjacent 3d10 Cu(1) into 3d9 and do not yield any hole doping of the CuO2 planes, so that $`p`$ is strictly $`0`$. It is confirmed by the fact that the spectrum has the narrowest width amoung all dopings and is not found significantly dependent on oxygen content up to $`y=0.15`$.

For higher $`y`$, if $`\mathrm{\Delta }\delta `$ and $`\mathrm{\Delta }A_{hf}`$ are the FWHM of the respective distributions assumed to be uncorrelated and gaussian for simplicity, this will lead to an additional broadening with FWHM $`\mathrm{\Delta }K(y,T)`$:

$$\mathrm{\Delta }K^2(y,T)=\mathrm{\Delta }\delta ^2(y)+\mathrm{\Delta }K_{hf}^2(y,T)$$
(2)
where $`\mathrm{\Delta }K_{hf}=8\mu _B^1\mathrm{\Delta }A_{hf}\chi (T)`$. We note that the $`y=0`$ compound is the more ordered as it is tetragonal with no twin boundaries and empty chains, hence $`\mathrm{\Delta }\delta `$ and $`\mathrm{\Delta }A_{hf}`$ should increase at higher $`y`$. As these broadenings add to the doping induced broadening, taking their $`y=0`$ estimates at higher dopings will lead to an overestimate of the doping distribution. As $`\mathrm{\Delta }K_{hf}`$ is T-dependent, we can evaluate $`\mathrm{\Delta }A_{hf}`$ and $`\mathrm{\Delta }\delta `$ separately. Indeed, at $`y=0.6`$, $`\chi (T)`$ doubles between 120 K and 300 K so that $`\mathrm{\Delta }K_{hf}`$ doubles as well, whereas the FWHM increases only by $`11\%`$ as seen in fig.2. As this T-dependence might as well be due to the doping distribution, this leads to an upper bound $`\mathrm{\Delta }A_{hf}/A_{hf}0.13`$. Such an upper bound would be explained in a naive hybridization computation by a random displacement of Cu by 0.04 Å from its ideal position, a typical value from Rietveld measurements. We will consider the two extreme cases:

(i) $`\mathrm{\Delta }A_{hf}/A_{hf}=0`$ corresponding to $`\mathrm{\Delta }\delta =32`$ ppm

(ii) $`\mathrm{\Delta }A_{hf}/A_{hf}=0.13`$ corresponding to $`\mathrm{\Delta }\delta =29`$ ppm.

In order to extract the actual doping distribution, we start from a distribution of oxygen doping $`P(y)`$. We convert it for different temperatures into a shift distribution through the phenomenological $`{}_{}{}^{89}K`$ vs $`y`$ variation deduced from Fig.1. We convolute this $`K`$ distribution with the $`\delta `$ and $`A_{hf}`$ gaussian distributions with FWHM either (i) or (ii). This process is iterated until an optimal $`P(y)`$ is found to fit the experimental spectra for all temperatures with no additional free parameter.

In the underdoped $`y=0.6`$ compound, this fitting procedure is limited by the fact that the NMR shift is nearly insensitive to the doping distribution for $`y<0.6`$ as seen in fig.1. However, the thermogravimetry procedure during deoxydation constrains the measure of the average oxygen level $`\overline{y}=0.60\pm 0.02`$. We thus assume $`P(y)`$ to be symetric around $`\overline{y}`$. The best fit is then obtained for a gaussian distribution $`P(y)=\mathrm{exp}\left((y0.6)^2/\sigma ^2\right)`$ with $`\sigma =0.05`$ for (i) and $`\sigma =0.1`$ for (ii). This oxygen distribution $`P(y)`$ and the corresponding hole distribution $`𝒫(p)`$ are plotted in fig.3 together with the usual $`T_c`$ diagram. The corresponding simulations fit perfectly the experimental ones at all temperatures (examples of fits are given as dotted lines in fig.2).

For YBaCuO7, the chains are completely full. Further oxidation is prohibited. Therefore the distribution $`P(y)`$ is expected to be non symmetric . To take this limitation into account, we model $`P(y)`$ as a convolution of a gaussian by $`\mathrm{exp}(\left|y1\right|/\lambda _\pm )`$ where $`\lambda _\pm `$ are allowed to differ for $`y>1`$ and $`y<1`$. The best fit is found for $`\sigma =0.01,`$ $`\lambda _{}=0.06,`$ and $`\lambda _+=0.025`$ (i) or $`0.004`$ (ii)$`P(y)`$ and the corresponding fit of the spectra are respectively plotted in fig.3 and fig.2. The high asymmetric shape towards $`y<1`$ found for $`P(y)`$ confirms that no large overdoping of the planes is produced locally, as expected. In order to illustrate the high accuracy of our method, we choose $`P(y)`$ slightly broader than our best fit, plotted as the dotted line in the upper panel of fig.3. The corresponding simulation plotted as a dashed line in the upper panel of fig.2 clearly fails to fit the experimental spectrum.

In summary, our results demonstrate that the maximum possible distribution of doping is quite sharp with typical width $`\mathrm{\Delta }p0.025`$ for optimal doping and $`\mathrm{\Delta }p0.01`$ for the underdoped sample. The distribution is much smaller than in the LaSrCuO and Bi2212 cuprates. In LaSrCuO, at optimal doping, Singer et al. find $`\mathrm{\Delta }p0.09`$ as figured by the arrow in fig.3 . Another NMR study  using 17O and 63Cu NMR allowed to evidence a short length scale spatial modulation in underdoped LaSrCuO compatible with Ref.. In Bi2212, the STM experiments reveal regions on the surface with the usual superconducting gap while others display a pseudogap . This has been interpreted to be due either to some disorder effect  or to superconducting and insulating-like underdoped regions . The corresponding doping distribution of width $`\mathrm{\Delta }p=0.085`$ extracted in Ref. is plotted in fig.3. Similar widths are found for an optimally and an underdoped compound in Ref.. Therefore, YBaCuO appears much more homogeneous than both LaSrCuO and Bi2212. In other cuprates, no key experiment was performed to probe locally the doping distribution, to our knowledge. However, in Hg1201 and Tl2201, typical 17O NMR widths are similar to those in YBaCuO and much smaller than in the Bi and La compounds . Following the analysis done above and the proportionality between 17O and 89Y shifts, we then expect the Hg and Tl family not to exhibit any strong doping distribution either.

In the LaSrCuO family, the large distribution seen by NQR occurs in the bulk of the material and might be associated with tilt of the oxygen octaedra, buckling of the planes, or stripes. In the Bi material, the STM results might be only a surface specificity. In YBaCuO, the doping distribution $`P(y)`$ measured here is not only much narrower but might even have only a macroscopic origin. As we used powders of cristallites with sizes smaller than 30 $`\mu m`$, any oxygen gradient within each cristallite, or a correlation between $`\overline{y}`$ and the cristallite size could lead to the observed $`P(y)`$. At O7, such oxygen gradients are explained by the fact that oxygen diffusion between chains becomes limited at low temperature. This naturally leads to the asymmetry of $`P(y)`$ towards $`y<1`$ seen in fig.3. This effect has been observed systematically in the various studies published by our group using 89Y, or 17O NMR spectra in many different YBCO7 powders. At O6.6, the observed broad $`P(y)`$ is very narrow when plotted versus $`p`$ (fig.3). Indeed a change in oxygen content does not strongly modify the actual hole doping of the planes near $`y=0.6`$ hence leading to the plateau of $`T_c`$. This is to be associated with the fact that extra oxygens at such composition occupy empty chains, and do not modify the hole content, in analogy with the situation at $`y=0`$. Therefore, the actual distributions would be much narrower on a submicron size region similar to the one sampled by STM. For both dopings, the key factor to the good homogeneity encountered is then probably the presence of the chains. The existing chain disorder does not lead to a sizeable distribution of hole content in the planes. Furthermore, those chains are probably sufficiently far from the planes not to stabilize a charge segregation. This specificity makes YBa2Cu3O$`_{6.6\text{ or }7}`$ some of the best prototypes of clean homogeneous cuprates.

In conclusion, the nanoscale static segregations observed so far cannot be considered as an intrinsic phenomenon common to all cuprates. Contrarily to the observations in Bi2212 and LaSrCuO, the maximum hole doping distribution $`\mathrm{\Delta }p`$ found in YBCO is small enough to conclude that our YBCO samples do not consist in interleaved regions with qualitatively different physical properties (metal versus insulator, with or without pseudogap, etc). The only remaining possibility for charge segregation or stripe-like scenarii to apply inYBaCuO is then a dynamical process where the phase separation would have a lifetime smaller than our timescale of observation. For the present experiment, it corresponds to the inverse spectral width, typically 1 m$`\mathrm{sec}`$. Such dynamics might have been detected in inelastic neutron scattering experiments, but the existing results are still under debate .

Fig. Legend :

fig.1 : 89Y NMR shift $`{}_{}{}^{89}K`$ for $`H_{ext}c`$ for different oxygen contents $`y`$ versus temperature. In the inset, $`T_c`$ for each sample is plotted versus hole doping with the same symbol as the one used for $`K`$.

fig.2 : 89Y NMR Fourier Transform spectra are plotted at different temperatures and dopings in YBa2Cu3O6+y for $`H_{ext}c`$ (black solid lines) with arbitrary normalization. In each panel, a doping scale gives the relation between the shift $`{}_{}{}^{89}K`$ and hole doping $`p`$. Dotted lines are simulated spectra (see text). At O7:$`T=120`$ K, an additional dashed line is plotted, which corresponds to the simulation made using the dotted distribution of the upper panel of fig.3.

fig.3 : The full black dots represent $`T_c`$ versus oxygen content (upper panel) and hole doping (lower panel) for YBa2Cu3O6+y. Distributions of oxygen content $`P(y)`$ and the corresponding distribution of hole content $`𝒫(p)`$ were used to fit the spectra for $`y=1`$ and $`y=0.6`$ samples (see text). The dark and light gray enveloppes are obtained with a minimal distribution of chemical shift and hyperfine coupling within assumptions (i) and (ii) respectively. The dashed distribution and the arrow shown in the lower panel represent the distributions measured in Bi2212  and LaSrCuO .

\*present adress : Dept of Physics, IIT Bombay 400076 India.


text:
# Impurity-Semiconductor Band Hybridization Effects on the Critical Temperature of Diluted Magnetic Semiconductors.

## I Introduction.

It has been recently discovered that some Mn doped semiconductors, such as Ga1-xMnxAs, undergo a ferromagnetic-paramagnetic transition at temperatures near 100K Matsukura et al. (1998); Ohno (1998). Motivated by this high Curie temperature, there have been further experimental studies in III-V magnetic semiconductors, and ferro-paramagnetic transition temperatures as high as 1000K have been obtained in (Ga,Mn)N Sonoda et al. . The optimal Mn concentration is near 5$`\%`$ in these materials, making magnetic ions to be rather diluted in the host semiconductor.

The existence of ferromagnetic semiconductors at room temperature opens the possibility of combining the spatial modulations of both carrier and spin density. Therefore, diluted magnetic semiconductors (DMS) have attracted much attention for their potential use in spintronic devices. Some spintronic devices have been already implemented Koshihara et al. (1997), and different heterostructures in order to optimize and control the critical temperature have been proposed da Silva et al. (2001); Fernández-Rossier and Sham (2001); Brey and Guinea (2000); Jungwirth et al. (1999).

In this work we study magnetic properties III-V magnetic semiconductors, in particular of GaMnAs. The high critical temperature DMS’s have a high concentration of randomly located manganese ions. From electron paramagnetic resonance and optical experiments Szczytko et al. (1999); Linnarsson et al. (1997), it is known that the Mn ions in these semiconductors have a $`S=5/2`$ local moment. The Mn ions go in substitutionally at the cation sites, Ga, and contribute holes to the semiconductor valence band. Although the local moment of the magnetic ions in GaMnAs indicates that each Mn should contribute with a hole, it is experimentally observed that the density of itinerant holes, $`p`$, is a small fraction ($`10\%`$) of the Mn concentration Ohno and Matsukura (2001). The reduction in the carrier density is probably due to the presence of antisite defects which, acting as deep donors partially compensating the Mn acceptors, reduce the hole density Ohno et al. (1996). A large number of antisite defects is expected to be present in these systems because, in order to minimize the mobility of the Mn ions and the existence of clusters, the semiconductors are fabricated using molecular beam epitaxy growing techniques at low temperatures, for which a large number of these defects is known to occur Ohno et al. (1996); Esch et al. (1997).

It is widely accepted that the ferromagnetism in III-V semiconductors is induced by the motion of the electrical carriers in the system. However, there is no agreement on the precise role of the carriers. They can be either itinerant, with a clear character of the host semiconductor valence band states Dietl et al. (2000, 2001); Konig et al. (2000); Jungwirth et al. , or rather they can be mostly located on Mn impurities, with a clear impurity band magnetic, polaron character Berciu and Bhatt (2000); Bhatt et al. ; Kaminski and Sarma ; Yang and MacDonald . On the other hand, the carrier induced ferromagnetic interaction between Mn2+ ions in magnetic semiconductors has some resemblance to the one observed in manganese perovskites Tokura (1999), double perovskites Kobayashi et al. (1998), and pyrochlores Shimikawa et al. (1996). In these systems, the phase diagram is different from that of conventional itinerant ferromagnets, showing first-order transitions, phase separation Guinea et al. (2000), and formation of magnetic polarons near the Curie temperature Majumdar and Littlewood (1998); Calderón et al. (2000).

In order to study these issues, we use a tight-binding model Hamiltonian in which, besides the Hund’s coupling, a site energy difference between magnetic and normal ions, $`W`$, is introduced. Monte Carlo (MC) simulations are then performed for this model, with particular emphasis on the $`W`$ dependence of the phase diagram. A mean field calculation within the coherent-potential-approximation (CPA) is then carried out and compared with MC results. The CPA treatment confirms the robustness of the MC results for the effect of $`W`$, in addition to allowing a simpler exploration of important issues such as the existence of phase separation.

Our main results are the following:

1) The hole charge on the Mn ions, controlled by the impurity potential $`W`$, has a profound effect on the transition temperature. Upon increasing $`W`$, the hole charge on Mn decreases, leading to a substantial decrease of the Curie temperature.

2) For much of the model’s parameter space, the MC critical temperature is higher than that obtained using the simplest Weiss mean field theory. The MC critical temperature is always higher than the corresponding CPA temperature. Nevertheless, all trends in the MC results, particularly the effect of the impurity potential, are reproduced in the CPA calculation. We believe this to be a robust feature.

3) The CPA calculation shows the existence of small pockets of phase segregation instability close to band edges.

The paper is organized in the following way. In section II we introduce and explain the model. Section III and IV describe the MC and CPA methods, respectively. In Section V we present our results, both MC and CPA, making contact with previous calculations and/or experiments. Section VI summarizes our work.

## II The model.

It is generally accepted Dietl et al. (2000, 2001); Konig et al. (2000); Jungwirth et al. ; Berciu and Bhatt (2000); Bhatt et al. ; Kaminski and Sarma  that the DMS’s are governed by the following Hamiltonian:

$$H=H_h+H_{ions}+H_{hions},$$
(1)
$`H_h`$ and $`H_{ions}`$ are the parts of the Hamiltonian describing the holes and the Mn ions, respectively. The term $`H_{hions}`$ represents the interaction between the magnetic ions and the carriers.

$`H_h`$ is the sum of the kinetic energy of the carriers, the hole-hole interaction energy, and the interaction of the holes with the disorder potential from randomness in the host semiconductor, in particular from the presence of antisite defects. In the actual magnetic semiconductors, the carrier density is of the order of $`10^{20}`$ cm$`^3`$, and, for these values, it is justified to neglect the effect of the carrier-carrier interaction. In this work we do no take into account the interaction of carriers with antisite disorder. This contribution is probably important in actual samples, but here we are mainly interested in the effect of the on-site energy difference and thermal fluctuations on the critical temperature. We describe the motion of the holes using a single one-band hole dispersion, instead of the more appropriated six bands envelope function formalism M. Albolfath and MacDonald (2001). This one-band model is probably insufficient to give a quantitative prediction for the Curie temperature, but it includes the main features of the hole system we want to study: the effect of thermal fluctuations, Mn-GaAs orbital hybridization, and charge localization effects. With all this, the Hamiltonian for the holes is written as

$$H_h=t\underset{}{}(C_{i,\sigma }^+C_{j,\sigma }+h.c.),$$
(2)
where the sum runs over all first neighbors pairs on a simple cubic lattice, $`t`$ is the tunneling amplitude, and $`C_{i,\sigma }^+`$ creates a hole with spin $`\sigma `$ at site $`i`$. This model has a hole bandwidth of $`12t`$, and the parameter $`t`$ has been chosen to give a bandwidth of the order of that of GaAs ($``$10$`eV`$), implying $`t`$0.84$`eV`$. The band edges of this model correspond to a parabolic dispersion with an effective mass $`m^{}=\mathrm{}^2/2a^2t`$, being $`a`$ the lattice parameter of the simple cubic lattice. The value of $`a`$ has been chosen to make $`a^3`$ the volume per GaAs unit, i.e. $`a`$=3.56Å. With these parameters, we obtain an effective mass $`m^{}`$=0.36$`m_0`$.

The term $`H_{ions}`$ in Eq.(1) describes the direct antiferromagnetic interactions between the magnetic moments of the Mn ions. Those are much smaller than the interactions with the carrier spins and, therefore, we neglect them.

Finally, the term $`H_{hions}`$ represents the coupling between the electrical carriers and the magnetic impurities. There are two contributions to this term: i) An antiferromagnetic exchange interaction between the spins $`𝐒_I`$ of the Mn2+ ions located at sites $`I`$ and the spins of the itinerant carriers. This interaction produces a long range ferromagnetic interaction between the Mn ions. And ii) an interaction between the carriers charge and the potential arising from the Mn dopants:

$$H_{hions}=\frac{J}{S}\underset{I}{}𝐒_IC_{I,\alpha }^+𝝈_{\alpha ,\beta }C_{I,\beta }+W\underset{I}{}C_{I,\sigma }^+C_{I,\sigma }.$$
(3)
Here the sum runs over the positions of the magnetic impurities, and $`𝝈`$ are the Pauli matrices. The origin of $`W`$ is Coulombic, coming from the different electronegativity of Mn and GaAs atoms, and from the screening of the impurities by the electrical carriers.The quantity $`JW`$ governs the hybridization of the magnetic impurities with the host bands, and the excess of charge on the impurities with respect the host semiconductor atoms. For instance, in the case of a perfect ferromagnetic ground state, the electrical charge is uniformly distributed for $`J`$=$`W`$. The value of $`J`$ in our model is connected with the exchange coupling $`J_{pd}`$ appearing in the continuum model Konig et al.  through $`J=SJ_{pd}/2a^3`$. The value of $`J_{pd}`$ is typically 0.05-0.1$`eVnm^3`$ Matsukura et al. (1998); Ohno (1998); Omiya et al. (2000); Okabayashi et al. (1998), yielding $`J`$1.4-2.8$`eV`$, which gives a value of $`J/t`$1.7-3.4. Local spin density calculations (LSDA) and supercell methods Sanvito et al. (2001) applied to GaMnAs predict a value of $`J`$ rather larger than the experimental one, partially because the LSDA calculations overestimate exchange coupling strengths. Lacking reliable experimental information on the value of $`W`$, we consider it as a parameter with values between $`0`$ and $`J`$. Summarizing the model: we restrict ourselves to the terms $`H_h`$ and $`H_{hions}`$, given by Eq. 2 and Eq. 3, respectively.

## III Description of the Monte Carlo Algorithm.

We have performed classical MC simulations on the Mn core spin angles in order to calculate the temperature, $`T`$, dependent magnetic phase diagram. The simulations are done in a $`N\times N\times N`$ cubic lattice with periodic boundary conditions, and $`N_{Mn}`$=$`xN^3`$ randomly located magnetic ions. The number of carriers is taken to be $`N_h`$=$`pN^3`$ .

The simulations are done as follows. First, we produce a realization of the disorder: the positions of the Mn ions. For this realization, we compute the expectation value of the hole spin orientation at every place,

$$𝐧_i=<C_{i,\alpha }^+𝝈_{\alpha ,\beta }C_{i,\beta }>.$$
(4)
At a given temperature, the Mn spins fluctuate thermically on a fictitious magnetic field not  proportional to the carriers spin density,

$$𝐁_f(I)=\frac{J}{S}𝐧_I.$$
(5)
We perform MC simulations on the Mn spin angles until the system is equilibrated, when we choose a particular core spin configuration, and recalculate the hole spin density. This cycle is repeated until convergence in physical quantities is attained. Calculating the hole spin density requires the diagonalization of the Hamiltonian Eq.(1) at each step of the cycle. $`𝐧_i`$ is evaluated by filling up the $`N_h`$ lowest energy levels of the Hamiltonian. This can be done because the Fermi temperature is much higher than any other temperature in the system. For instance, for $`x`$=10$`\%`$, $`p`$=1$`\%`$, and $`J`$=$`W=`$3$`t`$, the Fermi energy is near 0.5$`t`$, while the Curie temperature is $``$0.015$`t`$.

The CPU diagonalization cost grows as $`N^6`$, severely limiting the system’s sizes under study. By diagonalizing the Hamiltonian in real space, it is not possible to analyze systems bigger than $`N`$=6 Calderón and Brey (1998); Vergés et al. (2002). In order to study bigger systems, it is better to use an adapted basis rather than the atomic localized basis. In the simulation, we use the basis obtained by diagonalizing the hole Hamiltonian for a perfect ferromagnetic alignment of the Mn ion spins. For the updating of the electron spin density at each cycle of the simulation, it is not necessary to include all the eigenvectors of the extended basis. We have found enough to work with a reduced basis formed by the $`N_c`$ lowest energy eigenstates of this basis for each carrier spin orientation. The value of $`N_c`$ depends on the particular values of the parameters. In general, we find that the size of the required basis increases with the value of $`J`$-$`W`$, i.e. it increases with the degree of localization of the hole charge on the Mn ions. Typically we can obtain good convergence with relatively small values of $`N_c`$, most of the results presented in this paper have been obtained with $`N_c`$=57. Note that, in the extended basis we use, the square of the wave functions is not spatially uniform, as it would happen in a Bloch or plane wave basis, but it is modulated by the spatial fluctuations of the potential felt by the carriers. Using this extended basis and an appropriate cutoff $`N_c`$, we can perform MC simulations in systems with sizes up to N=16.

In the simulations, we calculate the thermal averages of different quantities: the absolute value of the Mn ions spin polarization,

$$M=\frac{1}{SN_{Mn}}|\underset{I}{}𝐒_I|,$$
(6)
the absolute value of the hole’s spin polarization,

$$m=\frac{1}{N_h}|\underset{i}{}𝐧_i|,$$
(7)
and the average value of the hole charge on the Mn ions,

$$Q_{Mn}=\frac{1}{N_{Mn}}\underset{I}{}\rho _I,$$
(8)
being $`\rho _i=<_\sigma C_{i,\sigma }^+C_{i,\sigma }>`$ the electrical hole charge on site $`i`$. In the previous expressions, the index $`I`$ runs over the Mn positions and $`i`$ over all lattices sites.

## IV CPA method.

We have performed mean field calculations for the present model within the framework of the Coherent-Potential-Approximation (CPA) Taylor (1967); Soven (1967); Elliott et al. (1974). In recent times, this formalism has come to be known as the Dynamical-Mean-Field-Approximation (DMFA) Georges et al. (1996), and we use both terms (CPA and DMFA) as synonymous expressions. Given the possible confusion from the use of different flavors and names under the common umbrella of mean field calculations, we devote this section to briefly describe our CPA formalism, pointing out the origin and nature of the approximations and, where appropriate, making contact with other approaches.

Tracing out the fermions at constant chemical potential and temperature, we end up with a probability distribution for the core spins given by the following Boltzmann weight:

$$𝒫(\left\{𝐒_i\right\})exp\left[\beta \mathrm{\Omega }(\{𝐒_i\})\right]$$
(9)
where $`\mathrm{\Omega }(\{𝐒_i\})`$ is the carriers grand canonical potential for a given core-spin configuration. The exact distribution $`𝒫(\{𝐒_i\})`$ minimizes the following thermodynamical potential, considered as a functional over spin distributions:

$$\left[𝒫\right]=\mathrm{\Omega }(\{𝐒_i\})+T\mathrm{log}(𝒫)$$
(10)
where averages are taken over core-spin distributions. This formulation offers a convenient starting point for approximate variational approaches.

The common basis of all mean field calculations is the minimization of Eq. 10 with a factorization ansatz for the spin distribution:

$$𝒫(\{𝐒_i\})=\underset{i}{}p\left(𝐒_i\right)$$
(11)
From the critical phenomena point of view, this factorization legitimizes the term mean field, implying that all criticality (if any) will be classical. Carrying out the minimization of Eq. 10 with the only simplification of the factorization could be termed the absolute mean field approximation. Unfortunately, this procedure leaves us with the task of calculating the grand canonical potential for fermions (i.e. the density of states) in a site-disordered background: a formidable problem requiring numerical treatments (see, for instance, Ref. Alonso et al. (2001) for such a treatment in the Double-Exchange model).

To cope with the previous difficulty, additional approximations are needed to obtain the average density of states for the disordered problem. We will describe two of them, namely, the Virtual-Crystal-Approximation and CPA treatment (see Ref. Elliott et al. (1974) for a review on disordered systems).

i) Virtual Crystal Approximation (VCA). The average density of states for the disordered system is replaced by that of the average Hamiltonian. That implies a translationally invariant system with site energies given by the following average (a $`2\times 2`$ matrix in spin space): $`\epsilon _{\alpha ,\beta }=x(W+J\frac{𝐒}{S}𝝈_{\alpha ,\beta })`$. The minimization of $``$ is now straightforward, leading to the following well-known expression for the transition temperature:

$$T_c^{VCA}=\frac{2}{3}xJ^2n_\sigma (\mu )$$
(12)
where $`n_\sigma (\mu )`$ is the density of states per site and spin at the Fermi level (degenerate fermions assumed). This scheme is often referred to as the Weiss mean field approach, and we will consider both expressions (Weiss mean field and VCA) as synonymous terms.

ii) Coherent Potential Approximation (CPA). The CPA is a time-honored approach for the site-disordered ”alloy” problem Taylor (1967); Soven (1967). It is known to produce good results, particularly in three dimensions Elliott et al. (1974). The averaged Green’s function is obtained from a translationally-invariant, effective medium, characterized by a (site-diagonal and energy dependent) self-energy $`\mathrm{\Sigma }_{\alpha ,\beta }(z)`$ in the following manner:

$$G(z)G^{^{CPA}}=(zH_{_{_{CPA}}})^1$$
(13)
where (Dirac notation used to emphasize the one-body character of the problem)

$$H_{_{CPA}}=t\underset{(i,j)}{}|i,\sigma j,\sigma \left|+\underset{i}{}\right|i,\alpha \mathrm{\Sigma }_{\alpha ,\beta }(z)i,\beta |$$
(14)
with sum over spin indices implied. Given a particular site-disorder distribution, the self-energy is obtained with a minimum scattering requirement: the substitution of an effective site by a real one should produce zero (local) scattering on the average Economou (1983). That is:

$$(\widehat{\epsilon }_i\widehat{\mathrm{\Sigma }})[\widehat{1}\widehat{G}_{ii}^{^{CPA}}(\widehat{\epsilon }_i\widehat{\mathrm{\Sigma }})]^1=0$$
(15)
where hats mean $`2\times 2`$ matrices in spin space, and the random quantities $`\widehat{\epsilon _i}`$ are given by:

$$\epsilon _{i(\alpha ,\beta )}=\{\begin{array}{cc}0,& probability=(1x)\\ W+J\frac{𝐒_i}{S}𝝈_{\alpha ,\beta },& probability=xp\left(𝐒_i\right)\end{array}$$
(16)
As usual, the minimization of the thermodynamic potential can be cast in the form of a self-consistency equation for the probability distribution:

$$p\left(𝐒_i\right)exp\left[\beta \delta \mathrm{\Omega }(𝐒_i)\right]$$
(17)
where $`\delta \mathrm{\Omega }(𝐒_i)`$ can be interpreted as the total change in the fermion’s grand canonical potential for embedding an impurity site in the effective medium:

$$\delta \mathrm{\Omega }(𝐒_i)=𝑑\omega f(\omega )\frac{1}{\pi }\mathrm{}\mathrm{log}det\widehat{M}(𝐒_i)$$
(18)
$`f(\omega )`$ is the Fermi-Dirac factor, $`\mathrm{}`$ means imaginary part, and $`\widehat{M}(𝐒_i)`$ is given by:

$$\widehat{M}(𝐒_i)=\widehat{1}\widehat{G}_{ii}^{^{CPA}}(\omega ^+)[\widehat{\epsilon }(𝐒_i)\widehat{\mathrm{\Sigma }}(\omega ^+)]$$
(19)
It is through eqs. 17 and 18 that the CPA can be connected with the DMFA (see, for instance, Ref. Auslander and Kogan (a)), an approach primarily intended for genuine many-body problems like the Hubbard model Georges et al. (1996). It has been recently applied to models of itinerant carriers with ion impurities, believed to be relevant for manganites and doped magnetic semiconductors (see, for instance, Ref. Furukawa (1995, 1999); Auslander and Kogan (a); Chattopadhyay et al. (2001)).

Before leaving this section, let us comment on some features of these approaches. It is clear that the VCA is simple and handy, but the approximation in the density of states is too drastic. For instance, it predicts a ferromagnetic transition for any finite density of states. This cannot be true, even at the mean field level: for $`W=0`$ and small $`J`$, the problem becomes the standard Ruderman-Kittel-Kasuya-Yosida (RKKY) Hamiltonian, whose absolute mean field treatment shows ferromagnetism only for carrier density (holes or electrons) below a critical value $`p_c=0.5038`$ (sc lattice), due to the oscillatory behavior of the spin-spin coupling. Furthermore, the VCA is not sensitive to the value of $`W`$, obviously a major limitation for this problem. In contrast, the CPA is free from these shortcomings. For instance, it is able to reproduce the absolute mean field results of the RKKY limit Auslander and Kogan (b). In addition, extensive numerical results available for the exact density of states in the closely related Double-Exchange problem Alonso et al. (2001), show virtually no differences with the CPA density of states Gómez-Santos . This reinforces our confidence in this method for the present problem.

## V Results.

### V.1 Monte Carlo Results.

In Fig. 1 we plot the Mn ion’s spin polarization, $`M`$, as a function of $`T`$, for $`J`$=3$`t`$, and different values of $`W`$. The results correspond to a Mn concentration $`x`$=10$`\%`$, a hole concentration $`p`$=$`0.7\%`$, and a unit cell size $`N`$=14. In Fig. 2 we plot the hole’s spin polarization for the same set of parameters. Due to the finite size effects (note that there are only 274 Mn ions) the polarizations $`M`$ and $`m`$ are different from zero at any temperature, and we define the critical temperature, $`T_c`$, as the point where the second derivate of the polarization with respect to $`T`$ changes sign. This way of obtaining $`T_c`$ implies uncertainties of around 5$`\%`$. Within this error margin, the $`T_c`$’s obtained from the temperature dependence of $`m`$ and $`M`$ coincide. This is the expected result as the carrier’s spin polarization is due to the Mn polarization.

Comparing Fig. 1 and Fig. 2, we see that the ion spin polarization decreases faster with temperature than that of the holes. This seems to be generally true for the range of hole and Mn densities of interest, and a value of $`J`$=3$`t`$. This behavior has been already found in Weiss mean field calculations Jungwirth et al. (1999). Different disorder realizations are superimposed in figures 1 and 2. We find that, for small values of $`J`$-$`W`$, different disorder realizations give practically the same polarization curves, whereas, for small values of $`W`$, fluctuations related with disorder become more important. This is because the larger $`J`$-$`W`$, the deeper the effective impurity level. This increases the charge localized in the Mn ions, making the fluctuations associated with impurity randomness more important.

The critical temperature decreases when $`W`$ increases. As we will show, this is related to the amount of hole charge localized on Mn ions. In Fig. 3 we plot the average hole charge on Mn as a function of $`T`$, for two different values of $`W`$. Fig. 3 shows that the charge on the Mn’s decreases as the effective impurity potential increases, that $`Q_{Mn}`$ decreases with temperature, and that there is a change in the sign of the second derivate of $`Q_{Mn}`$ with respect to temperature at $`T_c`$.

The delocalization of charge with increasing $`T`$, and its abrupt change at $`T_c`$, indicate a possible change in transport properties of the system near the ferro-paramagnetic transition: we expect that the electrical resistivity of the system will decrease when $`T`$ increases. For a value of $`W`$=$`J`$, the effective impurity potential is exactly zero in the perfect ferromagnetic state. Close to this value, the charge on each Mn ion is practically $`p`$, independently of $`T`$. However, for values of $`W`$ near zero, the charge on the Mn ions is very large, more than three times $`p`$. Note that, in the case of $`W`$=0, the critical value of $`J`$ for binding a hole of the appropriate spin to a Mn ion is $`J_c=3.95677t`$. For values of $`J>J_c`$, the holes form a parallel-spin impurity band split from the main general band Chattopadhyay et al. (2001). In that case, the physics becomes more similar to that of the Double-Exchange materials. In our model, the parameter $`W`$ allows us to modulate the hybridization between the Mn ions orbitals and the host semiconductor band. In the case of $`J`$=3$`t`$, we believe that a lower impurity band is not formed for any positive value of $`W`$.

In Fig. 4 we plot the critical temperature obtained from MC simulations as a function of $`J`$ for different values of $`W`$, and for concentrations $`p`$=0.7$`\%`$ and $`x`$=10$`\%`$. We have also plotted the critical temperature, Eq. 12, obtained in the Weiss mean field theory Dietl et al. (2000, 2001). The Weiss mean field $`T_c`$ increases quadratically with $`J`$ and, as explained before, the approximation treats the magnetic impurities in the virtual crystal approximation, leading to a $`T_c^{VCA}`$ independent of $`W`$. From the MC results, we observe that the MC critical temperatures are smaller than the $`T_c^{VCA}`$’s only for values of $`W`$ very close to $`J`$. This is because, for $`W`$=$`J`$, the effective impurity potential is zero, and the virtual crystal approximation for obtaining the Curie temperature seems to work well. In this case, $`W`$=$`J`$ and small $`J`$, the thermal fluctuation correction to $`T_c^{VCA}`$ is similar to that occurring in the standard Heisenberg model, $`25\%`$. This becomes clear in Fig. 5, where we plot $`M(T)`$ as obtained in MC simulation and compare it with the predictions of Weiss mean field. Both curves are very similar at low $`T`$, and only near the Curie temperature do they separate, due to thermal fluctuation effects. For small values of $`W`$, the MC critical temperatures are larger than the $`T_c^{VCA}`$’s. This is due to the charge located on the Mn ions which increases the effective magnetic field felt by the Mn spins.

If the effective impurity potential $`J`$-$`W`$ increases, holes would eventually form a localized impurity band, leading to a decrease of $`T_c`$. As commented above, we do not get close enough to this limit, whose MC description would require a large basis beyond our present computational capabilities.

We have also studied the dependence of the critical temperature on the Mn concentration, $`x`$. This is shown in Fig. 6, where we plot $`T_c`$ as a function of $`x`$ for $`J`$=3$`t`$, $`W`$=2$`t`$, and $`p`$=0.7$`\%`$. Within the error margin, the critical temperature increases linearly with $`x`$, as expected for low impurity concentration.

Finally, we have studied the dependence of $`T_c`$ on the hole charge density $`p`$. In Fig. 7 we plot the MC results corresponding to the case $`J`$=3$`t`$, $`x=10\%`$, and two different values of $`W`$. For comparison, we have also plotted the Weiss mean field result. The Weiss mean field predicts a $`p^{1/3}`$ dependence, at low hole density. The MC results could be fitted perfectly to this dependence. Again we obtain that the Weiss mean field $`T_c`$’s are higher that the MC results only for values of $`W`$ near the value of the exchange coupling $`J`$.

Before finishing our discussion on the MC results, a comparison with previous simulations is in order. Recently, Schliemann and coworkers Schliemann et al. (2001) have presented some MC results for (III,Mn)V semiconductors. They work in a model similar to ours, but with the site energy difference set to zero, $`W`$=0. They describe the motion of the carriers using a plane wave basis with a cutoff $`k_c`$, and the exchange potential created by the magnetic impurities is modeled with a Gaussian of thickness $`a_0`$. They obtain that the critical temperature is substantially smaller than the estimated using Weiss mean field theory. This is opposite to our results: for $`W`$=0 our MC $`T_c`$ is higher than the mean field value. As the exchange constant and the densities (Mn and carriers) used in Ref. Schliemann et al. (2001) are similar to ours, there is a clear contradiction between our results and those of reference Schliemann et al. (2001). A possible explanation for this discrepancy can be found in the wavelength cutoff used in ref(Schliemann et al. (2001)), considerable larger than the impurity Gaussian width, and, perhaps, unable to trap carriers charge. If this were the case, their results could best correspond to our calculation for $`W`$=$`J`$, for which $`T_c`$ is smaller than $`T_c^{VCA}`$. In reference (Schliemann et al. (2001)), the wavelength cutoff is near 25Å, whereas the value of $`a_0`$ is smaller than 5Å.

A comparison with experimental results is strongly dependent on the values of the parameters. For instance, in the case of $`t`$=0.85$`eV`$, $`J`$=3$`t`$, $`x`$=10$`\%`$, and $`p`$=0.7$`\%`$, the critical temperature varies between 80K for $`W`$=$`J`$, and 350$`K`$ for $`W`$=0. This range would agree rather well with experimental data, assuming values of $`W`$ near $`J`$. Given the oversimplified nature of our model, this comparison with experiments cannot be taken too far away. In any case, our results indicate clearly the relevant role of value of $`W`$ on the magnetic transition: decreasing its value, if possible experimentally, would increase considerably the Curie temperature.

### V.2 CPA results.

In this section, we present results for the phase diagram obtained with the CPA method explained in IV. CPA results are parameterized by the density of states of the bare system, chosen to be that of the simple cubic lattice. Although unnecessary from the CPA point of view, we have assumed degenerate fermions. This is both physically correct for the relevant range of parameters, and makes the comparison with MC results more meaningful. We have verified that the inclusion of fermion temperature does not change results in any significant way.

In fig. 8 we present the critical temperature for $`J=3t`$ as a function of carrier concentration in the relevant region, for several values of $`W`$. As in the MC calculations (see Fig. 4), the most important effect of $`W`$ is the depression of $`T_c`$. Notice, though, that the overall temperature scale of the CPA calculation is lower than that of the MC results (compare with Fig. 7). This reverses the trend found in most statistical mechanics models, where the MC (exact) transition temperatures are higher that the mean field counterparts. We ignore whether this effect is a genuine peculiarity of our model or rather the outcome of approximations. We observe that $`T_c`$ follows the $`p^{1/3}`$ dependence with carrier concentration expected to hold close to a 3-d band in the Weiss mean field approximation, but this agreement between CPA and VCA is restricted to the small carrier concentration region (see below). In addition, the very apparent effect of $`W`$ on $`T_c`$ is completely missed in the VCA, as mentioned before.

$`T_c`$ as a function of Hund coupling $`J`$ is presented in fig. 9, for several values of $`W`$. The prominent effect of $`W`$ in reducing the critical temperature is very clear. Comparing with the MC data of Fig. 4, we see qualitatively similar curves, but, again, the MC temperatures are higher than the CPA ones. It is worth comparing the CPA results of fig. 9 with those of the VCA (Eq. 12), which predicts a dependence of the form: $`T_c^{VCA}xJ^2`$, irrespective of the size of $`J`$ and $`x`$. In the CPA treatment, this dependency only holds for small values of both $`J`$ and $`x`$, and significant deviations from the $`J^2`$ law are already evident in fig. 9, particularly for $`W=0`$. We have also studied the dependence of $`T_c`$ with impurity density, as shown in the inset of fig. 9 for small values of $`x`$. As in the MC calculations, an approximately linear regime is obtained in this limit. This linearity suggests that the effective magnetic interaction between core-spins could be well approximated by a two-body interaction in this regime. Notice that, in spite of the VCA prediction (Eq. 12), this linear dependence does not hold for higher impurity density, or even in the concentration range of fig. 9 for $`W=0`$, where non linear effects are already explicit. This behavior with $`J`$ and $`x`$ has been observed before for the present problem in a CPA calculation with $`W=0`$ and a model semicircular density of states Chattopadhyay et al. (2001). Its survival upon inclusion of disorder $`W0`$ and for the real density of states (sc lattice) confirms that it is a robust feature of the model.

It has been argued earlier that the effect of $`W`$ on $`T_c`$ can be correlated with changes in carrier density on the impurity sites. Being the magnetic coupling mediated by mobile carriers, decreasing their density on the impurity sites (with the repulsive term $`W`$) should lead to a decrease in $`T_c`$. We can test this behavior analyzing the limit of an almost filled band, with electrons as carriers (remember that our bare fermions are holes). The effect of $`W`$ would then be the opposite and, therefore, $`T_c`$ should increase. To verify this idea while gaining a better perspective of the problem, we have plotted in fig. 10 the phase diagram in the whole range of carrier concentration, for three values of $`W`$. For $`W=0`$, the problem ($`T_c`$) is particle-hole symmetric, and the ferromagnetic instability is restricted to the band edges. Although quantitatively far from the small $`J`$ limit, this confinement to band edges is qualitatively similar to the RKKY regime, where the same behavior can be traced back to the oscillatory nature of the magnetic coupling. As anticipated, fig. 10 confirms the asymmetric effect of $`W`$ on both types of carriers: it decreases $`T_c`$ at the lower edge (holes) and increases $`T_c`$ at the upper edge (electrons). Notice that, for electrons, $`W`$ acts as an attractive impurity, increasing the carrier density in the impurity sites, therefore increasing $`T_c`$. It is worth mentioning that fig. 10 also contains the CPA results for $`W<0`$: hole density $`p`$ and $`W>0`$ is equivalent to having chosen $`W<0`$ and hole density $`2p`$ (i.e. electron density $`p`$). From this point of view, the physics close to the upper end of fig. 10 upon increasing the value of $`W`$ is controlled by the development of a magnetic impurity band, with magnetic behavior becoming closer to that of the Double-Exchange model Chattopadhyay et al. (2001). Another interesting effect of increasing $`W`$ shows up in fig. 10: the spread of the magnetic instability from the lower edge into the band. We cannot offer a simple explanation for this fact.

Finally, we have investigated the phase diagram looking for the presence of phase separation (PS). PS is a common feature in magnetic models with itinerant carriers Nagaev (1997); Arovas and Guinea (1998); Dagotto et al. (1998); Arovas et al. (1999); Moreo et al. (1999). It seems to appear close to the magnetic transition temperature, when $`T_c`$ changes abruptly with carrier concentration Guinea et al. (2000). In order to search for PS, we explore the existence of global minima in the thermodynamic potential, in addition to the usual local minimum employed to obtain the curve of $`T_c`$. We observe small regions, close to the band edges,where the systems is unstable versus PS. This is shown in the insets of fig. 10 for two values of $`W`$. Comparing results for $`W=0`$ and $`W=2t`$, we see that the PS region decreases with increasing $`W`$, due to the associated reduction in the dependence of $`T_c`$ with carrier concentration. In the limit of small electron concentration, similar PS areas (not shown here) are found with opposite behavior with $`W`$.

As a final remark, we note that our calculation becomes qualitatively similar to a recent CPA analysis Chattopadhyay et al. (2001), carried out in the absence of impurity potential: $`W=0`$. Quantitatively, however, it is different (even for $`W=0`$) because of our using of the real density of states (sc lattice) as opposed to a model semicircular one. Nevertheless, close to the band edges, the differences between both calculations for $`W=0`$ amount to a different choice of lattice regularization, or, equivalently, a rescaling of the Hund coupling $`J`$.

## VI Summary.

We have studied a simple model for diluted magnetic semiconductors. In this model the carriers motion is described by a tight binding Hamiltonian and the hole-Mn ion coupling has two contributions: an antiferromagnetic exchange interaction between the Mn spin and the carrier spin, $`J`$, and an electrostatic interaction, $`W`$, between the carriers charge and the potential arising from the Mn dopants. We have studied this model performing Monte Carlo simulations and mean field CPA calculations. Our main results are the following:

1) The hole charge on the Mn ions has a relevant effect on the Curie temperature. The effective impurity potential $`JW`$ controls the hybridization of the Mn impurities with the host bands and the hole charge on the impurities. By increasing $`JW`$, if possible experimentally, the charge on the Mn ions and the critical temperature would increase considerably.

2) The Monte Carlo critical temperature is always higher than the corresponding CPA temperature. We also find that the Monte Carlo critical temperature is in general higher than the obtained from the simplest Weiss mean field theory. Only for values of $`JW`$ near cero the Weiss $`T_c`$ is higher than the Monte Carlo one.

3) The CPA calculation shows the existence of small pockets of phase segregation instability close to band edges.

## VII Acknowledgments.

We thank Prof. F. Guinea and Dr. J. Fernández Rossier for discussions. Finantial support is acknowledged from Grants No PB96-0085 (MEC, Spain) and CAM-07N/0008/2001 and CAM-07N/0015/2001 (Madrid, Spain)


text:
# Twisted Eguchi-Kawai Reduced Chiral Models

## 1 Introduction

A renewed interest in reduced matrix models has been motivated by the recent discovery of their relevance to noncommutative field theories and to superstring matrix models, see, e.g., Refs.  and references therein. Reduced matrix models were introduced for the study of the large-$`N`$ limit of $`SU(N)`$ gauge theories, see, e.g., Ref.  for a review. Eguchi and Kawai  pointed out that, as a consequence of the large-$`N`$ factorization, one can construct one-point theories equivalent to lattice $`SU(N)`$ gauge theories in the limit $`N\mathrm{}`$. The original proposal was just the replacement of all link variables $`U_\mu (x)`$ of the Wilson formulation of lattice gauge theories with four $`SU(N)`$ matrices according to the reduction rule $`U_\mu (x)U_\mu `$. This simple reduction procedure was shown to fail in the weak-coupling region, because of the spontaneous breaking of a symmetry that is crucial for the equivalence . Modifications of the original Eguchi-Kawai model have been then proposed . The most promising one is represented by the so-called twisted Eguchi-Kawai (TEK) reduction , which replaces $`U_\mu (x)D(x)U_\mu D(x)^{}`$, where $`D(x)=_\mu \mathrm{\Gamma }_\mu ^{x_\mu }`$ and $`\mathrm{\Gamma }_\mu `$ are appropriate $`SU(N)`$ matrices obeying the ’t Hooft algebra.

Since the proposal of the application of the TEK reduction prescription , the numerical checks  of the validity of the procedure have not been, to our opinion, completely exhaustive. We believe that, also in view of the renewed interest for matrix models, and particularly for TEK reduced versions, it is important to perform an accurate test of the conjectured equivalence with the corresponding lattice matrix model. This point may eventually be of relevance for the rôle of TEK models in noncommutative field theories or in superstring matrix (toy) models, where the theory should be analyzed in a limit different from the planar one.

In this paper we study the application of the twisted Eguchi-Kawai (TEK) reduction procedure to two-dimensional principal chiral models on the lattice . Two-dimensional chiral models and four-dimensional gauge theories present interesting analogies, see, e.g., Refs. . For example, they are both asymptotically free unitary matrix theories, and their large-$`N`$ limit is represented by a formal sum of planar graphs. In the large $`N`$ limit, one can apply the TEK reduction procedure to both lattice theories. In both cases the validity of the prescription is demonstrated by showing that the Schwinger-Dyson equations for the reduced and the non-reduced theory become formally the same in the limit $`N=\mathrm{}`$. In two-dimensional lattice chiral models a stringent test of the TEK reduction procedure should be simplified, for essentially three reasons: a better knowledge of the large-$`N`$ limit of the lattice chiral model, an easier way to define observables in spin models, and a simpler reduced model constituted by only one unitary matrix.

It has been recently realized that TEK reduced model may provide a nonperturbative formulation of noncommutative field theories . In particular, one may establish a correspondence between the TEK reduced model of the lattice principal chiral theory and a noncommutative U(1) principal chiral model . This is formally achieved in a continuum limit that keeps the noncommutativity parameter $`\theta _{\mu \nu }`$ fixed , where $`\theta _{\mu \nu }`$ scales as $`a^2N`$ and $`a`$ is the lattice spacing. Assuming the equivalence of the TEK reduced and the commutative principal chiral model in the planar limit, one may rewrite this condition as a double scaling limit involving the lattice coupling and $`N`$ . Therefore, our study of TEK reduced model in the planar limit is also useful to clarify the correspondence between TEK model and noncommutative field theory, and therefore the possibility of achieving a nonperturbative lattice formulation of a noncommutative field theory.

The plan of the paper is as follows. In Sec. 2 we introduce the original field theory and its reduced counterpart, and show on what grounds one expects that it coincides with the non-reduced theory in the large $`N`$ limit. In Sec. 3 we present our numerical results obtained by performing Monte Carlo simulations. Finally, Sec. 4 contains comments and conclusions.

## 2 TEK reduced chiral models

### 2.1 TEK reduced models

Twisted Eguchi-Kawai reduced models are based on the idea that, when $`N\mathrm{}`$, the $`\mathrm{SU}(N)`$ group becomes so large that it may accommodate the full Poincarè group as a subgroup, and allow us to find representations of the translation and rotation operators among the elements of $`\mathrm{SU}(N)`$. As a consequence, one may reformulate the full theory in terms of a finite number of matrix field variables defined at a single space-time site (or on the $`d`$ links emerging from the site in the case of a lattice gauge theory) and of the above-mentioned representations of the translation group. This reformulation is called “twisted Eguchi-Kawai” reduced version of the theory .

Consider a lattice field theory with fields $`\varphi (x)`$ in the adjoint representation of $`SU(N)`$. The twisted reduction prescription reads:

$$\varphi (x)D(x)\varphi D^{}(x),$$
(1)
where

$$D(x)=\underset{\mu }{}(\mathrm{\Gamma }_\mu )^{x_\mu },$$
(2)
and $`\mathrm{\Gamma }_\mu `$ are $`SU(N)`$ matrices obeying the so called t’Hooft algebra:

$$\mathrm{\Gamma }_\mu \mathrm{\Gamma }_\nu =Z_{\mu \nu }\mathrm{\Gamma }_\nu \mathrm{\Gamma }_\mu ,$$
(3)
where $`Z_{\mu \nu }`$ is an element of the center of the group, $`_N`$,

$$Z_{\mu \nu }=\mathrm{exp}\left(\frac{2\pi i}{N}n_{\mu \nu }\right)$$
(4)
and $`n_{\mu \nu }`$ is an integer valued antisymmetric $`d\times d`$ matrix (in $`d`$ dimensions). $`\mathrm{\Gamma }_\mu `$ is the matrix implementing the translation by one lattice spacing in the $`\mu `$ direction. The reduced action is simply obtained by substituing (1) into the action of the original lattice field theory:

$$S_{TEK}(\varphi ,n_{\mu \nu })=\frac{1}{\mathrm{vol}}S[D(x)\varphi D^{}(x)],$$
(5)
and the partition function reads:

$$Z_{TEK}=𝒟\varphi \mathrm{exp}(S_{TEK}),$$
(6)
for a fixed value of $`Z_{\mu \nu }`$. The expectation value of any functional of the reduced field $`\varphi `$ is given by:

$$O(\varphi )_{TEK}=\frac{1}{Z_{TEK}}𝒟\varphi O(\varphi )\mathrm{exp}(S_{TEK})$$
(7)
The correspondence between correlation functions of the reduced model and of the original field theory is as follows. Let $`f(x)`$ be any invariant functional of the field $`\varphi (x)`$. Then

$$f[\varphi (x)]_{\mathrm{field}\mathrm{theory}}=f[D(x)\varphi D^{}(x)]_{TEK}$$
(8)
The actual check of validity of the reduction procedure is based on the comparison of the Schwinger-Dyson equations of the reduced and original models. This procedure however requires some attention, since the limit of infinitely many degrees of freedom within the group itself allows the possibility of spontaneous breakdown of some of the symmetries which would be preserved for any finite value of $`N`$. Moreover, large $`N`$ is a thermodynamical limit: $`N`$ must go to infinity before any other limit is considered, and sometimes the limiting procedures do not commute.

### 2.2 Principal chiral models

The two-dimensional principal chiral field theory is defined through the following action

$$S=\frac{1}{T}d^2x\mathrm{Tr}\left[_\mu U(x)_\mu U^{}(x)\right],$$
(9)
where $`U(x)`$ are unitary matrix variables. The action is naturally invariant under the transformations $`UV_LU`$ and $`UUV_R`$ where $`V_L,V_RU(N)`$. As amply discussed in the literature (see, e.g., Refs.  and references therein), two-dimensional chiral models and four-dimensional lattice gauge theories manifest deep analogies, such asymptotic freedom, large-$`N`$ planar limit, etc…

A corresponding lattice principal chiral model can be obtained by the usual substitution of the derivative with a finite difference, obtaining

$$S=\beta N\underset{x}{}\underset{\mu =1,2}{}\mathrm{Tr}\left[U_xU_{x+\mu }^{}+U_{x+\mu }U_x^{}\right],$$
(10)
where $`\beta N1/T`$. In view of a large-$`N`$ analysis one may consider both $`SU(N)`$ and $`U(N)`$ models, since they are expected to reproduce the same statistical theory in the limit $`N\mathrm{}`$, which is meant at fixed $`\beta `$. Moreover, $`SU(N)`$ and $`U(N)`$ chiral models have the same continuum limit at any finite $`N2`$.

For a recent review of results concerning chiral models see, e.g., Ref. . Here we only mention the presence of a peculiar large-$`N`$ phase transition  which manifests itself with a power-law singularity in the $`N\mathrm{}`$ limit of the specific heat, i.e. $`C_{\mathrm{}}|\beta _c\beta |^\alpha `$, and with a peak in the specific-heat at finite $`N`$ that becomes sharper and sharper with increasing $`N`$ . The analysis of the $`N=\mathrm{}`$ strong-coupling series  provides the estimates $`\beta _c=0.3060(4)`$ and $`\alpha =0.27(3)`$, while the extrapolations to $`N=\mathrm{}`$ of the position of the specific-heat peak, determined by Monte Carlo simulations for $`SU(N)`$ and $`U(N)`$, gives $`\beta _c=0.3957(3)`$.

### 2.3 The reduced version: symmetries and equivalence

Implementing (1) into (10) gives the TEK reduced lattice chiral model action:

$$S_{TEK}=\beta N\underset{\mu =1,2}{}\mathrm{Tr}[U\mathrm{\Gamma }_\mu U^{}\mathrm{\Gamma }_\mu ^{}+h.c.].$$
(11)
In the present case, $`n_{\mu \nu }`$ of equation (4) is generically, in $`d=2`$, of the form

$$n_{\mu \nu }\left(\begin{array}{cc}0& M\\ M& 0\end{array}\right),M.$$
(12)
For a given $`N`$ and $`M`$ the solution to () is provided, up to global $`SU(N)`$ transformations, by the $`N\times N`$ shift and clock matrices

$$\begin{array}{c}S_{i,j}^{(M)}\delta _{i+M,j},\\ \\ C_{i,j}e^{\frac{2\pi i}{N}(i1)}\delta _{i,j}.\end{array}$$
(13)
The two matrices $`\mathrm{\Gamma }_\mu `$ are given respectively by $`S`$ and $`C`$.

While the original field theory is invariant under the $`SU(N)_L\times SU(N)_R`$ group, the reduced theory is invariant under the following two symmetries:

$`UzUz_N,`$ (14)
$`UD(x)UD^{}(x).`$ (15)
The first symmetry corresponds to the contraction of the original symmetry of the model to the center of the algebra of the symmetry group, while the second one is reminiscent of the translational invariance of the original theory.

In order to show the equivalence of the reduced theory with the original field theory in the large $`N`$ limit, the starting point is given by the Schwinger-Dyson equations, i.e. the equation obtained from a generic correlation function in the original theory (the reduced case is the obvious translation using Eq. (1))

$$G^{(n)}=\mathrm{Tr}\left(\underset{i=1}{\overset{n}{}}U_{x_i}^{k_i}\right),$$
(16)
by the following steps:

1. substitute $`U_{x_i}^{k_i}`$ with $`\lambda ^aU_{x_i}^{k_i}`$, where $`\lambda ^a`$ is a $`SU(N)`$ generator;
2. implement the change of variables $`U(1+i\alpha \lambda ^b)U`$ in the integral (implicitly taken in the average), using the Haar measure invariance;
3. contract the $`a,b`$ indices and extract the linear part in the parameter $`\alpha `$.

By comparing the two sets of equations obtained in the original and in the reduced theory, one finds manifestly corresponding terms and also some extra terms in the reduced theory. These extra terms exactly cancel if the symmetry $`UD(a)UD^{}(a)`$ of the reduced model is not spontaneously broken . Therefore, the two models possess the same Schwinger-Dyson equations if the above mentioned symmetry is not broken in any regime, and in particular in the continuum limit. The point is then whether this implies the equivalence of the two models, i.e. that they have the same expectation values of invariant quantities. Since the large $`N`$ equivalence of the twisted reduced models is based on the above arguments even for the lattice gauge theory case, it is an important task to verify, if the invoked symmetry holds, whether the Schwinger-Dyson equation based demonstration of the equivalence is indeed conclusive and sufficient or not.

The equivalence of the lattice chiral model and the reduced TEK model in the $`N=\mathrm{}`$ limit is not actually questioned in the strong-coupling region, i.e. for sufficiently small values of $`\beta `$. In this region the symmetries of the TEK reduced model are unbroken and therefore the equivalence of Schwinger-Dyson equations holds. One may then argue that this is sufficient to yield the equivalence of the reduced and nonreduced models. In the strong-coupling domain, i.e. within the convergence radius of the strong-coupling expansion, expectation values are analytic in the coupling $`\beta `$, and their boundary value at $`\beta =0`$ can be easily computed. As a consequence, it is formally possible to solve the Schwinger-Dyson equations in terms of strong-coupling series by iterative methods. Since the $`N=\mathrm{}`$ boundary values at $`\beta =0`$ of the lattice chiral model and the reduced TEK model coincide, equal Schwinger-Dyson equations imply the formal equivalence within the convergence radius of the strong coupling expansion 111Notice that the convergence radius of the strong coupling expansion does not necessarily coincide with the large-$`N`$ transition point $`\beta _c`$, but it may be smaller.. Problems may occur in the weak-coupling side of a large-$`N`$ phase transition.

It has been also argued that  the reduced TEK model is, in a sense, equivalent to a theory on a box of size $`L=N`$. In the large-$`N`$ limit the finite-$`N`$ corrections should be $`O(1/N^2)`$, just as in the $`SU(N)`$ lattice chiral theory. Since $`N^2=L^2`$, finite-$`N`$ corrections can be seen as finite volume corrections. Therefore, one may expect that the asymptotic large-$`N`$ regime is reached for $`N>>\xi `$, where $`\xi `$ is the $`N=\mathrm{}`$ correlation length at the given value of $`\beta `$. Note that in TEK reduced models the large-$`N`$ and thermodynamic limits are connected and approached simultaneously, while in the standard lattice models one first performs the $`L\mathrm{}`$ limit and then the large-$`N`$ limit. In a general field theory there are no a priori reasons why these two limits should commute.

## 3 Numerical results

In order to study the large-$`N`$ limit of the TEK model, we have performed Monte Carlo simulations for several values of $`N`$, up to $`N=200`$, and for values of $`\beta `$ in the strong and weak coupling regime. We used a Metropolis algorithm to update the $`SU(N)`$ matrix $`U`$. Trial matrices were selected by multiplying the actual matrix $`U`$ by a random SU(2) matrix embedded in $`SU(N)`$, choosing randomly among the $`N(N1)/2`$ $`SU(2)`$ subgroups. More precisely, once the $`SU(2)`$ subgroup is randomly chosen, we performed ten Metropolis hits with an approximate acceptance of 50%. Each $`SU(2)`$ updating requires $`O(N)`$ operations. The number of $`SU(2)`$-subgroup updatings per run was $`O(10^9)`$. We should also mention that for large values of $`N`$, $`N60`$ say, we observed some problem of thermalization when using completely random configurations (for example, constructed by multiplying $`N(N1)/2`$ completely random SU(2) matrix embedded in a $`N\times N`$ and associated to different subgroups) as starting point of our simulations. This occurred in both the strong- and weak-coupling regions, worsening with increasing $`\beta `$. So, as starting point we used either moderately random matrices or the unity matrix. It is well known that a simple Metropolis algorithm does not provide a particularly efficient method to simulate a statistical system. Our choice was essentially due to the fact that the reduced action is quadratic in the matrix variable $`U`$. Moreover, it does not lend itself to a linearization by introducing new matrix variables, as in the case of the reduced TEK gauge theory .

In our simulations we measured the internal energy

$$E=1\frac{1}{N}\mathrm{Tr}U\mathrm{\Gamma }_\mu U^{}\mathrm{\Gamma }_\mu ^{},$$
(17)
the quantity corresponding to the two-point function $`G(x)=\mathrm{Tr}U(0)U(x)^{}`$ using the equivalence (8), and in particular its low-momentum components, such as the (zero-momentum) magnetic susceptibility

$$\chi =\left|\mathrm{Tr}U\right|^2.$$
(18)
In order to check for possible spontaneous breaking of the symmetries (14) and (15), we considered the two quantities

$`S_1={\displaystyle \frac{1}{N}}\mathrm{Tr}U,`$ (19)
$`S_2={\displaystyle \frac{1}{N}}\mathrm{Tr}U\mathrm{\Gamma }_1U^{}\mathrm{\Gamma }_2^{}.`$ (20)
Indeed, $`S_1`$ is not invariant under the transformation (14) (but it is invariant with respect to (15)), while $`S_2`$ is not invariant under the transformation (15) (but it is invariant with respect to (14)).

Information on the large-$`N`$ limit of lattice chiral models in the strong coupling region can be obtained from the $`N=\mathrm{}`$ strong-coupling expansion and extrapolation of Monte Carlo results for the $`SU(N)`$ and $`U(N)`$ unitary groups. $`N=\mathrm{}`$ strong-coupling series are reported in Refs. . For example, the internal energy is known to 17th order and the magnetic susceptibility to 15th order. Monte Carlo simulations for relatively large values of $`N`$ (up to $`N=30`$) can be found in Refs. . It is useful to consider both $`SU(N)`$ and $`U(N)`$ data because their expectation values aproach the $`N=\mathrm{}`$ values from opposite sides, and allow us to have a better control of the large-$`N`$ limit.

As already mentioned, the equivalence of reduced and nonreduced models is not questioned in the strong-coupling region. However, Monte Carlo simulations in the strong coupling region are still useful to study the convergence of the large-$`N`$ limit of TEK, which is expected to be controlled by powers of $`N^2`$.

After verifying a general agreement between the large-$`N`$ TEK results and the $`N=\mathrm{}`$ strong-coupling series for small $`\beta `$, we have performed longer runs for $`\beta =0.28`$ in order to obtain an accurate test. This value of $`\beta `$ is not too far from the large-$`N`$ transition point $`\beta _c0.306`$. The large-$`N`$ limit of $`E`$ and $`\chi `$ of the lattice chiral model can be estimated from the analysis of the corresponding strong-coupling series , obtaining 222 Following Ref. , we used Dlog-Padé approximants to resum the series of $`dE/d\beta (\beta _c\beta )^\alpha `$ and $`d\mathrm{ln}\chi /d\beta (\beta _c\beta )^\alpha `$ with $`\beta _c0.306`$ and $`\alpha 0.27`$. The internal energy and the magnetic suceptibility are then obtained by integrating with respect to $`\beta `$. $`E_{\mathrm{}}=0.65370(4)`$ and $`\chi _{\mathrm{}}=7.02(1)`$. Extrapolations of the Monte Carlo results reported in Refs. , using the fact that the corrections are expected to be powers of $`N^2`$, provide consistent results, i.e. 333 For example, the simplest extrapolation using a straight line $`a+bN^2`$ passing by the data for the largest values of $`N`$ (i.e. $`\chi (N=30)=7.032(5)`$, $`\chi (N=21)=7.069(5)`$ for $`SU(N)`$ and $`\chi (N=21)=6.972(8)`$, $`\chi (N=15)=6.924(7)`$ for $`U(N)`$) gives $`\chi (N)=7.00(1)+32(6)N^2`$ for $`SU(N)`$ and $`\chi (N)=7.02(2)22(5)N^2`$ for $`U(N)`$.  $`\chi _{\mathrm{}}=7.00(2)`$ for $`SU(N)`$, $`\chi _{\mathrm{}}=7.02(2)`$ and $`E_{\mathrm{}}=0.65373(4)`$ for $`U(N)`$. Moreover, we mention that the large-$`N`$ limit of the second-moment correlation length is $`\xi _{\mathrm{}}1.53`$. In Table 1 we report $`E`$ and $`\chi `$ as obtained by our simulations for several values of $`N`$, up to $`N=150`$. In Figs. 1 and 2 we plot them versus $`1/N^2`$, which is the order of the expected asymptotic corrections. There we also show the MC data obtained from simulations of $`SU(N)`$ and $`U(N)`$ lattice chiral models. The figures show clearly that the correct $`N=\mathrm{}`$ limit is approached. To be more quantitative, fitting the data for $`N30`$ allowing for a $`O(N^2)`$ correction, we obtained

$`E0.65377(3){\displaystyle \frac{5.46(7)}{N^2}},`$ (21)
$`\chi 7.02(2)+{\displaystyle \frac{2.24(4)\times 10^3}{N^2}},`$ (22)
with acceptable $`\chi ^2/d.o.f1`$. This analysis shows that the correct large-$`N`$ limit is approached with corrections of $`O(N^2)`$, as expected. Note, however, that the $`O(1/N^2)`$ corrections are much larger than those found for the lattice chiral models, by about a factor $`10^2`$ for $`\chi `$, see the footnote 3.

As already mentioned, lattice unitary chiral models present a peculiar large-$`N`$ phase transition, which is characterized by a power-law singularity of the specific heat at $`N=\mathrm{}`$ . At finite $`N`$, numerical simulations show that the peak in the specific-heat becomes sharper and sharper with increasing $`N`$ . A peak in the specific heat is also observed in the TEK reduced model: it becomes sharper and sharper with increasing $`N`$. The $`N\mathrm{}`$ extrapolation of its location provides the estimate $`\beta _c0.306`$, in agreement with the results for the lattice chiral models.

These results provide an accurate check that for $`\beta <\beta _c`$ the TEK model shows the expected behavior, i.e. it converges to the same large-$`N`$ limit of lattice chiral model, with corrections of $`O(N^2)`$. We should also observe that, due to the large corrections, the convergence to the $`N=\mathrm{}`$ limit is rather slow, much slower than the convergence of the lattice chiral models. Comparing with the numerical effort to reach a similar accuracy in simulations of lattice chiral models, the reduced TEK model does not appear to be convenient to numerically investigate the large-$`N`$ limit of chiral models, at least in the strong-coupling domain.

Let us now turn to the weak-coupling region, i.e. for $`\beta >\beta _c`$. In this case information on the large-$`N`$ limit of lattice chiral models comes only from the MC simulations of Refs. . We performed runs for $`\beta =0.31`$ and $`\beta =0.32`$, for which the $`N=\mathrm{}`$ correlation length are estimated to be $`\xi 4`$ and $`\xi 6`$, respectively. MC simulations up to $`N=200`$ for these values of $`\beta `$ should provide a convincing check of the large-$`N`$ convergence, which should be controlled by the ratio $`N/\xi `$ according to the arguments mentioned in Sec. 2.3. The results of the simulations for $`E`$, $`\chi `$, and the two quantities $`S_1`$ and $`S_2`$, cf. Eqs. (19) and (20), are reported in Table 2. They should be compared with the the extrapolations to $`N\mathrm{}`$ of the MC results reported in Ref. , that are $`\chi _{\mathrm{}}=34.1(2)`$ and $`E_{\mathrm{}}0.519`$ for $`\beta =0.31`$, and $`\chi _{\mathrm{}}65`$ and $`E_{\mathrm{}}0.485`$ for $`\beta =0.32`$. Our MC data are also plotted in Figs. 3 and 4, where we also show the MC results for $`SU(N)`$ and $`U(N)`$ lattice chiral models, taken from Refs. . While the data for the energy seem somehow to approach the corresponding $`E_{\mathrm{}}`$ (although it is not evident that the corrections are $`O(N^2)`$), those for $`\chi `$ do not show any evidence of convergence to $`\chi _{\mathrm{}}`$. $`\chi `$ appears to increase with $`N`$ up to $`N=200`$. Similar results are obtained for $`\beta =0.32`$. We have also verified that the reduced two-point function does not exponentially decrease at large distance, explaining the behavior of its space integral $`\chi `$. Finally, note that the results for $`S_1`$ and $`S_2`$, see Table 2, are always consistent with zero, so that they do not give any indication of breaking of the symmetries (14) and (15).

Our results apparently contradict earlier numerical works , which were based on a much smaller statistic, and where the comparison was essentially restricted to the internal energy and the quantity $`S_1`$.

## 4 Conclusions

We have investigated the application of the twisted Eguchi-Kawai reduction procedure to two-dimensional principal chiral models on the lattice, as general testcase of the TEK reduction procedure.

In the strong coupling region, i.e. for $`\beta <\beta _c`$, the large-$`N`$ limit of the reduced theory appears to reproduce the results available for the original field theory in the large-$`N`$ limit, and that the TEK reduced model approaches the large-$`N`$ limit with power corrections of $`N^2`$. On the other hand, the validity of the reduction prescription apparently breaks down for $`\beta >\beta _c`$. In particular, it is not supported by the MC results for the magnetic susceptibility, that do not show evidence of convergence up to $`N=200`$ and for $`\beta `$ values corresponding to large-$`N`$ correlation lengths $`\xi 4`$ and $`\xi 6`$ (estimated extrapolating MC data for lattice chiral models). Of course, we cannot exclude that the correct asymptotic regime is eventually set for larger values of $`N`$, $`N10^2`$ say. Although we note that already for smaller values of $`\beta \beta _c`$ the approach to the large-$`N`$ limit is rather slow, in order to explain our results we should suppose a strong slowing down of the convergence with increasing $`\beta `$. In this regard, we should also add that the equivalence of the theories in the continuum limit requires a uniform convergence to $`N=\mathrm{}`$ when approaching the continuum limit. Another possibility is that our numerical simulations have been trapped by nontrivial extrema of the action, which has nothing to do with the correct limit related to the chiral field theory. However, our numerical checks definitely exclude it at least for $`N60`$.

We have also checked for possible spontaneous breakings of the symmetries (14) and (15) of the TEK reduced model, by computing the quantities $`S_1`$ and $`S_2`$, cf. Eqs. (19) and (20), which are respectively not invariant with respect to the symmetry (14) and (15). Their expectation values appear to vanish, thus suggesting that the symmetries of TEK are not spontaneously broken, although one cannot exclude breaking of symmetries associated with more complicated order parameters. We recall that one of such symmetries, i.e. the symmetry (15), is essential for the demostration of the equivalence of the $`N=\mathrm{}`$ Schwinger-Dyson equations. Assuming that the Schwinger-Dyson equations become formally equal for $`N=\mathrm{}`$, the question to be addressed is whether this is sufficient to identify a theory. Our numerical results suggest that this could not be the case, under the hypothesis that one can recover the $`N=\mathrm{}`$ theory from the finite $`N`$ behaviour of the reduced theory in its $`N\mathrm{}`$ limit. In other words, the $`N=\mathrm{}`$ TEK reduced theory may be still equivalent to the large-$`N`$ limit of lattice models, but the limit is singular, so that it cannot be reached from finite values of $`N`$. This issue is of relevance in any application of the TEK reduction procedure, as in the case of the lattice gauge theory.

The apparent failure of the TEK reduced theory to reproduce the planar limit of lattice chiral models does not necessarily prejudice the possibility of considering the TEK theory system in a different limit, possibly relevant to the fields where twisted Eguchi-Kawai reduced models were recently applied: noncommutative field theories and superstring matrix models. Indeed, TEK reduced chiral models with symmetry group $`U(N)`$ are formally equivalent to a U(1) noncommutative theory . This equivalence requires a nontrivial continuum limit that is different from the planar limit we have considered in this work, since one needs to keep the noncommutativity parameter $`\theta _{\mu \nu }`$ fixed, and $`\theta _{\mu \nu }`$ scales as $`a^2N`$ where $`a`$ is the lattice spacing. Thus, the relevant continuum limit should be realized for $`a0`$ and $`N\mathrm{}`$ keeping $`a^2N`$ fixed. The lattice spacing is somehow controlled by the coupling $`\beta `$, which should be appropriately tuned to realize the limit $`a0`$. The relation between $`a`$ and $`\beta `$ is therefore essential to define the continuum limit associated with the noncommutative field theory. However, the apparent failure of TEK model to reproduce the planar limit of the principal chiral theory does not allow us to use the known asymptotic behavior of the commutative theory in the continuum limit, for which asymptotic freedom and renormalization group imply $`a\mathrm{exp}(8\pi \beta )`$ in the continuum limit $`a0`$ and for $`N=\mathrm{}`$. Of course, this point deserves further investigation.

###### Acknowledgments.

We thank Luigi Del Debbio, Andrea Pelissetto, and Paolo Rossi for useful discussions.