On the reliability of recent Monte Carlo studies of dilute systems of localized spins interacting with itinerant carriers

In this paper, we discuss magnetic properties of dilute systems of localized spins interacting with itinerant carriers. More precisely, we compare recent available Monte Carlo results with our two step approach (TSA) calculations. The TSA consists first on the determination of the magnetic couplings and then on a proper treatment of the resulting effective dilute Heisenberg Hamiltonian. We show important disagreement between the Monte Carlo results (in principle exact) and our TSA calculations. We analyze the origins and shed light on the reasons of those dissensions. In contrast to one could expect, we demonstrate that the available MC calculations suffer from severe numerical shortcomings. More precisely, (i) the finite size effects appear to be huge in dilute systems, (ii) the statistical sampling (disorder configurations) was far too small, and (iii) the determination of the Curie temperature was too rough. In addition, we provide new arguments to explain a recent disagreement between the Monte Carlo simulations and TSA for the model study of the well known III-V diluted magnetic semiconductor Ga$_{1-x}$Mn$_{x}$As. We hope that this work will motivate new systematic large scale Monte Carlo calculations.


I. INTRODUCTION
During the last decade, the study of the effects of disorder on magnetic properties (Curie/Néel temperatures, magnetic excitations spectrum, spin susceptibility,...) of disordered/dilute magnetic materials became a central field of research for both theoreticians and experimentalists. It is now clear that even for a qualitative understanding of the underlying physics, disorder, dilution and thermal fluctuations should be handled in a proper way. Among widely studied materials one finds oxides manganite as for example La 1−x Sr x MnO 3 or La 1−x Ba x MnO 3 , which exhibit the fascinating giant magneto-resistance phenomenon 1,2,3,4 . The diluted magnetic semiconductors (DMS) as the well known III-V compound Ga 1−x Mn x As 5, 6,7,8 or the II-VI Zn 1−x Cr x Te 9 have attracted a considerable attention. The interest for DMS is partly motivated by their technological potential for spintronics. One can also mention the so-called d 0 materials as the wide band gap oxides HfO 2 , ZrO 2 , CaO, ZnO, or the irradiated graphite layers,... 10,11,12,13,14 . These new materials may become an alternative for spintronic applications. The Heusler alloys as Ni 2+x Mn 1−x Sn 15,16 or the double perovskites as Sr 2 FeMoO 6 17,18 are also materials in which disorder clearly plays a crucial role.
Ab-initio based studies have provided the most efficient and reliable tool, to allow theoretical and quantitative study of the magnetic properties without adjusting parameters. As an example, it was possible to study in great details and quantitatively the magnetic properties of the III-V diluted magnetic semiconductor Ga 1−x Mn x As, both in the presence or absence of additional compensating defects 19,20,21 . In order to test the interest of a particular material for technological applications, the ab initio approach provides the most appropriate tool. However, because it is material specific and very complex, it does not allow to draw general conclusions which could be valid for a whole family of materials. The study of relevant minimal model is crucial and allows to fill this gap. The model approach is essential to understand the influence of a particular physical free parameter. It can also help for the determination of potential spintronic candidates. The minimal one band model that describes itinerant carriers (holes/electrons) interacting with localized moments reads, In the first term, t ij =t if i and j are nearest neighbor, otherwise it is 0. The random variable p i is one if the site is occupied by a magnetic impurity, otherwise it is 0, we denote x the concentration of impurities. S i denotes the localized impurity spin at site i. For most of the materials S is large enough, we consider here the case of classical spins only. Moreover in ferromagnets the quantum fluctuations are not relevant. In the case of manganites p i = 1 (one Mn per site) J i is the Hund coupling (J H ) between the itinerant e g carriers and the localized t 2g magnetic moments (S = 3/2). In diluted magnetic semiconductors, J i is the local p-d coupling (J pd ). The additional last term describes the effects of the disorder induced by the substitution of a cation by another one. Note that this term includes in particular the electrostatic contribution which originates from the difference of charge between the substituted cation and the original one.
In this manuscript, all the calculations are performed on a simple cubic lattice. We also neglect in eq.(1) the last term. The case of finite on-site potential will be discussed elsewhere 22 . We expect the more realistic case of a multi-band Hamiltonian to lead to similar conclusions to those that will be presented in this manuscript.

II. THE TWO STEPS APPROACH PREDICTIONS
Reliable calculations requires several crucial conditions to be fulfilled. In other words, one should treat (i) the local coupling in a non perturbative way, (ii) the thermal and transverse spin fluctuations properly, and (iii) the disorder and dilution effects should be treated exactly (beyond an effective medium approach). The last condition allows for the localization of the itinerant carriers. This will strongly affect both magnetic and transport properties. The full Monte Carlo method is a priori the best procedure, since localized spins degrees of freedom and itinerant carriers are treated exactly, simultaneously and on equal footings. Note that the word full is added to avoid confusion with the Monte Carlo treatment of the effective Heisenberg Hamiltonian in the two step approach procedure presented below. An alternative approach which fulfill these requirements is the Two Step Approach (TSA). TSA is in spirit similar to what is widely used in the ab-initio based studies. It consists in two different stages. First, for a given configuration of disorder, at T = 0 K, one diagonalizes the Hamiltonian given by Eq. (1).This provides the spin resolved one particle Green's Function of the itinerant carrier (GF) G σ ij (ω) (where σ =↑ or ↓). In this first step the appropriate magnetic texture for the localized spins is used. For the ferromagnetism, all localized spins are parallel. Note that, the calculations are done without any approximations at T = 0 K, thus vertex corrections and the multiple scattering of the carriers on the magnetic impurities are included exactly. From G σ ij (ω) one then calculate the magnetic couplings between the localized moments (integrate out the carriers degrees of freedom). We end up with the disordered/dilute Heisenberg Hamiltonian which reads, where the magnetic couplings are given by the well known expression 23,24 In the present case ∆ i =∆ j =J pd (exchange splitting). The second step of the TSA consists in diagonalizing reliably the effective disordered Heisenberg Hamiltonian (Eq. (2)). Note that the diagonalization can be performed either by a Monte Carlo treatment (as mentioned above) or within the semi-analytical method called Self-Consistent Local RPA (SC-LRPA) 20 . Once again, to avoid any confusion, we underline and distinguish between the full MC study (no mapping to Heisenberg Hamiltonian) and the Monte Carlo that can be used to diagonalize Hamiltonian (2) in the TSA. It is important to recall that TSA have already been successfully used for the study of various realistic materials. In these studies, the first step of the TSA was provided by ab initio calculations (realistic exchange integrals). The agreement between calculated Curie temperatures and the experimental data was excellent 25,26,27,28 . Concerning the second step, it has also been shown that for a given set of exchange integrals, the SC-LRPA is as accurate as Monte Carlo calculations 29,30,31 . As it will be seen in the following, one great advantage of the TSA is the very large system sizes that can be handled, but are still inaccessible within the full Monte Carlo method. Thus, in this manuscript the second step of the TSA is performed within SC-LRPA. x = 1) and non disordered system (V = 0), in the double exchange limit (JS=∞).

A. Results for non dilute systems
Before discussing the case of dilute magnetic systems, it is important to discuss the non dilute case, e.g x = 1. A comparison between the full Monte Carlo treatment and the two steps approach has been done in the limit JS= ∞ (double exchange (DE) limit) 30 . In this case, the nearest neighbor coupling between localized moments reads J ij S 2 = t ij c † i c j , the spin index σ is irrelevant in this case, the fermion operators correspond to the majority spin. The comparison was performed in the presence of Anderson disorder, e.g the variables ǫ i in Eq.(1) were chosen randomly in the interval [−V/2,V/2], V denotes the strength of the disorder. Even in the presence of disorder an excellent agreement between the full Monte Carlo calculations and the two step approach was achieved, the difference in the Curie temperature was within 10 %. We now show that finite size effects are very small for the non dilute case. This will explain the agreement found between TSA calculations performed on large systems and full Monte Carlo simulations done with much smaller clusters. In Fig. 1, the nearest neighbor magnetic coupling (J 1 ) between impurities as a function of 1/N (N is the number of sites) is plotted. It is clearly observed that J 1 depends very weakly on the system size. The finite size effects are negligible at any carrier density (p h = γx).
In Fig. 2, we have plotted the Curie temperature as a function of the carrier density per magnetic moment (γ). We recall that in the clean limit the local SC-RPA reduces to the standard RPA calculation. In the case of classical spins the Curie temperature is given by +cos(q y a)+cos(q z a))). We clearly observe that the finite size effects are completely negligible even for very small systems. For the smallest system (N = 4 3 ), the Curie temperature is already accurate within less than 2%. Thus, in the case of clean systems and in the double exchange regime (JS=∞), the finite size effects are indeed completely neg-ligible. Similar conclusions have been reached in the case where the on-site disorder (V = 0 in the previous Hamiltonian)is also included 22 .We have even obtained that the average over only few configurations of disorder is enough to provide an accurate value. Thus, one can understand the success of the full Monte Carlo simulations to provide accurate Curie temperature for the study of the non dilute DE model even though in these studies the systems sizes are typically of the order of N ≈ 5 3 sites. The restriction to such small clusters is in fact due to the very large CPU time and memory costs of the standard full Monte Carlo simulations. However, as it will be seen in the following subsection, for dilute systems, the use of too small system sizes as well as too weak statistical sampling affect the results in a dramatic way.

B. Results for dilute systems
For simplicity and to allow a direct comparison with the available full Monte Carlo data, we set the impurity concentration to x = 0.065. We will vary both, the amplitude of J (its sign is irrelevant for classical spins) and the carrier density. Fig. 3 shows the DOS for various coupling JS. We observe that the impurity band (IB) splits from the valence band for values of JS ≥ W . In the strong coupling regime, JS = 2W , the IB is completely separated from the valence band, the magnetic couplings are expected to be dominated by the double exchange mechanism. Fig. 4 shows the couplings as a function of the impurity distance, for two different carrier concentrations and different values of JS. Note that, the cal- culations were performed on a sufficiently large cluster (N = 20 3 ) to neglect the finite size effects. The discussion of the finite size effects will be done in a forthcoming paragraph. First, we observe in Fig. 4 that the couplings are all ferromagnetic and of relatively short range. For x = 0.065, the typical average distance between magnetic impurities is d ≈ 2.4a. It corresponds to the distance between impurities on the ordered underlying super-lattice.
For the smallest values of JS, the couplings are small at any distances, this should lead to very small Curie temperatures. As JS is increased, the nearest neighbor coupling increases strongly until it saturates whilst the range of the couplings is reducing (clearly seen in the insets). For very large JS (JS ≫ W), the couplings are essentially reduced to the nearest neighbor coupling only, as expected in the double exchange limit. The limiting value is expected to be independent of the carrier density and corresponds roughly to the kinetic energy of a sin-gle hole hopping between two sites. Note that there are also some additional contributions coming from one hole localized on three sites, 2 holes on three sites and so on. One can already say that, in the large JS limit, the Curie temperature is expected to vanish since the average distance between impurities is larger than the range of the couplings. In other words, the system is below the percolation threshold. For JS ≈ 2W (see DOS in fig.1), the couplings at distances near d are tiny, thus one expects very small Curie temperatures.
In Fig. 5, we have plotted the Curie temperature as a function of JS. The carrier density is the same as used in Fig. 4. We again recall that the Monte Carlo treatment of this Hamiltonian would lead to similar results as it was already shown in several publications 19,29 . For both values of γ, we observe that T C increases until it reaches a maximum, then it decreases until it vanishes beyond a critical value which depends on the magnetic impurity concentration. Note that the maximum is located at JS ≈ 0.5 for γ = 0.125 and JS ≈ 0.6 for γ = 0.25. The highest T C are respectively T max C (γ = 0.125) ≈ 3.10 −4 W and T max C (γ = 0.25) ≈ 5.10 −4 W . For small values of J (perturbative limit) one notices that within SC-LRPA T C ∝ J 2 , as expected from Eq.(3). In Fig. 6, we evaluate the role of the disorder and the importance of the thermal fluctuations. For that purpose, we compare T C calculated within SC-LRPA to the Mean Field Virtual Crystal Approximation (MF-VCA) value. We remind that the MF-VCA value is T MF −V CA C = 2 3 x i =0 J 0i . Within this approximation both thermal and transverse fluctuations and disorder are treated at the lowest order. In fact the disorder (dilution)is neglected and just appears as a trivial prefactor (x) in the expression of T C . In the limit of the very small J (inset of Fig. 6), we observe a good agreement between SC-LRPA and MF-VCA. Within the perturbative limit, and for relatively weak densities of carriers only, one could indeed expect such a behavior (similar results are observed for γ = 0.25 of Fig. 5a). As we increase JS, we observe important quantitative and qualitative differences, as Let us now show that T C is in fact mainly controlled by the typical coupling namely J(d), where d was defined previously. For x = 0.065 we remind that the average distance between impurities is d = d ≈ 2.4a. In Fig. 7 we have plotted the variation of this coupling as a function of JS/W . We clearly see that the shape of the calculated curve is very similar to that obtained in Fig. 5a. We find that the ratio R = TC J(d)S 2 is almost independent of JS and R ≈ 4. On the other hand, the non monotonic behaviour of J(d) is qualitatively different from the monotonic increase of J 1 as seen in the inset. The limiting value of J 1 (for JS → ∞) is expected to be close to that of the clean system as plotted in Fig.1. For large JS, J 1 dominates and leads to the unrealistic MF-VCA critical temperatures (see Fig. 6).

III. COMPARISON WITH AVAILABLE FULL MONTE CARLO CALCULATIONS
In this section, we propose to compare our results with recent full Monte Carlo simulations results available in the literature. We will see that the quoted full Monte Carlo results are questionable since suffering from severe numerical shortcomings. To our knowledge, the only available full Monte Carlo study of the diluted Hamiltonian given by eq. ison with our calculations (performed on a single band model). The MC data of ref. 32 are shown both in the insets of Fig. 5a and b. The MC critical temperatures were obtained for systems with N = 5 3 to 6 3 sites and the average over disorder was performed over 7 configurations of disorder only. Let us notice that a system of size N = 5 3 (resp. 6 3 ) contains only 8 (resp. 14) magnetic impurities; thus for a density of carriers per impurity γ = 0.25, such small systems respectively contain only 2 and 3 holes in the whole cluster. As seen in Fig. 5, several crucial differences appear immediately with respect to our results. First, we notice that for both carrier densities (γ = 0.125 and 0.25) the Monte Carlo simulations predict critical T C much higher than those obtained within the two steps approach. Indeed, if one refers for example to Fig. 5a, one observes that the MC values are 15 to 60 times larger! In particular, the highest T C are typically at least 20 times larger than our results.
For small values of JS, one notices that SC-LRPA and the Monte Carlo simulations do not predict the same behavior (quadratic in J in our case). For large JS, whilst our calculated T C is strongly suppressed and vanishes for values of JS ≥ J c S where J c S ≈ 2 W , the MC value decreases much slowly and seems to saturate at large JS (see Fig. 5a On the basis of what was previously discussed, we now address the reasons which may explain the huge differences with TSA and why the Mean Field VCA becomes a lower bound with respect to the full Monte Carlo results. To be more specific, it will be shown that the available Monte Carlo results suffer from several combined effects: (i) strong finite size effects, (ii) insufficient sampling (too few disorder configurations) and (iii) inaccurate and approximative procedure for the determination of the Curie temperature.

IV. ORIGINS OF THE DISSENSIONS BETWEEN TSA AND MONTE CARLO SIMULATIONS
In this section we will discuss both the importance of the statistical sampling and finite size effects. In Fig.8 a  and b (γ = 0.25 and 0.125) we have plotted the variation of J ij S 2 as a function of the distance between impurities, for different system sizes (N = L 3 where L varies from 4 to 16). We observe that the couplings are larger for the smallest systems and that finite size effects are especially strong at distances near the typical and relevant one, d. Thus, we already expect the Curie temperature to be strongly size dependent. In Fig.9 the average Curie tem- perature as a function of JS/W is shown for various sizes. The hole density (per carrier) is set to γ = 0.125. Note that each averaged T C is obtained using the associated couplings, e.g calculated for the corresponding system. Note also that the average was properly performed, we have used a sufficiently large number of disorder configurations (few thousands for the smallest systems and few hundreds for the largest). We analyze now the importance of the statistical sampling. For a given value of JS, we observe beyond N = 14 3 an insignificant variation of the averaged Curie temperature: beyond this size the thermodynamic limit is properly described. Because within standard full Monte Carlo method the systems are usually restricted to relatively small sizes 32,33,34,35,36 of the order of N = 4 3 , let us do the comparison between the calculations done with the smallest and largest systems. For instance, we consider JS = 0.7 W . We see that the Curie temperature for L = 5 is T C = 10.5 10 −4 W whilst in the thermodynamic limit T C = 2.0 10 −4 W (500 % difference!). We now show that the sampling over disorder is also crucial and should be done properly, especially in dilute systems. Note for instance that, Franceschetti et al. 37 , in their ab initio based study, have shown important fluctuations of the Curie temperature from sample to sample. In Fig.10 we have plotted the distribution of the Curie temperatures obtained using different system sizes. To facilitate the discussion, we have plotted on Fig.11 the averaged Curie temperatures, T L C = +∞ −∞ T C P L (T C )dT C , and the widths at mid height of the distributions (∆T C /W ) as a function of 1/N. For the smallest system (N = 4 3 ) one observes a very important spreading of the critical temperatures distribution: they can vary from at least one order of magnitude! From Fig.11, for this case, ∆T C /W is rather close to the corresponding average value. For these small systems, one understands easily that 10 configurations of disorder are definitely insufficient to provide a reliable average. In other words, the calculated averaged Curie temperature over a small number of configurations (few tens) could easily lead to critical temperatures 10-20 times larger than that calculated properly (averaged over few thousands of configurations). This implies, for example, that the Curie temperatures shown in Fig.9 could be easily 50 times larger than the calculated ones in the thermodynamic limit if the average would have been taken over only few configurations. As the system size increases, one gradually observes a decrease in ∆T C , which then rapidly tends to zero (thermodynamic limit), whilst T L C tends towards a constant value (see Fig.11). One can already consider that the thermodynamic limit value of T C is reached for L ≥ 16. In the light of the results presented in this manuscript, a criterion to get reasonable results would consist in considering system sizes where L is at least five to six times larger than the typical distance between impurities. This criterion is only rough, in the case of long range magnetic couplings, such as RKKY exchanges, the finite size effects can be even more drastic. In addition, only ten configurations of disorder for small systems are definitely insufficient to get the average value of T C in a reliable manner. It is obvious that these numerical requirements are very difficult to fulfill within standard Monte Carlo calculations but they are definitely essential. In addition, in the presence of inhomogeneities, temperatures are also expected to be even more sensitive to both finite size effects and statistical sampling 38 . Moreover, the way of extracting the Curie temperature is essential. Within SC-LRPA, the problem does not arise. The critical temperatures are directly given by a semi-analytical equation solved in a self-consistent way, no extrapolation from the magnetization curve is used. Concerning the Monte Carlo simulations of reference 32 , the method used to extract T C from the magnetization curve is not very accurate (see figure 4.a and 4.b of ref. 32 ). Potential errors related to this method are added to the uncertainties coming from both finite size effects and from the insufficiency of statistical sampling. A more accurate way to obtain the Curie temperatures in MC simulations is based on the method of Binder cumulants 39 , but this method represents an additional cost in terms of computing time since it needs a finite size effect analysis. Unfortunately, in spite of the small system sizes and the small number of configurations of disorder considered in ref. 32 , the Monte Carlo simulations already needs a huge resource in term of CPU time and memory.
In a recent paper based on Monte Carlo simulations, Yildirim et al. have studied the diluted Kondo model applied to the case of Ga 1−x Mn x As by including the realistic band structure of the host material 34 (see also ref. 41 for a comment). They have calculated the Curie temperature as a function of the local coupling J for systems containing x = 8.5% of magnetic impurities and p h = 0.75 hole per impurity. We underline here that their numerical Monte Carlo approach suffers from the same numerical shortcomings as those previously discussed. The system considered in this letter contain typically 20 localized spins only and the average are done using just five configurations of disorder. Additionally, the value of J considered in ref. 34 , corresponds to the perturbative regime. Thus, on the basis of the present study, we argue that within this limit, the couplings should exhibit RKKY oscillations which should lead in the thermodynamic limit to small Curie temperatures or eventually to no ferromagnetic phase 42 . Unfortunately, and as previously noticed, the smallness of the cluster considered in the MC calculations hide the effects of the asymptotic RKKY tail, and thus leading to finite and large Curie temperatures. We argue that by improving the statistics and increas-ing the size of the systems the Curie temperature should vanish when the frustration will become effective. Furthermore, we argue that the large MC critical temperatures found in the large JS regime (see Fig.1 of ref. 34 ) are in fact a numerical artefact due again to the insufficient statistical sampling and finite size effects. Indeed, in this regime double exchange mechanism dominates and thus very small or vanishing Curie temperatures are expected, see the corresponding densities of states in Fig.4 of 34 where the IB is separated from the valence band.

V. CONCLUSION
To conclude, in this work we have shown that the study of the diluted magnetic systems requires a rigorous numerical treatment. We have shed light on the origins of the disagreements between MC and TSA for the study of diluted magnetic systems. In contrast to the non dilute systems, the finite size effects as well as the importance of the statistical sampling appear to be crucial. We have shown that available Monte Carlo simulations for the diluted model Hamiltonian (1) suffer from severe numerical insufficiencies. Although Monte Carlo simulations are in principle exact, but because the calculations are unfortunately restricted to relatively small systems and weak statistical sampling, the obtained Curie temperatures are often largely overestimated. In addition, limiting regimes (strong and weak couplings) are not properly described. It would be of great interest to check both the importance of statistical sampling and the finite size effects within new large scale Monte Carlo studies. Among MC simulations which allow to reach relatively large system sizes, one can quote the Hybrid Monte Carlo method (HMC) 43 , the Polynomial Expansion Monte Carlo (PEMC) 44 or its fastest counterpart, namely the Truncated Polynomial Expansion Monte Carlo (TPEMC) method 45,46,47 . In Particular TPEMC handles higher systems up to 20 3 sites within a reasonable CPU time (see table II of ref. 46 ).