text:
# Monte Carlo simulations of an impurity band model for III-V diluted magnetic semiconductors

## I Introduction

III-V diluted magnetic semiconductors (DMS) have become a very active area of research due to their interesting magnetic and transport properties.  Thus far Ga1-xMnxAs has received the greatest amount of attention, due to the observation that it becomes ferromagnetic with a Curie temperature as high as 110K, when $`x=0.053`$. More recently, ferromagnetism above room temperature has been observed in (Ga,Mn)N.

The observed ferromagnetism is widely accepted to be due to a charge carrier mediated coupling between the Mn spins. Several models have been proposed to explain the detailed phenomenology of these compounds.  In particular, we have initiated an effort to understand the effects of disorder in Mn positions on the properties of these compounds. An additional source of disorder is the large degree of compensation seen in these materials, which has been attributed to As antisite defects,  because of which, the carrier density is significantly less than the Mn density.

To this end, we study the low Mn concentration regime, near and below the metal-insulator transition (MIT) at $`x0.035`$, where disorder effects would be expected to be the most pronounced. In this limit, we model the charge carriers in terms of an impurity band comprised of states around each Mn acceptor, which is taken to be the source of charge carriers mediating ferromagnetism. Evidence for the relevance of an impurity band has been provided by a number of recent experiments, such as an scanning tunneling microscope study  which showed the existence of an impurity band in (Ga,Mn)As samples with $`x=0.0050.06`$. Angle resolved photoemission spectroscopy in a sample with $`x=0.035`$ has also revealed a well formed impurity band, and confirmed that the Fermi energy lies in its vicinity.

A complete, detailed description of the impurity band in the presence of compensating defects is an extremely difficult enterprise,  hence we have proposed a simple effective Hamiltonian that captures (at least qualitatively) the relevant impurity band physics. This Hamiltonian was studied comprehensively at the mean-field level, and the magnetic properties were observed to have a number of surprising properties. The mean-field magnetization curves have very unusual, concave upward shapes, unlike the magnetization curves of conventional ferromagnets. Some of these features have been seen in experimental measurements, especially for samples with a low carrier density and high degree of compensation.  This unusual shape of the magnetization curves was identified to be a direct effect of positional disorder of the Mn ions.

The mean field calculation suggested that there is considerable inhomogeneity in the magnetization of individual Mn spins at temperatures below $`T_c`$, particularly for small values of $`x`$. Experimentally, disorder appears to be relevant even in the metallic phase – Barkhausen jumps were observed in a sample with $`x=0.047`$,  indicating the presence of frozen-in magnetic disorder. Theoretically, disorder has also been credited as leading to an instability in collinear Mn ground states in the metallic phase,  via spin waves, which can lead to a reduction in the saturation magnetization at low temperatures as observed in annealed samples with $`x=0.05`$.

Increased disorder was also shown to lead to an increased value of the critical temperature $`T_c`$ in the mean field study. It is well known that the mean-field approximation underestimates the effect of thermal fluctuations and therefore overestimates the value of the critical temperature. As a result, it is important to check the mean field results against Monte Carlo simulations which properly account for thermal fluctuations, to test to what extent the phenomenology found in the mean-field study is maintained when these fluctuations are included.

In this paper, we report the results of such a Monte Carlo study. We find that the magnetization curves retain their unusual shape even when thermal fluctuations are included and that, as expected, the critical temperature is lowered from the mean field value. The extent of this decrease is most pronounced for low values of $`x`$ and high compensation. The observation that disorder can lead to a higher critical temperature than a purely ordered case  is also confirmed, although the change is more modest than in the mean field study. The spatial inhomogeneity observed in the mean field studies  is confirmed and found to be little changed by thermal fluctuations. A correlation between regions of larger local charge density and larger local magnetization is also established.

This paper is organized in the following manner: in section II we describe the model Hamiltonian and the values chosen for various parameters of the system. We also discuss the use of finite size scaling to determine the ferromagnetic transition and describe the quantities that we calculate. In section III, we introduce a perturbative scheme to speed up the Monte Carlo simulations, describe our implementation and discuss testing on a toy model for which we can compare with exact results. We present the results of our Monte Carlo simulations for the impurity band model of DMS in sections IV and V. Section VI summarizes our conclusions and discusses the implications of our results for experiments and modelling of III-V DMS.

## II Model

### A Hamiltonian

When manganese is introduced into GaAs, the Mn impurities have been shown to substitute on the Ga FCC sub-lattice of the zinc-blende structure of the undoped semiconductor for small values of $`x`$.  However, at larger values of $`x`$ ($``$ 0.07), Mn ions can form MnAs clusters, which have the NiAs structure.  Based on these experimental findings, we assume a zinc-blende structure and Mn substituting only on Ga sites for the low Mn concetrations we study. Each Mn impurity has a spin-$`\frac{5}{2}`$ from its half-filled 3$`d`$ shell. The nominal valence-II of Mn implies that when it substitutes for the valence-III Ga, it acts as an acceptor. Thus an isolated Mn can bind a hole in an impurity level that we characterize by a hydrogenic orbital with Bohr radius $`a_B`$, with wavefunction $`\varphi _i(𝐫)\mathrm{exp}[|𝐫𝐑_i|/a_B]`$. In Mn doped GaAs there are substantial central cell corrections  which we phenomenologically incorporate by adjusting the Bohr radius $`a_B`$.

Whilst the true carriers in this system are holes with spin $`\frac{3}{2}`$, we consider the case of electron doping, so the carrier spin is $`\frac{1}{2}`$. This leads to some differences  (in particular, the frustration effects recently claimed for hole doping  are absent in our model). However, the effects of disorder and impurity potentials, which are our main focus here, should not be qualitatively changed.  The carriers (electrons) are then in an impurity band below the conduction band minimum of the host semiconductor.

We study the model introduced by Berciu and Bhatt,  for which the Hamiltonian is

$``$ $`=`$ $`{\displaystyle \underset{i,j}{}}t_{ij}c_{i\sigma }^{}c_{j\sigma }+{\displaystyle \underset{i,j}{}}J_{ij}𝐒_i\left(c_{j\sigma }^{}{\displaystyle \frac{1}{2}}𝝈_{\alpha \beta }c_{j\beta }\right)`$ (2)
$`g^{}\mu _B𝐇{\displaystyle \underset{i}{}}c_{i\sigma }^{}{\displaystyle \frac{1}{2}}𝝈_{\alpha \beta }c_{i\beta }g\mu _B𝐇{\displaystyle \underset{i}{}}𝐒_i.`$
The random positions of the Mn impurities are labelled by $`𝐑_i`$, and the spin of the Mn impurity at $`𝐑_i`$ is $`𝐒_i`$. In the electron formalism for the charge carriers, $`c_{i\sigma }^{}`$ is the creation operator for a carrier with spin $`𝝈`$ in the bound state associated with the $`i^{th}`$ impurity. The hopping matrix is given by $`t_{ij}=t(|𝐑_i𝐑_j|),`$ where $`t(r)=2(1+r/a_B)e^{r/a_B}\mathrm{Ry}`$,  and the Rydberg (Ry) is the binding energy of a hole. We assume that the Mn spins are strongly localized and hence the exchange integral is given by $`J_{ij}=J_0\left|\varphi _{ij}\right|^2=J_0e^{2|𝐑_i𝐑_j|/a_B},`$ which is proportional to the charge density at the $`i^{th}`$ Mn site of the carrier in the hydrogenic wavefunction around the $`j^{th}`$ Mn site. The external magnetic field is given by $`𝐇`$, the Landé g-factors for the Mn spins and carriers are $`g`$ and $`g^{}`$ respectively, and $`\mu _B`$ is the Bohr magneton. Our simulations are for zero magnetic field.

We consider finite size cubic samples with $`L`$ cubic unit cells of the GaAs structure per side, with periodic boundary conditions, and $`N_d`$ Mn impurities. Thus $`x`$ is given by $`N_d/(4L^3)`$, and is related to the concentration of Mn impurities, $`n_{Mn}`$ through $`n_{Mn}=4x/a^3`$, where $`a=5.65`$ Å is the GaAs lattice constant. The number of carriers is $`N_h=pN_d`$, with $`p`$ between 0.1 and 0.3, as indicated by experimental studies of samples grown by molecular beam epitaxy.  As mentioned previously, these low values of $`p`$ are due to compensation processes, in which As antisites are believed to play an important role.

The only difference between our model and that in Ref. , is that instead of studying spin-$`\frac{5}{2}`$ Mn spins, we treat the Mn spins as classical Heisenberg vector spins. This should be a reasonable approximation since $`S=5/2`$ is a large spin and a Quantum Monte Carlo calculation would not gain much due to the uncertainties in the materials parameters. \[At the mean field level, this approximation has the effect of lowering $`T_c`$ by a factor of $`1+\frac{1}{S}`$ relative to the $`T_c`$ for quantum spins\]. One quantity which differs significantly with classical spins is the specific heat at low temperatures,  which we have not studied.

The values we use for numerical parameters are $`a_B`$ = 7.8 Å, 1 Ry = 112.4 meV, and $`J_0`$ = 15 meV, as discussed in Ref. . With these parameters we have an impurity band whose bottom lies around 200-300 meV (2-3 Ryd) below the host conduction band and the Fermi energy varies from around 13 to 55 meV (depending on $`x`$ and $`p`$) above the bottom of the impurity band for the parameter range considered. These values are of the same order of magnitude as the observed splitting and bandwidth of the impurity band in angle resolved photoemission experiments.

Our goal in this work is partly to understand the deviations from mean-field theory in the model Eq. (2), therefore we also perform mean field calculations using Langevin functions to represent the polarization of classical spins (in the quantum case one uses Brillouin functions), with which we compare our Monte Carlo results.

### B Method of Simulation

Consider a system with classical spin degrees of freedom and fermionic degrees of freedom, such as described by the Hamiltonian (2). The assumption of classical spins means that one can parameterize the spin at each site by its $`z`$-component and azimuthal angle, i.e. $`𝐒_i=(S_i^z,\varphi _i)`$, and

$`S_i^x`$ $`=`$ $`S\sqrt{1{\displaystyle \frac{S_{i}^{z}{}_{}{}^{2}}{S^2}}}\mathrm{cos}\varphi _i,`$ (3)
$`S_i^y`$ $`=`$ $`S\sqrt{1{\displaystyle \frac{S_{i}^{z}{}_{}{}^{2}}{S^2}}}\mathrm{sin}\varphi _i.`$ (4)
For any given configuration of classical spins, $`\{𝐒\}=\{S^z,\varphi \}`$, the Hamiltonian (2) can be diagonalized to give

$$(\{𝐒,a_n^{},a_n\})=\underset{n}{}E_n(\{S^z,\varphi \})a_n^{}a_n,$$
(5)
where the states $`n`$ are a diagonal basis and

$$a_n^{}=\underset{i\sigma }{}\psi _{n\sigma }(i)c_{i\sigma }^{},a_n=\underset{i\sigma }{}\psi _{n\sigma }^{}(i)c_{i\sigma },$$
(6)
are linear combinations of the $`c_{i\sigma }^{}`$ and $`c_{i\sigma }`$ operators. Here $`\psi _{n\sigma }(i)`$ is the eigenfunction for the $`n^{\mathrm{th}}`$ energy level at site $`i`$ for spin $`𝝈`$. Using the grand-canonical ensemble, the partition function is

$`Z`$ $`=`$ $`\left[{\displaystyle \underset{i=1}{\overset{N_d}{}}}{\displaystyle _1^1}𝑑S_i^z{\displaystyle _0^{2\pi }}𝑑\varphi _i\right]{\displaystyle \underset{n_1=0,1}{}}\mathrm{}{\displaystyle \underset{n_{N_0}=0,1}{}}`$ (8)
$`\times e^{\beta _{k=1}^{N_0}(E_n(\{S^z,\varphi \})\mu )n_k},`$
where $`n_k=0,1`$ is the occupation number of level $`k`$ (there are $`N_0=2N_d`$ levels), and $`\beta =1/(k_BT)`$ is the inverse temperature. Summing over fermion degrees of freedom leads to the result

$`Z`$ $`=`$ $`\left[{\displaystyle \underset{i=1}{\overset{N_d}{}}}{\displaystyle _1^1}𝑑S_i^z{\displaystyle _0^{2\pi }}𝑑\varphi _i\right]{\displaystyle \underset{n=1}{\overset{N_0}{}}}\left(1+e^{\beta (E_n(\{S^z,\varphi \})\mu )}\right).`$ (9)
This can be cast in a similar form to that used for a spin-only system

$`Z=\left[{\displaystyle \underset{i=1}{\overset{N_d}{}}}{\displaystyle _1^1}𝑑S_i^z{\displaystyle _0^{2\pi }}𝑑\varphi _i\right]e^{\beta _c(\{S^z,\varphi \})},`$ (11)
where the carrier free energy $`_c(\{S^z,\varphi \})`$ for a given configuration of classical spins is

$`_c(\{S^z,\varphi \})={\displaystyle \frac{1}{\beta }}{\displaystyle \underset{n=1}{\overset{N_0}{}}}\mathrm{log}\left(1+e^{\beta (E_n(\{S^z,\varphi \})\mu )}\right).`$ (12)
The chemical potential $`\mu `$ is determined from the condition

$`N_c`$ $`=`$ $`{\displaystyle \frac{}{(\beta \mu )}}\mathrm{log}Z,`$ (13)
$`=`$ $`{\displaystyle \frac{1}{Z}}\left[{\displaystyle \underset{i=1}{\overset{N_d}{}}}{\displaystyle _1^1}𝑑S_i^z{\displaystyle _0^{2\pi }}𝑑\varphi _i\right]{\displaystyle \underset{n=1}{\overset{N_0}{}}}f(E_n)e^{\beta _c(\{S^z,\varphi \})},`$ (14)
where

$$f(E_n)=\frac{1}{e^{\beta (E_n\mu )}+1},$$
(16)
is the Fermi distribution function and $`N_c`$ is the expectation value for the number of carriers.

Equations (11) and (12) imply that we can use Monte Carlo techniques to evaluate various thermodynamic quantities, except that we must use the carrier free energy rather than the internal energy in the Metropolis algorithm. An important point to note is that since the calculation is in the grand canonical ensemble, both the temperature and chemical potential must be kept fixed during the simulation. This differs from a recent Monte Carlo study in which $`N_c`$ was held constant at each Monte Carlo step by varying $`\mu `$ rather than holding it constant, in our simulations we hold $`\beta `$ and $`\mu `$ fixed for each run and then we average over disorder using samples with equal $`N_c`$, since we wish to average samples which have the same $`x`$ and $`p`$.

### C Magnetic and Thermodynamic Quantities

Our Monte Carlo simulations allow us to compute various magnetic and thermodynamic quantities after equilibration. We perform equilibrium thermal averages (indicated by $`\mathrm{}`$) for each sample (realization of disordered Hamiltonian) and then average over many different realizations of disorder. (The disorder average is indicated by an overbar $`\overline{\mathrm{}}`$). We collect data for both global (bulk) and local quantities. We first consider the global quantities we study. We compute the moments of the average magnetization per Mn spin, $`M^q`$, and the average magnetization per carrier, $`m^q`$, at each Monte Carlo step, where the $`q^{\mathrm{th}}`$ moments are given by

$`M^q`$ $`=`$ $`\left({\displaystyle \frac{\sqrt{|_i𝐒_i|^2}}{SN_d}}\right)^q,`$ (17)
$`m^q`$ $`=`$ $`\left({\displaystyle \frac{1}{2}}{\displaystyle \frac{\sqrt{\left|_{i,\alpha \beta }c_{i\alpha }^{}𝝈_{\alpha \beta }c_{i\beta }\right|^2}}{N_h}}\right)^q.`$ (18)
In the above equation $`S=\frac{5}{2}`$. The average spin per Mn, $`S_{Mn}`$, and the average fermion spin, $`s_c`$ are thus

$$S_{Mn}=\overline{M},s_c=\overline{m}.$$
(19)
Note that with our definition, $`S_{Mn}`$ is normalized, i.e. for a fully polarized spin state $`S_{Mn}=1`$. The negative sign for the fermion magnetization is due to the antiferromagnetic interactions between Mn spins and carriers which leads to oppositely oriented polarizations. The Mn susceptibility is given by

$$\chi _{Mn}=\beta \overline{\left[M^2M^2\right]},$$
(20)
while the carrier susceptibility is

$$\chi _h=\beta \overline{\left[m^2m^2\right]}.$$
(21)
For each measured quantity, we calculate the statistical errors in the conventional manner from the variance. One of the limitations of studying a magnetic model on small lattices, as we are forced to do here, is that it is not easy to identify the position of a thermodynamic transition from considering quantities such as the magnetization. However, from finite size scaling theory, quantities that are dimensionless are particularly useful to identify the transition temperature. One such quantity, known as the Binder cumulant, is given by

$$G(L,T)=\frac{1}{2}\left(53\overline{\left(\frac{M^4}{M^2^2}\right)}\right).$$
(22)
$`G(L,T)`$, which measures the ratio of the fourth moment to the square of the second moment of the magnetization for a finite system, is defined such that in the paramagnetic phase $`G(L,T)`$ decreases with $`L`$, and tends to zero as $`L\mathrm{}`$, while in the ferromagnetic phase it increases with increasing size $`L`$ and tends to unity in the thermodynamic limit. Near the transition temperature $`T_c`$, being dimensionless, it has the finite size scaling form $`G(L,T)=G[L^{\frac{1}{\nu }}(TT_c)]`$,  where $`\nu `$ is the exponent of the diverging spin-spin correlation length $`\xi (TT_c)^\nu `$. Consequently, at $`T_c`$, $`G(L,T_c)`$ is independent of $`L`$; $`T_c`$ can be identified by a simultaneous crossing of $`G(L,T)`$ vs. $`T`$ curves for different $`L`$. Because $`G`$ is dimensionless, and depends only on the ratio $`L/\xi `$, rather than both $`L`$ and $`\xi `$, this method is found to be more reliable in determining $`T_c`$ than analysis of the onset of magnetization, or peaks in the magnetic susceptibility in finite sized samples.

Local quantities we calculate are the local charge density at each Mn site:

$$\rho _i=\underset{j,\sigma }{}|\varphi _{ij}|^2c_{j\sigma }^{}c_{j\sigma },$$
(23)
and the local magnetization, which we define to be the average projection of the spin at site $`i`$ on the total magnetization $`𝐌`$

$$M_i^{\mathrm{local}}=𝐒_i𝐌.$$
(24)
We are interested in the individual distributions $`𝒫(\rho _i)`$ and $`\stackrel{~}{𝒫}(M_i^{\mathrm{local}})`$ and the joint distribution $`𝒫_J(\rho _i,M_i^{\mathrm{local}})`$ to characterize the local environment at different Mn sites.

## III Perturbative Monte Carlo (PMC)

In principle, to perform Monte Carlo simulations on a model with fermion and classical degrees of freedom, one needs to diagonalize the fermion part of the Hamiltonian after each “spin flip” (more precisely a spin rotation), i.e. a new choice of the classical variable. This leads to new eigenvalues which are used to compute the change in the carrier free energy $`_c(\{S^z,\varphi \})`$. This is computationally time consuming, and hence one would prefer a quicker, approximate method which is still reasonably accurate. One such approach is the hybrid Monte Carlo (HMC) algorithm used recently on the DMS problem,  and also on the double exchange model.  We have developed an algorithm that works in a similar spirit to HMC, that we describe below.

We select a spin, at site $`i`$ and allow it to perform a small rotation:

$`S_i^z`$ $``$ $`S_i^z+\delta S_i^z,`$ (25)
$`\varphi _i`$ $``$ $`\varphi _i+\delta \varphi _i,`$ (26)
where $`\delta S_i^z[\lambda /2,\lambda /2]`$, and $`\delta \varphi _i[\lambda \pi ,\lambda \pi ]`$, are restricted to a small region on the surface of the unit sphere ($`\lambda 1`$), in a way that leads to uniform sampling. We use perturbation theory to compute the change in the carrier eigenenergies due to the rotation of the spin at site $`i`$:

$$E_n(\{S^z,\varphi \})E_n(\{S^z,\varphi \})+\delta E_n$$
(27)
where

$$\delta E_n=\frac{1}{2}\delta 𝐒_i\underset{j}{}J_{ij}\underset{\alpha \beta }{}\psi _{n\alpha }^{}(j)𝝈_{\alpha \beta }\psi _{n\beta }(j).$$
(28)
As a result, $`\delta E_nO(\lambda )`$ and therefore perturbation theory is essentially exact as $`\lambda 0`$.

We now use Eq. (12) to compute the change in carrier free energy associated with the spin rotation, and the Metropolis criterion to decide whether to accept the spin rotation. We perform such an update for each spin in the system using the perturbation scheme. After a complete sweep through the system, we compute the eigenenergies $`E_n`$ and eigenfunctions $`\psi _n(i)`$ corresponding to the new spin configuration using exact diagonalization, and start a new Monte Carlo sweep. We found that this approach was quicker than if we diagonalized after every spin flip, generally by a factor of 3-4.

### A The chemical potential

One important issue is the choice of the chemical potential $`\mu `$ for the desired average number of charge carriers, $`N_c=N_h`$. To determine the chemical potential we consider two replicas of the system, one starting from a fully polarized (ferromagnetic) configuration, the other a purely random (paramagnetic) configuration. After every few Monte Carlo steps (in practice 5 MC steps worked well), we use the condition

$$N_h=\underset{n}{}\frac{1}{e^{\beta (E_n\mu )}+1},$$
(29)
to update the value of $`\mu `$ for each replica. When the magnetization $`M`$ and chemical potential of each replica agree to within 2 %, we continue to calculate the chemical potential for each replica after every 5 Monte Carlo steps, but their dynamics are determined by the mean chemical potential, $`\overline{\mu }`$. This average chemical potential is free to vary up until some equilibration time, and then the chemical potential is fixed as $`\mu =\overline{\mu }`$ where the average is over the time during the equilibration period for which $`\overline{\mu }`$ is used to calculate the carrier free energy. The equilibration time used depended on the temperature, $`x`$ and $`p`$, but was generally between 20000 and 40000 Monte Carlo steps. (The magnetization generally equilibrated within 2000 Monte Carlo steps, while the remainder of the equilibration was required to obtain an accurate value of the chemical potential). The fixed chemical potential is used for the remainder of the run, during which data is collected. We typically use 20000 to 40000 Monte Carlo steps to collect data. This procedure was found to obtain a chemical potential that yielded values of $`N_c`$ for the sample that were generally within 2% of $`N_h`$ (not surprisingly it was found to be more effective in larger samples where the relative sizes of the fluctuations are correspondingly smaller).

In several other studies of models with fermions coupled to spins  the electron occupation numbers have been taken to be those corresponding to $`T=0`$ (i.e. the Fermi distribution is replaced by a Heaviside distribution). In our simulations, we allow the filling to change as a function of temperature (i.e. we do not assume degenerate electrons). We do, however truncate the number of states we include – their number varies as a function of temperature, such that states of high energy which have $`10^5`$ or lower filling are discarded (the cutoff is sufficiently small that the results are not sensitive to its value).

### B Testing the PMC algorithm

We tested our perturbative Monte Carlo algorithm on a simple model of fermions coupled to classical spins which admits an exact solution. The Hamiltonian, shown below in Eq. (30), corresponds to fermions hopping along a 1-dimensional chain with $`N`$ sites and periodic boundary conditions. At every lattice site $`i`$, there is a classical spin $`𝐒_i`$, and the fermions are equally strongly coupled to all spins on the chain. The model Hamiltonian is thus

$``$ $`=`$ $`t{\displaystyle \underset{ij,\sigma }{}}c_{i\sigma }^{}c_{j\sigma }+{\displaystyle \frac{J}{N}}{\displaystyle \underset{i}{}}𝐒_i{\displaystyle \underset{j,\alpha \beta }{}}c_{j\alpha }^{}{\displaystyle \frac{1}{2}}𝝈_{\alpha \beta }c_{j\beta },`$ (30)
where energy is measured in units of $`t=1`$ and only nearest neighbour hopping is allowed. Note that the exchange scales as $`J/N`$ to ensure an extensive energy. The eigenvalues for this model can be calculated exactly as discussed in Appendix A. Exact results for any number of classical spins, $`N`$, and $`N_h`$ fermions can be calculated, particularly for the magnetization per spin, $`S`$, the Binder cumulant $`G(L)`$ and the fermion magnetization $`s_c`$. A comparison of the exact results for $`S(T)`$ and those calculated using PMC are shown in Fig. 1 for $`t=J`$, $`N=20`$, $`N_h=3`$ and using 40000 Monte Carlo steps.

In Fig. 2 we plot the relative error in the fermion magnetization against the number of Monte Carlo steps $`(t)`$ for the same parameters used for Fig. 1 at $`T/J=0.0064`$ where the magnetization is about $`0.5`$. At long times the convergence is of the form $`t^{\frac{1}{2}}`$, where $`t`$ is the number of Monte Carlo steps, as expected for Monte Carlo simulations. We have also made a comparison between the PMC method and diagonalization after every spin rotation for the model of DMS that we are really interested in \[Eq. (2)\], which gave agreement to within the 2 % statistical error bars. There is also some contribution to the systematic error from our approach, but we choose $`\lambda `$ to be small (generally $`\lambda 0.03`$), so that this error has little effect on our results.

## IV Results: Global Quantities

Four different values of $`x`$ and $`p`$ were chosen for the Monte Carlo study, representative of the concentrations and compensation seen in the experimental materials. These were: $`x=0.01`$, $`p=0.1`$; $`x=0.01`$, $`p=0.3`$; $`x=0.03`$, $`p=0.1`$ and $`x=0.03`$, $`p=0.3`$. Note that using Mott’s criterion as in Ref. , the $`x=0.01`$, $`p=0.1`$ sample is in the insulating phase, whilst the remainder of the cases considered are somewhere on the metallic side of the MIT. In the case of the most metallic sample $`x=0.03`$, $`p=0.3`$ it may be inappropriate to neglect the host band states for a quantitative description of the magnetic properties, but we have included this also within the context of an impurity band model for comparison with the other cases.

We carried out simulations on lattices of linear size between $`L=7`$ and $`L=14`$, which contained between $`40`$ and $`110`$ Mn spins and between $`4`$ and $`21`$ carriers. The sizes that were considered and their labels are tabulated in Table I. We averaged up to 700 samples per data point, depending on the size and temperature. Typically at least 30-40 samples were averaged for each data point. In this section we present our results for global quantities such as the magnetization, the Binder cumulant (to determine $`T_c`$) and the magnetic susceptibility. We also present our results comparing ordered and disordered samples. We focus on local quantities in the next section.

### A Magnetization

In Fig. 3 we plot the average magnetization per Mn spin, $`S_{\mathrm{Mn}}(T)`$, as a function of temperature for all four combinations of $`x`$ and $`p`$ using samples containing between 60 and 70 Mn spins in each case (thus finite size effects are similar in all four curves). Clearly, the critical temperature, $`T_c`$ increases with increasing $`x`$ and $`p`$. The magnetization curves for cases B2 and C2, which have the same hole concentration $`px`$ have similar numerical values, although the curve for $`x=0.01`$, $`p=0.3`$ appears to have a lower $`T_c`$. The most important feature is the unusual shape of the magnetization curves. The magnetization decreases rapidly from full polarization at $`T=0`$, leading to linear or concave upward shapes, similar to those found in the mean field approximation  and for insulating II-VI DMS,  and very unlike the strongly convex upward magnetization curves seen in conventional ferromagnets such as iron.

The carrier magnetizations, $`s_c(T)`$ for the same sizes shown in Fig. 3 are shown in Fig. 4. The curves mirror the feature of the Mn magnetizations that increasing $`x`$ and $`p`$ lead to greater polarization at higher temperatures – there is also a clearer distinction between curves for B2 and C2, where $`x=0.01`$, $`p=0.3`$ polarizes at noticeably lower temperatures than $`x=0.03`$, $`p=0.1`$. One major difference from the results of the mean field study is that the curves for the carrier magnetizations appear to be much more like the Mn magnetization curves, whereas in the mean field study  the carrier curves remained almost fully polarized until $`T`$ was close to $`T_c`$. This feature is in closer agreement with experiments than the mean field results.

### B Critical Temperature

We use the Binder cumulant curves, $`G(L,T)`$ (see Eq. (22)) to identify the critical temperature $`T_c`$ for each carrier and Mn concentration. As described in Sec. II C, $`T_c`$ is indicated by the simultaneous intersection of $`G(L,T)`$ curves for sufficiently large sizes, $`L`$. Fig. 5 shows $`G(L,T)`$ curve for $`x=0.03`$ and $`p=0.3`$ for both sample sizes (D1 and D2 of Table I) over a wide temperature range ($`T/J`$ = 0.02 to 0.6). As expected $`G(L,T)`$ decreases with increasing size $`L`$ at high $`T`$ whereas the reverse is true at low $`T`$. At the transition temperature $`T_c`$, $`G(L,T)`$ is expected to be independent of $`L`$, which is indicated by the crossing point of the two curves, implying $`T_c/J=0.45`$. The solid curves are spline fits to the data appropriately weighted by the error bars.

Figs. 6 \- 8 show the $`G(L,T)`$ data for the other three concentrations studied, each in the vicinity of the crossing point. Several hundred samples were generally averaged for each data point as indicated in the figure captions. The solid curves are spline fits to the data over the entire temperature range studied (typically $`T/T_c`$ = 0 to 3). For the one case where we studied more than two system sizes ($`x=0.01`$, $`p=0.1`$), which we exhibit in Fig. 6, all curves are consistent with a single intersection point, although a small size dependence cannot be ruled out.  Such a dependence would represent corerections to finite size scaling, arising from the relatively small sizes of the samples studied. (The effective linear size of the spin system is $`N_d^{\frac{1}{3}}`$, which varies from 3.4 to 5). The small sizes, necessitated by the need to repeatedly diagonalize the fermion Hamiltonian, also limited the dynamic range available in our study. Because of this, the curves for the sizes studied do not splay out very dramatically around $`T_c`$. This, in turn, necessitates calculation of $`G(L,T)`$ to high precision, which requires long runs and averaging over many samples. Despite these drawbacks, we find $`G(L,T)`$ to be a more reliable estimator of $`T_c`$ than e.g. peaks in $`\chi (T)`$ for the largest size studied. Based on Figures 5 to 8, we estimate $`T_c`$ as shown in Table II.

In Table II we compare the Curie temperature $`T_c`$ as determined in the Monte Carlo simulations here with the results of the mean-field approximation using Langevin functions for classical vector spins.

As expected, the mean-field approximation overestimates $`T_c`$. While the reduction due to fluctuations is only about $`20\%`$ for the largest $`x`$ and $`p`$ studied, it can be much more significant (a factor of 5 or more) at densities at and below the metal insulator transition density. (We mention in passing that results in the mean-field approximation for quantum spins using Brillouin functions are higher than for classical spins by $`6080\%`$).

Since there are many models for DMS which are solved at the mean field level, having an understanding of the behaviour of these models when fluctuations are considered is very important if one wants to make quantitative fits to data. For a model where Mn ions interact with carriers in an unperturbed host band, a similar reduction in $`T_c`$ was found in Monte Carlo simulations.  One important difference between the results here and those found in Ref.  however, is that here mean field theory becomes more accurate with increasing carrier concentration at a given $`x`$, whereas the opposite appears to be the case in that model.

The Monte Carlo simulations show clearly a strong decrease of $`T_c`$ at low carrier density $`(px)`$, in agreement with experiment, and rectify the unphysically large $`T_c`$ obtained in the position dependent mean-field treatment in this limit. From Table II it appears clear that the main dependence of $`T_c`$ comes from the carrier density $`(px)`$, with a relatively weaker dependence on $`(x/p)`$.

### C Magnetic Susceptibility

The Mn and carrier susceptibilities \[Eqs. (20) and (21)\] are shown in Fig. 9 and Fig. 10 respectively, as a function of temperature, for the same sizes and dopings that were used for Mn and carrier magnetizations in Fig. 3 and Fig. 4.

At high temperatures, well above $`T_c`$, one expects

$$\frac{1}{\chi }\frac{3k_BT}{n(g\mu _B)^2S^2}\left(1\frac{\theta }{T}\right),$$
(31)
for classical spins (the factor of $`S^2`$ is replaced by $`S(S+1)`$ for quantum spins). We have investigated the behaviour of the inverse susceptibility in the temperature range $`T=2T_c`$ to $`6T_c`$, and found that the effective spin $`S_{\mathrm{eff}}1.61.7S`$, which suggests spin clusters (polarons) exist in the system well above $`T_c`$. Fitting data close to $`T_c`$ ($`T_c`$ to 2 $`T_c`$) Eq. (31) gives $`S_{\mathrm{eff}}2S`$, implying a significant enhancement of the Curie constant over the high $`T`$ value. An even larger enhancement has been seen in (Ga,Mn)N.

### D Effects of Disorder

One of the major results of the mean field study was that samples with maximal disorder in the position of Mn spins were found to have a higher $`T_c`$ than those with less disorder.  This was explained as being due to carriers being able to lower their total energy more in regions with higher Mn density, by polarizing Mn spins and then hopping between these sites. However, it is expected that the mean field factorization, which assumes that the carrier spin is either directed parallel or antiparallel to the overall magnetization, tends to align “islands” that might not otherwise be aligned until lower temperatures. This would be more likely to lead to a larger decrease in $`T_c`$ for disordered samples when thermal fluctuations are considered, since global magnetization around $`T_c`$ is more likely to be destroyed than in ordered samples.

We have tested whether the finding of disorder enhancing $`T_c`$ is robust, by investigating two different cases. The first is for $`x=0.01`$, $`p=0.1`$, where we compared the weakly disordered case and fully disordered case, and the second is for $`x=0.03`$, $`p=0.3`$, where we compared the fully ordered case and the fully disordered case.

In both cases we observed that the more ordered sample had a lower $`T_c`$. However, the enhancement is much smaller than obtained within the mean field treatment, as the arguments given above would suggest.

For $`x=0.01`$, $`p=0.1`$, we determined the value of $`T_c/J0.037`$ for the fully disordered case. We define mild disorder to correspond to the situation where Mn spins are chosen to be displaced from a fully ordered Mn lattice to one of the 12 nearest neighbour sites on the fcc sublattice. We considered $`N_d=64`$, $`N_h=7`$ and $`N_d=125`$, $`N_h=13`$, and then used the Binder cumulant to estimate $`T_c`$.

Note that the number of Mn spins in our simulations are slightly different from those used for other cases, since an ordered lattice of Mn needs to be commensurate with the underlying lattice. Whilst we do not have a particularly accurate determination of $`T_c`$, it is clear in Fig. 11 that the $`G(L)`$ curves cross at a temperature well below the $`T_c`$ of the disordered case ($`T_c/J=0.037`$).

In Fig. 12 we compare the magnetization calculated for $`x=0.03`$, $`p=0.3`$ with the completely random case D2 and the completely ordered case with $`N_d=64`$, $`N_h=19`$. The $`x`$ and $`p`$ value were chosen since they give the largest $`T_c`$ and it was hoped that the effects would be more noticeable. It is clear that even when thermal fluctuations are included, the more ordered sample has a lower value of $`T_c`$. (The magnetization at high $`T`$ is a finite size effect).

## V Results: Local quantities

We now discuss our results for various local quantities such as the local magnetization and local charge densities for various Mn and carrier concentrations.

### A Charge densities

In Fig. 13 we show a histogram of the local carrier charge densities $`\rho _i`$ \[see Eq. (23)\] at all Mn sites $`i`$, for $`x=0.01`$, $`p=0.1`$ (case A3) and averaged over 89 samples at the temperature $`T/J=0.01`$, compared with the corresponding distribution for $`x=0.0093`$ and $`p=0.1`$ obtained using the mean field approximation with Brillouin functions close to $`T=0`$.  Clearly the two distributions agree very well overall. The two peaks of the histogram suggest two populations of Mn sites. The peak at $`\rho _i0.7`$ is from sites which have a high probability of trapping a carrier, whilst the broader peak at much lower values of $`\rho _i`$ is due to sites which have very little probability of having a carrier on them. For these sites, $`\rho _i`$ comes mainly through tailing from nearby sites that have a high charge carrier density \[see Eq. (23)\]. This shows that the charge carriers reside primarily on sites where there is a higher than average local Mn concentration, due to the strong interactions with these Mn spins. Because of the inhomogeneous charge distribution, spins of Mn atoms on sites devoid of charge carriers have very small antiferromagnetic couplings to the charge carriers, and therefore remain essentially free down to low temperatures. This explains the unusual shape of the magnetization curves.

If we compare the histograms for the charge densities with the data obtained in the mean field model, we notice two important features. Firstly, the histogram is much less temperature dependent for $`TT_c`$. This is probably because we are generally looking at lower temperature scales than in the mean field case. Secondly, these distributions have large widths typical of highly disordered systems. This large width was first found in the mean-field case, and is the motivation for a simplified phenomenological model for DMS, based on a two-component picture.  In this simplified model, one can divide the Mn spins into strongly and weakly interacting components depending on the temperature. This simplified model has been shown to be adequate to reproduce the results of the full distribution at the mean field level, and can also explain experimental results on a qualitative basis.

In Fig. 14 we show a histogram of the total charge densities $`\rho _i`$ at all Mn sites $`i`$, for $`x=0.01`$, $`p=0.3`$, $`N_d=69`$, $`N_h=21`$ and averaged over 41 samples at temperature $`T/J=0.1`$. The most noticeable change as the temperature is varied is that the height of the peak at large $`\rho _i`$ decreases, corresponding to carriers becoming less localized at higher temperatures.

In Fig. 15 and 16 we show the histogram of the total electron charge densities $`\rho _i`$ at each Mn site $`i`$, for $`x=0.03`$, $`p=0.1`$, $`N_d=41`$, $`N_h=4`$ at a temperature $`T/J=0.14`$, averaged over 315 samples, and for $`x=0.03`$, $`p=0.3`$, $`N_d=41`$, $`N_h=12`$ at a temperature $`T/J=0.4`$, averaged over 76 samples. (In both cases the temperatures are around 90% of $`T_c`$). Unlike the case of lower $`x`$, these distributions have only one peak, however similarly to the lower $`x`$ case, the distribution is virtually independent of $`T`$ below $`T_c`$.

Besides the change in the shape of the distribution, the width of the distribution is also significantly smaller at the higher hole concentration. This is consistent with the observation that the eigenstates at the Fermi energy are found to be delocalized at $`x=0.03`$, whilst they appear to be localized at $`x=0.01`$.

### B Local magnetizations

We now consider the local magnetizations. In Fig. 17 we show a histogram of the local magnetizations $`M_i^{\mathrm{local}}`$ \[see Eq. (24)\] at all Mn sites $`i`$, for $`x=0.01`$, $`p=0.1`$, $`N_d=110`$, $`N_h=11`$ at three temperatures, one well below $`T_c`$, one around $`T_c/2`$ and one around $`3T_c/2`$. The number of samples to generate the histogram was 10 at $`T/J=0.0053`$, 188 at $`T/J=0.02`$ and 58 at $`T/J=0.06`$. At the lowest temperature, there is a peak at around $`M_i^{\mathrm{local}}=0.7`$, and then a very broad tail that stretches to local magnetizations that are antiparallel to the overall magnetization. For intermediate temperatures ($`T=0.02JT_c/2`$), there is evidence of two populations of Mn spins illustrated by the two peaks in the histogram, whilst for high temperatures (above $`T_c`$), there is a peak centered very close to zero local magnetization with some weight for small local magnetizations. This is consistent with the two component picture described in the section on the local charge distribution – at temperatures well below $`T_c`$ most of the Mn spins are strongly coupled, leading to a peak and then a tail of more weakly coupled spins. At intermediate temperatures, there are two comparable populations of weakly and strongly coupled Mn spins, whilst at high temperatures, the local magnetizations are all small and there is no long-range magnetic order.

Fig. 18 shows a similar histogram of the local magnetizations $`M_i^{\mathrm{local}}`$, for $`x=0.01`$, $`p=0.3`$, $`N_d=69`$, $`N_h=21`$ at temperatures $`T/J=0.01`$, 0.033, 0.06, and 0.08 averaged over 12, 20, 70 and 106 samples respectively. The feature of two populations of Mn spins also appears to be present in this case, for $`T/J`$ between $`0.03`$ and $`0.06`$ ($`T_c/J=0.12`$).

Fig. 19 and 20 show the temperature evolution of the corresponding histograms for $`x=0.03`$, $`p=0.1`$, ($`N_d=41`$, $`N_h=4`$) and $`x=0.03`$, $`p=0.3`$ ($`N_d=41`$, $`N_h=12`$) for temperatures ranging from well below to just below $`T_c`$. For the former, with $`T_c/J=0.18`$ the data are for $`T/J=0.01`$, 0.033, 0.08, and 0.14, averaged over 49, 48, 48 and 334 samples respectively, whilst for the latter ($`T_c/J=0.45`$), we show data for $`T/J=0.04`$, 0.08, 0.19 and 0.40 averaged over 16, 22, 26 and 75 samples respectively.

At each temperature there is typically a peak with some breadth, but unlike the two cases with $`x=0.01`$, there is no double-peak structure seen. This can be understood as increasing $`x`$ leading to smaller relative density variations for the Mn spins, and hence smaller fluctuations in the range of local environments, resulting in narrower distributions for local charge and magnetization.

### C Correlation of charge and magnetization

Large local magnetization of the local Mn spins are correlated with large local charge densities of the carriers at the Mn site. To illustrate this point, in Fig. 21 we show the correlation of large local magnetizations and charge densities in a single sample with 110 Mn spins and 11 carriers at a temperature $`T=0.01J0.25T_c`$. Open circles correspond to sites with low charge density $`\rho _i<0.1`$ and low magnetization $`M_i^{\mathrm{local}}<0.4`$. Solid filled circles correspond to sites with $`\rho _i>0.1`$ and $`M_i^{\mathrm{local}}>0.4`$ Sites with $`\rho _i<0.1`$ but $`M_i^{\mathrm{local}}>0.4`$ or with $`\rho _i>0.1`$ and $`M_i^{\mathrm{local}}<0.4`$ are shown as half-filled circles. Less than 14% of the sites are half-filled indicating a strong correlation between $`\rho _i`$ and $`M_i^{\mathrm{local}}`$.

Fig. 22 shows the joint distribution function for local magnetization and local charge density for $`x=0.01`$, $`p=0.1`$. The data inhabits a narrow band with the peak corresponding to localized carriers also corresponding to Mn spins with large local magnetizations.

In Fig. 23 we plot the joint distribution function for local magnetization and local charge density for $`x=0.03`$, $`p=0.1`$. The distribution is similar to the previous case in that there is also a strong correlation between charge density and local magnetization. However, the data occupy a smaller region of the $`M_i^{\mathrm{local}}`$, $`\rho _i`$ plane (which is not surprising since both distributions are narrower than in the insulating case).

## VI Discussion and Conclusions

In this study we have performed a Monte Carlo study of an impurity band model for III-V DMS. In order to do this we are simulating a model of fermions coupled to classical degrees (spins) and this requires diagonalization of the fermion problem for every classical spin configuration. This is very time-consuming computationally. To speed up the procedure, we introduce quantum mechanical perturbation theory coupled with Monte Carlo, which we call the PMC scheme. The method and tests are described in detail in Sec. III. Our model restricts charge carriers to an impurity band formed from the isolated acceptor impurity states introduced by the Mn ions, interacting via antiferromagnetic exchange with Mn spins, which we treat classically. This model is based on the picture for the low doping limit of localized carriers. Experimental evidence suggests that such an impurity band exists in the vicinity of the metal-insulator transition (MIT).  The advantage of this type of model in comparison to other models that have been suggested for III-V DMS which start from a valence band point of view for the carriers, is that it naturally incorporates the physics associated with the MIT, which may be important even in the metallic region. We have considered parameters appropriate for (Ga,Mn)As, although qualitative features of the results may well apply for other III-V DMS.

The impurity band model has previously been studied using a mean-field approximation,  in which it was found that the magnetic properties are very unlike those of a conventional ferromagnet. The magnetization curves $`S_{Mn}(T)`$ were found to be concave upwards for a significant portion of the temperatures below $`T_c`$, unlike conventional convex upward curves, and it was predicted that the magnetization should be inhomogeneous. The role of disorder in the material was also examined by comparing the magnetization curves for ordered arrays of dopants and impurities placed with varying degrees of randomness. It was found that randomly placed impurities led to a higher $`T_c`$ than that found for an ordered lattice of impurities. However, in the mean field solution, the electron spin is chosen to be either parallel or antiparallel to the total magnetization. Further, mean field approximations neglect temporal fluctuations, and are known to overestimate $`T_c`$. These overestimations could be significantly different in the ordered and disordered cases.

The results of the Monte Carlo simulations confirm several of the results of the mean field study, while differing on some others. Firstly, the shapes of the magnetization versus temperature curves, $`S_{Mn}(T)`$ are found to be unusual, as in the mean field study, but more linear in $`T`$ compared with the strongly concave upward curves in the mean field approximation. Further, the carrier and Mn spin magnetization follow each other more closely, unlike the mean field, where the carriers remain polarized to much higher temperatures than the Mn spins. We believe this difference is a direct consequence of the neglect of temporal fluctuations in the mean field approximation. The Monte Carlo data appear to be much more in accord with experimental results, which suggest that both Mn and carrier spins depolarize with roughly the same temperature dependence.  In fact, the magnetization curves observed in this study bear a striking resemblence to those obtained with magnetic circular dichroism results in Ref. . Whilst magnetic susceptibilities were not considered in the mean field study, in a Monte Carlo study applicable to II-VI semiconductors, a peak was found in the ordered phase, well below $`T_c`$, in addition to the singularity at $`T_c`$.  Our susceptibility data show a single peak; however, the peak is at a temperature significantly below the $`T_c`$ obtained from the $`G(L,T)`$ curves. This is expected, since our sample sizes are small and the two peaks are not separated as a result. The relative proximity of the low $`T`$ peak to $`T_c`$ compared to the II-VI case, makes this more difficult to resolve. Nevertheless, the explanation of the peak being due to free or partially-free spins that are outside the percolating magnetic cluster appears to be viable in both cases.

To study the inhomogeneities in the magnetic behaviour identified in the mean field study, we have studied distributions of local quantities – charge density and magnetization. We have considered their joint distribution as well, to investigate the correlation between the two. The calculations of the local charge density at low temperatures for the lowest density ($`x=0.01`$ and $`p=0.1`$), which appears to be insulating, are in very close accord with those obtained previously at the mean field level.  The generic feature that appears to be present is that the local charge density has significant dispersion. For $`x=0.01`$ we find a two-peaked structure in the local charge density distribution, corresponding to some sites having quasi-localized carriers and others having very low charge density, whilst for $`x=0.03`$ there is still a broad distribution of charge densities, but there appears to be only one peak, rather than two, indicating that the carriers are in general delocalized. The distributions of local magnetizations appear to reflect the same physical picture – for $`x=0.01`$ at temperatures approximately $`T_c/2`$, there are two peaks, which appears to correspond to two populations of Mn spins – one which is still strongly magnetized, and one which is very slightly magnetized, whereas for $`x=0.03`$ this double-peaked structure is not observed. The $`x=0.01`$ behaviour was previously predicted from mean field calculations,  and formed the basis for a phenomenological model that has been shown to have the capacity to describe experimental magnetization curves.  Finally, the joint distribution of local charge density and local magnetization indicate that there is a strong correlation between higher charge density and large local magnetization.

One of the main aims of the Monte Carlo study was to discover how important fluctuations are in determining the critical temperature $`T_c`$ at various $`x`$ and $`p`$. As shown in Table II mean field theory appears to be more precise in this model for increasing hole concentration, which is in contrast with Monte Carlo simulations on a valence band model where the opposite trend was observed.

A second item of interest is the quantification of the dependence of $`T_c`$ on the Mn concentration $`x`$ and the carrier density $`px`$. For individual Mn coupled to free carriers, several models exist involving a combination of Fermi energy $`E_F(px)^{\frac{2}{3}}`$, and exchange energy $`E_{ex}Jx`$. The generic dependence of $`T_c`$ is found to be of the form $`T_cx^\alpha (px)^\beta `$, where mean field estimates  give $`\alpha =1`$, with $`\beta `$ varying between $`\frac{1}{3}`$ and 1, whilst an analysis involving collective spin wave excitations with RKKY interactions yields $`\alpha =2`$, $`\beta =\frac{1}{3}`$ for weak coupling ($`E_{ex}E_F`$) and $`\alpha =\frac{1}{3}`$, $`\beta =1`$ for strong coupling ($`E_{ex}E_F`$). That such a dependence exists for an impurity band is not clear; nevertheless in the range of Mn and carrier concentrations studied ($`x=0.01`$ \- 0.03, $`px=`$ 0.001 - 0.009), our $`T_c`$ can be fit by an expression of the form $`T_cx^\alpha (px)^\beta `$ with $`\alpha =0.5\pm 0.15`$, $`\beta =0.85\pm 0.15`$. As can be seen in Fig. 24, a double logarithmic plot of $`T_c/J`$ versus $`x(px)^{\frac{5}{3}}`$ yields a straight line with a slope of $`\frac{1}{2}`$. Since the dynamic range in carrier concentration is larger, the data restrict the range of allowed $`\beta `$ more than $`\alpha `$. It is interesting that both the dependence on carrier density for fixed Mn concentration ($`x`$), $`T_cp^\beta `$ and on Mn concentration for fixed degree of compensation ($`p`$), $`T_cx^{\alpha +\beta }`$, yield exponents $`\beta =\frac{5}{6}`$, $`\alpha +\beta =\frac{4}{3}`$, which lie within the range predicted by various treatments for free carriers, for which $`\beta `$ varies from $`\frac{1}{3}`$ to 1 and $`\alpha +\beta `$ from $`\frac{2}{3}`$ to 2. Determination of the separate dependences on carrier and Mn concentration in experiment would be useful in clarifying the applicability of various theoretical models.

Experimental magnetization measurements also suggest that there is some degree of variation of local environment for Mn spins, as suggested here. When SQUID magnetometer measurements and anomalous Hall effect measurements are compared, the magnetization curves that are observed have very different characteristic behaviour. The SQUID measurements tend to be non-conventional magnetization curves as observed here, since they sample all Mn spins. However, the Hall effect measurements are made in transport and give curves that look quite conventional. It has been suggested  that this can be understood by the fact that the carriers that are measured in transport preferentially interact with the Mn spins that are strongly coupled to carriers, so that the magnetic properties of only one population of Mn spins are sampled in transport. The correlation between local charge density and magnetization seen here supports this scenario.

There are still many experimental questions that remain unresolved about the nature of the magnetic behaviour in (Ga,Mn)As. A more systematic understanding of the dependence of $`T_c`$ on $`x`$ and $`p`$ will help determine which models are more appropriate for which regions of the phase diagram. In particular, more accurate determinations of $`p`$ appear to be one of the most important ingredients. Recent steps in this direction where the carrier concentration has been controlled with electron doping are a start,  but this may introduce other complications due to dopant centers being spatially distinct from the spins, unlike the case for Mn, where the two are in the same place. Local probes of the material like nuclear magnetic resonance will also help to uncover to what extent the magnetic environment is inhomogeneous, and how this depends on $`p`$ and $`x`$.

Despite the qualitative agreements with experimental data, there are a number of effects that are left out in the treatment of DMS by our model. Within the mean field approximation, a number of hopping integrals have been considered within a tight-binding description of the impurity band, and it appears that the most important feature in determining $`T_c`$ is the density of states at the Fermi energy.  Other effects that have been left out in this model are carrier-carrier interactions, valence band states, spin-orbit effects and direct Mn-Mn interactions which could lead to frustration.  Whilst the direct Mn-Mn exchange should not be important for the $`x`$ values considered here, since there are relatively few Mn spins that are close enough for their antiferromagnetic exchange to be important, they may become significantly more important at higher $`x0.1`$, and affect $`T_c`$.

In conclusion, we have performed Monte Carlo simulations on an impurity band model for III-V DMS. We have confirmed many of the features that were seen at the level of mean field simulations – unusual magnetization curves and inhomogeneous magnetization and charge density below $`T_c`$. The unusual magnetization curves are also in qualitative agreement with many experiments. We have made a comparison of $`T_c`$ determined using each method and find that for larger values of Mn concentration and carrier concentration, the mean field determination of $`T_c`$ become more accurate. We also find a power-law relation between $`T_c`$ and carrier concentration and Mn concentration, that could be compared with experiments.

## VII ACKNOWLEDGEMENTS

We thank Xin Wan for helpful suggestions both in the early stages of this work and during preparation of the manuscript. This research was supported by NSF grant DMR-9809483. M. B. acknowledges support from the Natural Sciences and Engineering Research Council of Canada.

## A Exact solution for the toy model

Consider the Hamiltonian of Eq. (30)

$$=t\underset{i,j}{}c_{i\sigma }^{}c_{j\sigma }+\frac{J}{N}𝐒\underset{j=1}{\overset{N}{}}c_{j\alpha }^{}\frac{1}{2}𝝈_{\alpha \beta }c_{j\beta },$$
(A1)
where

$$𝐒=\underset{i=1}{\overset{N}{}}𝐒_i,$$

is the total spin of the system.

We parametrize $`𝐒=S(\mathrm{sin}\theta \mathrm{cos}\varphi ,\mathrm{sin}\theta \mathrm{sin}\varphi ,\mathrm{cos}\theta )`$ and perform the canonical transformation

$`d_i`$ $`=`$ $`\mathrm{cos}{\displaystyle \frac{\theta }{2}}c_i+\mathrm{sin}{\displaystyle \frac{\theta }{2}}e^{i\varphi }c_i,`$
$`d_i`$ $`=`$ $`\mathrm{sin}{\displaystyle \frac{\theta }{2}}e^{i\varphi }c_i+\mathrm{cos}{\displaystyle \frac{\theta }{2}}c_i.`$
It is then straightforward to show that

$$=t\underset{i,j}{}d_{i\sigma }^{}d_{i\sigma }+\frac{JS}{2N}\underset{i=1}{\overset{N}{}}\left(d_i^{}d_id_i^{}d_i\right).$$

As a result, we can diagonalize the Hamiltonian in the $`\stackrel{}{k}`$-space

$$=\underset{k}{}E_{k\sigma }(S)d_{k\sigma }^{}d_{k\sigma },$$
(A2)
where

$$E_{k\sigma }(S)=2t\mathrm{cos}ka+\frac{1}{2N}J_oS\sigma ,$$
(A3)
with $`\sigma =\pm 1`$ and the cyclic boundary conditions imply that $`ka=2\pi n/N`$, with $`n=0,1,\mathrm{},N1`$. There is a single lowest eigenvalue corresponding to $`k=0`$, and then degenerate eigenvalues corresponding to left and right moving modes.

Thus the grand-canonical partition function is \[see Eq. (8)\]

$$Z_N=\underset{i=1}{\overset{N}{}}\left[𝑑\mathrm{\Omega }_i\right]\underset{k\sigma }{}\left(1+e^{\beta (E_{k\sigma }(S)\mu )}\right),$$
(A4)
where we use the simplified notation $`𝑑\mathrm{\Omega }_i=_0^\pi \mathrm{sin}\theta _id\theta _i_0^{2\pi }𝑑\varphi _i`$ for the integral over solid angle. One can avoid the $`2N`$ multiple integrals over individual spin angles, and replace them by an expression of the general form

$$Z_N=_0^{2\pi }𝑑\varphi _0^\pi 𝑑\theta \mathrm{sin}\theta _0^N𝑑SS^2F_N(S,\theta ,\varphi )$$

$$\times \underset{k\sigma }{}(1+e^{\beta (E_{k\sigma }(S)\mu )}),$$
(A5)
where $`\theta ,\varphi `$ define the orientation of the total spin $`𝐒`$, and we use the fact that the magnitude of the total spin, $`S`$, varies between 0 and $`N`$. Comparing Eqs. (A4) and (A5), the definition of the “weight” $`F_N(𝐒)`$ is

$$F_N(𝐒)=\underset{i=1}{\overset{N}{}}\left[𝑑\mathrm{\Omega }_i\right]\delta \left(𝐒\underset{i=1}{\overset{N}{}}𝐒_i\right),$$
(A6)
from which it is straightforward to derive the recurrence relation

$$F_N(𝐒)=𝑑\mathrm{\Omega }_NF_{N1}(𝐒𝐒_N)\theta (N1|𝐒𝐒_N|),$$

where the Heaviside function $`\theta `$ insures that the argument of $`F_{N1}`$ cannot have a magnitude larger than $`N1`$. It is apparent that, in fact, $`F_N(𝐒)=F_N(S)`$. This can easily be seen from the definition (A6) as well, since one can choose to define the angles $`\mathrm{\Omega }_i`$ with respect to the system of coordinates in which $`𝐒=S𝛀=S𝐳`$, and the result cannot depend on the particular orientation $`\theta ,\varphi `$ of the total spin $`𝐒`$. Using the new variable $`y=\mathrm{cos}\theta _N`$, and performing the trivial integral over $`\varphi _N`$, the recurrence relation can be rewritten as:

$$F_N(S)=2\pi _1^1𝑑yF_{N1}(\sqrt{1+S^22Sy})$$

$$\times \theta (N1\sqrt{1+S^22Sy}).$$
(A7)
Defining

$$F_N(S)=(4\pi )^{N1}f_N(S),$$
(A8)
we obtain the recurrence formula

$$f_N(x)=\frac{1}{2x}_{|x1|}^{min(x+1,N1)}𝑑zzf_{N1}(z).$$
(A9)
This is supplemented by the “initial” condition

$$f_2(x)=\frac{1}{2x},$$
(A10)
which can be obtained through direct integration from Eq. (A6).

From the recurrence relation Eq. (A9), one can determine the general solution for the weight function $`f_N(x)`$ for any integer $`N`$. This is a piecewise function, given by

$$f_N(x)=\frac{1}{2^{N1}(N2)!x}\underset{k=0}{\overset{m}{}}(1)^kC_N^k(N2kx)^{N2}$$

on the subinterval $`N2(m+1)xN2m`$ of the support interval $`[0,N]`$, with $`m`$ an integer $`0m<N/2`$ and where $`C_N^k=N!/(k!(Nk)!)`$ is the appropriate binomial coefficient.

The partition function is thus

$$Z_N=(4\pi )^N_0^N𝑑SS^2f_N(S)\underset{k\sigma }{}\left(1+e^{\beta (E_{k\sigma }(S)\mu )}\right),$$

where the chemical potential $`\mu `$ is fixed from the condition for the average number of fermions:

$`n_c={\displaystyle _0^N}𝑑SS^2f_N(S){\displaystyle \underset{q\alpha }{}}n_{q\alpha }{\displaystyle \frac{\underset{k\sigma }{}\left(1+e^{\beta (E_{k\sigma }(S)\mu )}\right)}{Z_N}},`$ (A11)
where $`n_{q\alpha }=[e^{\beta (E_{q\alpha }(S)\mu )}+1]^1`$ are the occupation numbers of the fermionic levels. Since an analytical expression is available for the weight function $`f_N(S)`$, these integrals can be evaluated numerically. Once the chemical potential $`\mu `$ is known, any other expectation value, such as the total spin magnitude

$$S=\frac{1}{Z_N}_0^N𝑑SS^3f_N(S)\underset{k\sigma }{}\left(1+e^{\beta (E_{k\sigma }\mu )}\right)$$

or the total fermionic spin

$$|s|=_0^N𝑑SS^2f_N(S)\underset{q\alpha }{}\frac{\alpha }{2}n_{q\alpha }\frac{\underset{k\sigma }{}\left(1+e^{\beta (E_{k\sigma }\mu )}\right)}{Z_N}$$

can be computed for any given temperature $`T`$.


format: text

text:
# References

LC-TH-2002-005

Logarithmic Fingerprints of Virtual Supersymmetry Partially supported by EU contract HPRN-CT-2000-00149

M. Beccariaa,b, F.M. Renardc and C. Verzegnassid,e
aDipartimento di Fisica, Università di Lecce, Via Arnesano, 73100 Lecce, Italy.

bINFN, Sezione di Lecce

c Physique Mathématique et Théorique, UMR 5825

Université Montpellier II, F-34095 Montpellier Cedex 5.

d Dipartimento di Fisica Teorica, Università di Trieste, Strada Costiera 14, Miramare (Trieste)

e INFN, Sezione di Trieste

Abstract

We consider the high energy behaviour of the amplitudes for production of leptons, quarks, Higgs bosons, sleptons, squarks, gauge bosons, charginos and neutralinos at lepton colliders. We concentrate our discussion on the terms arising at one loop which grow logarithmically with the energy, typically $`[aln\frac{s}{m^2}ln^2\frac{s}{m^2}]`$. We show that in each of the above processes the coefficient ”a” reflects in a remarkable way the basic gauge and Higgs structure of the underlying interactions. A comparison with experiments at future colliders should thus provide a clean way to test the validity of the MSSM structure.

It is by now well-known that the electroweak radiative corrections to standard processes increase strongly with the center of mass energy $`\sqrt{s}`$. This arises already at the one loop level due to the presence of large double(DL) and single(SL) logarithms $`\frac{\alpha }{\pi }ln^2\frac{s}{m^2}`$ , $`\frac{\alpha }{\pi }ln\frac{s}{m^2}`$, . In the TeV range such terms reach the several tens of percent and should be easily measurable (and analyzable) at future lepton colliders  whose experimental accuracy should be at the few permille level.

The relevance of these large logarithmic effects at high energy colliders has been stressed recently in the process $`e^+e^{}f\overline{f}`$ for both the SM  and the Minimal Supersymmetric Standard Model (MSSM) cases , and in the process of production of scalar pairs in the MSSM . The processes $`e^+e^{}\gamma \gamma ,\gamma Z,ZZ`$  has been analyzed in the same spirit as well as $`\gamma \gamma f\overline{f}`$  measurable at photon-photon colliders. The logarithmic structure of the processes $`e^+e^{}W^+W^{},H^+H^{},H_a^0H_b^0,\chi ^+\chi ^{},\chi ^0\chi ^0`$ has also been analysed and we extract some results from a publication in preparation .

In order for these analyses of the MSSM cases to be already meaningful in the one TeV range, a necessary condition is that of a light SUSY scenario in which all the SUSY masses relevant for the considered process are supposed to be not heavier than a few hundred GeV. This condition can be rescaled for higher energies provided that $`M_{SUSY}<<\sqrt{s}`$. However, the large size of the effects in the several TeV range would require a treatment (even approximate) of the higher order effects in order to obtain a good theoretical prediction (at the one percent level). Such attempts have already started  and it has been claimed that DL as well as universal and angular dependent SL can be exponentiated, see also  for an application to the MSSM case. This means that the one loop structure is the basic ingredient on which we can concentrate our discussion.

Comparing the SM and the MSSM logarithmic effects in these various processes, we were impressed by a number of recurrent, impressively simple differences. The purpose of this short note is that of presenting in a systematic way these differences and of discussing their intuitive, deep physical origins.

The logarithmic terms at the one loop level

At the one loop level the logarithmic terms appearing in $`e^+e^{}`$ processes can be separated into three categories, Parameter Renomalization (PR) terms, Universal terms and Angular dependent terms :

(1) the PR terms

These are the well-known terms generated by the gauge boson or gaugino self-energy contributions. They are process dependent, but can be computed in a straightforward way from the Born amplitude i.e.,

$$\frac{1}{4\pi ^2}[g_1^4\beta _1\frac{dA^{Born}}{dg_1^2}+g_2^4\beta _2\frac{dA^{Born}}{dg_2^2}][ln\frac{s}{\mu ^2}]$$
(1)
where $`\beta _1=\frac{5}{9}N_{fam}\frac{1}{24},(\frac{5}{6}N_{fam}\frac{1}{4})`$ and $`\beta _2=\frac{1}{3}N_{fam}+\frac{43}{24},(\frac{1}{2}N_{fam}+\frac{5}{4})`$ are the SM, (MSSM) values of the usual RG functions associated to the gauge couplings $`g_{1,2}`$.

(2) Universal terms

Also called ”Sudakov” terms, these terms appear to be typically of the form $`[aln\frac{s}{m^2}ln^2\frac{s}{m^2}]`$. They factorize the Born amplitude in a process independent and angular independent fashion. They are specific of the quantum numbers of each external particle. In the high energy range they are the dominant terms.

The single logs correspond to collinear singularities and the double logs to a coincidence of collinear and soft singularities in the one loop diagrams.

Although the results were first obtained by canonical computations of self-energy, triangle and box diagrams, in the covariant $`\xi =1`$ gauge, or in an axial gauge, the simplest way to obtain these log terms is through the splitting function formalism . With the splitting functions

$$\frac{1+x^2}{1x},\frac{1x}{2},\frac{2x}{1x},\frac{1}{2},2[\frac{x}{1x}+\frac{1x}{x}+x(1x)],\frac{x^2+(1x)^2}{2},x(1x),$$

for $`fgf`$, $`fsf`$, $`sgs`$, $`sf\overline{f}`$, $`ggg`$, $`gf\overline{f}`$, $`gs\overline{s}`$ (where $`f`$,$`s`$,$`g`$ represent fermions, scalars and gauge bosons), respectively, and the addition of the parameter renormalization terms, one immediately recovers the results of the diagrammatic computations.

(3) Angular dependent terms

They are just residual parts of the soft-collinear singularity arising from box diagrams involving gauge boson exchanges which generate $`ln^2|x|`$ terms ($`xt\frac{s}{2}(1cos\theta )`$ or $`u\frac{s}{2}(1+cos\theta )`$, being the Mandelstam parameters). After having extracted the universal angular independent part $`ln^2s`$, one remains with an angular dependent terms of the type

$$2ln\frac{s}{m^2}ln\frac{|x|}{s}+ln^2\frac{|x|}{s}$$

whose contribution constitutes an additional process-dependent single log. There are only few such terms, which are all of pure SM origin and have been computed for all considered processes.

Among these 3 types of terms the richest structure is in the Sudakov part (2), on which we now concentrate. We write this contribution as

$$A^{1loop}=[1+\frac{\alpha }{\pi }c]A^{Born}$$
(2)
This factorization applies to each external line of the process, $`c`$ being a coefficient that depends on the nature of the external particle, on the type of interaction and on the energy.

When the external particle is one member of two mixed states ($`i=1,2`$; see the examples below), the above equation has to be written in a matrix form:

$$A_i^{1loop}=\underset{j}{}[\delta _{ij}+\frac{\alpha }{\pi }c_{ij}]A_j^{Born}$$
(3)
We now list a number of typical examples.

chiral lepton or quark ($`f_{L,R}`$)

in the SM

$`c(f_{L,R})={\displaystyle \frac{1}{8}}[3ln{\displaystyle \frac{s}{M_W^2}}ln^2{\displaystyle \frac{s}{M_W^2}}][{\displaystyle \frac{I_f(I_f+1)}{s_W^2}}+{\displaystyle \frac{Y_f^2}{4c_W^2}}]_{L,R}`$
$`+[ln{\displaystyle \frac{s}{M^2}}][{\displaystyle \frac{1}{32s_W^2M_W^2}}]\{[m_t^2+m_b^2][\delta _{f,t_L}+\delta _{f,b_L}]+2m_t^2\delta _{f,t_R}+2m_b^2\delta _{f,b_R}\}`$ (4)
in the MSSM

$`c(f_{L,R})={\displaystyle \frac{1}{8}}[2ln{\displaystyle \frac{s}{M_W^2}}ln^2{\displaystyle \frac{s}{M_W^2}}][{\displaystyle \frac{I_f(I_f+1)}{s_W^2}}+{\displaystyle \frac{Y_f^2}{4c_W^2}}]_{L,R}`$
$`+[ln{\displaystyle \frac{s}{M^2}}][{\displaystyle \frac{1}{32s_W^2M_W^2}}]\{[2m_t^2(1+\mathrm{cot}^2\beta )+2m_b^2(1+\mathrm{tan}^2\beta )][\delta _{f,t_L}+\delta _{f,b_L}]`$
$`+[4m_t^2(1+\mathrm{cot}^2\beta )]\delta _{f,t_R}+[4m_b^2(1+\mathrm{tan}^2\beta )]\delta _{f,b_R}\}`$ (5)
One recognizes the Yukawa part appearing for heavy quarks, where $`\beta `$ is the mixing angle between the vacuum expectation values of the up and down Higgs chiral superfield (in standard notation $`\mathrm{tan}\beta =v_u/v_d`$). The scale $`M`$ which appears in the single logs is in principle the value of the highest mass running inside the corresponding loop. In the SM case it should be the top quark mass; in the MSSM case it will be a heavy squark or a chargino mass. We shall come back to this point in the final discussion.

transverse $`W_T^\pm ,\gamma ,Z_T`$ in the SM and in the MSSM

$$c(W)=\frac{1}{4s_W^2}[ln^2\frac{s}{M_W^2}]$$
(6)
$$c_{\gamma \gamma }=\frac{1}{4}[ln^2\frac{s}{M_W^2}]c_{ZZ}=\frac{c_W^2}{4s_W^2}[ln^2\frac{s}{M_W^2}]c_{\gamma Z}=\frac{c_W}{4s_W}[ln^2\frac{s}{M_W^2}]$$
(7)
neutral Higgs and charged or neutral Goldstones in the SM

$`c(H_{SM})=c(G^0)=c(G^\pm )`$ $`=`$ $`{\displaystyle \frac{(1+2c_W^2)}{32s_W^2c_W^2}}[4ln{\displaystyle \frac{s}{M_W^2}}ln^2{\displaystyle \frac{s}{M_W^2}}]`$ (8)
$`+{\displaystyle \frac{3}{16s_W^2M_W^2}}[m_t^2+m_b^2][ln{\displaystyle \frac{s}{M^2}}]`$
The charged and neutral Goldstone states are equivalent at high energy to the longitudinal $`W_L^\pm `$ and $`Z_L`$ components.

sleptons or squarks in the MSSM ($`\stackrel{~}{f}_{L,R}`$)

Same expression as for leptons and quarks in the MSSM.

charged and neutral Higgs bosons and Goldstones in the MSSM

A first $`2\times 2`$ matrix describes the ($`H^\pm ,G^\pm W_L^\pm `$) set, as well as the ($`A^0,G^0Z_L`$) set:

$$c_{11}=\frac{(1+2c_W^2)}{32s_W^2c_W^2}[2ln\frac{s}{M_W^2}ln^2\frac{s}{M_W^2}]+\frac{3}{16s_W^2M_W^2}[m_t^2cot^2\beta +m_b^2tan^2\beta ][ln\frac{s}{M^2}]$$
(9)
$$c_{22}=\frac{(1+2c_W^2)}{32s_W^2c_W^2}[2ln\frac{s}{M_W^2}ln^2\frac{s}{M_W^2}]+\frac{3}{16s_W^2M_W^2}[m_t^2+m_b^2][ln\frac{s}{M^2}]$$
(10)
$$c_{12}=\frac{3}{16s_W^2M_W^2}[m_t^2\mathrm{cot}\beta m_b^2\mathrm{tan}\beta ][ln\frac{s}{M^2}]$$
(11)
A second matrix describes the ($`H^0,h^0`$) set

$`c_{H^0H^0}`$ $`=`$ $`{\displaystyle \frac{(1+2c_W^2)}{32s_W^2c_W^2}}[2ln{\displaystyle \frac{s}{M_W^2}}ln^2{\displaystyle \frac{s}{M_W^2}}]`$ (12)
$`+{\displaystyle \frac{3}{16s_W^2M_W^2}}[m_t^2\mathrm{sin}^2\alpha (1+\mathrm{cot}^2\beta )+m_b^2\mathrm{cos}^2\alpha (1+\mathrm{tan}^2\beta )][ln{\displaystyle \frac{s}{M_W^2}}]`$
$`c_{h^0h^0}`$ $`=`$ $`{\displaystyle \frac{(1+2c_W^2)}{32s_W^2c_W^2}}[2ln{\displaystyle \frac{s}{M_W^2}}ln^2{\displaystyle \frac{s}{M^2}}]`$ (13)
$`+{\displaystyle \frac{3}{16s_W^2M_W^2}}[m_t^2\mathrm{cos}^2\alpha (1+\mathrm{cot}^2\beta )+m_b^2\mathrm{sin}^2\alpha (1+\mathrm{tan}^2\beta )][ln{\displaystyle \frac{s}{M^2}}]`$
$$c_{H^0h^0}=\frac{3\mathrm{cos}\alpha \mathrm{sin}\alpha }{16s_W^2M_W^2}[m_t^2(1+\mathrm{cot}^2\beta )+m_b^2(1+\mathrm{tan}^2\beta )][ln\frac{s}{M^2}]$$
(14)
where $`\alpha `$ is the mixing angle between the neutral CP even physical Higgs bosons; at tree level, $`\alpha `$ is a simple combination of $`\mathrm{tan}\beta `$ and the masses of the neutral (CP even and odd) physical Higgs bosons.

charginos $`\chi _i^+`$ in the MSSM

$`c_{\chi _i^+\chi _j^+}={\displaystyle \frac{1}{4s_W^2}}[ln^2{\displaystyle \frac{s}{M_W^2}}](Z_{1i}^+Z_{1j}^+P_L+Z_{1i}^{}Z_{1j}^{}P_R)+`$
$`{\displaystyle \frac{(1+2c_W^2)}{32s_W^2c_W^2}}[2ln{\displaystyle \frac{s}{M_W^2}}ln^2{\displaystyle \frac{s}{M_W^2}}](Z_{2i}^+Z_{2j}^+P_L+Z_{2i}^{}Z_{2j}^{}P_R)+`$
$`[ln{\displaystyle \frac{s}{M^2}}]({\displaystyle \frac{3}{16s_W^2M_W^2}})[m_t^2(1+\mathrm{cot}^2\beta )Z_{2i}^+Z_{2j}^+P_L+m_b^2(1+\mathrm{tan}^2\beta )Z_{2i}^{}Z_{2j}^{}P_R]`$ (15)
The mixing matrix elements $`Z_{1i}^\pm `$ correspond  to the charged gaugino components and $`Z_{2i}^\pm `$ to the charged higgsino components.

neutralinos $`\chi _i^0`$ in the MSSM

$`c_{\chi _i^0\chi _j^0}={\displaystyle \frac{1}{4s_W^2}}[ln^2{\displaystyle \frac{s}{M_W^2}}](Z_{2i}^NZ_{2j}^N)+`$
$`[{\displaystyle \frac{(1+2c_W^2)}{32s_W^2c_W^2}}][2ln{\displaystyle \frac{s}{M_W^2}}ln^2{\displaystyle \frac{s}{M_W^2}}](Z_{4i}^NZ_{4j}^N+Z_{3i}^NZ_{3j}^N)+`$
$`[ln{\displaystyle \frac{s}{M^2}}]({\displaystyle \frac{3}{16s_W^2M_W^2}})[m_t^2(1+\mathrm{cot}^2\beta )Z_{4i}^NZ_{4j}^N+m_b^2(1+\mathrm{tan}^2\beta )Z_{3i}^NZ_{3j}^N]`$ (16)
The mixing matrix elements $`Z_{2i}^N`$ correspond  to the neutral gaugino $`(\stackrel{~}{W_3})`$ components (there is no contribution from the Bino $`\stackrel{~}{B}`$), and $`Z_{3i}^N,Z_{4i}^N`$ to the neutral higgsino components.

Discussion of the SM terms

The $`[ln^2\frac{s}{M^2}]`$ terms, in which $`M=M_W`$ or $`M_Z`$, arise from the coincidence of soft and collinear singularities. They only appear (because of helicity conservation vertices) in SM gauge terms. They correspond to the term $`1/(1x)`$ in the splitting function with emission of a gauge boson. The minus sign is fixed by unitarity (positivity of the transition probability).

These SM gauge terms contain also a single $`[ln\frac{s}{M^2}]`$ part arising from the remaining part of the gauge splitting functions. For a fermion line one obtains the combination $`[3ln\frac{s}{M^2}ln^2\frac{s}{M^2}]`$ and for a scalar line $`[4ln\frac{s}{M^2}ln^2\frac{s}{M^2}]`$. The factor $`3`$ or $`4`$ can be traced back to the spin nature of the gauge vertices $`f\overline{f}g`$, $`s\overline{s}g`$ (where $`g`$ is a gauge boson) and more technically to the Lorentz transformation from the c.m. frame to the collinear frame of the usual $`1+\mathrm{cos}^2\theta `$ and $`\mathrm{sin}^2\theta `$ distributions of fermion or scalar pairs.

We have also obtained SL of Higgs origin due to Yukawa couplings to heavy quarks. They appear both in heavy quark production processes and in other final states involving Higgs and Goldstones (where the heavy quarks contribute virtually). The minus sign in front of the SL is also a consequence of unitarity.

Discussion of the SUSY terms and of the complete MSSM terms

Additional single logarithmic terms $`[ln\frac{s}{M^2}]`$ arise from diagrams involving scalar couplings of supersymmetric particles (sfermions, charginos, neutralinos, charged or neutral Higgses). They have also two different origins. First, a gauge origin, with the same gauge couplings as in SM terms, so that the complete MSSM combination is now $`2lnln^2`$. Secondly, a Yukawa origin, but in this case the extended Higgs structure generates new contributions depending on the parameter $`tan\beta `$, and in the case of external $`H^0,h^0`$, also on the mixing angle $`\alpha `$. In the MSSM (and except for the $`H^0,h^0`$ case) $`tan\beta `$ is the only new SUSY parameter which enters the asymptotic expressions, and leads to $`m_t^2cot^2\beta `$ and $`m_b^2tan^2\beta `$ terms. The minus sign in front of the SL is also a consequence of unitarity (positivity of the corresponding transition probability). As already mentioned, the scale $`M`$ which appears in these single logs is in principle the value of the highest mass running inside the corresponding loop, but at logarithmic accuracy, provided that one is especially interested in the slope in $`logs`$, as suggested in Ref., the choice of $`M`$ is harmless.

A list of benchmark features

We now underline the benchmark features which arise from the above results.

A first feature is that, for lepton and quark production the SM ”fermion-gauge” combination

$`[3ln\frac{s}{M_W^2}ln^2\frac{s}{M_W^2}]`$ becomes $`[2ln\frac{s}{M_W^2}ln^2\frac{s}{M_W^2}]`$

in the MSSM when gaugino terms are added.

The second feature is that, for slepton and squark production the SM ”scalar-gauge” combination

$`[4ln\frac{s}{M_W^2}ln^2\frac{s}{M_W^2}]`$ also becomes $`[2ln\frac{s}{M_W^2}ln^2\frac{s}{M_W^2}]`$

when the corresponding gaugino terms are added.
So the combination $`[2ln\frac{s}{M_W^2}ln^2\frac{s}{M_W^2}]`$ appears to be the typical MSSM combination of the whole fermion-sfermion supermultiplet.

The transformation of $`[4ln\frac{s}{M_W^2}ln^2\frac{s}{M_W^2}]`$ into $`[2ln\frac{s}{M_W^2}ln^2\frac{s}{M_W^2}]`$ also occurs when going from SM Higgs $`H_{SM}`$ or Goldstone $`G^{\pm ,0}W^\pm ,Z`$ production, to MSSM charged or neutral Higgses $`H^\pm ,H^0,h^0,A^0`$ or Goldstones.

For transverse gauge boson lines, only the quadratic term $`[ln^2\frac{s}{M_W^2}]`$ appears, both in the SM and in the MSSM. This should provide a test of the assumed ”minimal” gauge structure of the MSSM, i.e. of absence of additional (higher) gauge bosons.

In the case of chargino and neutralino production, the $`[2ln\frac{s}{M_W^2}ln^2\frac{s}{M_W^2}]`$ combination can also be found in the Higgsino components, and the $`[ln^2\frac{s}{M_W^2}]`$ term in the gaugino components. This leads to an additional potential check of the assumed supersymmetric nature of the interactions of these particles which can be achieved by a measurement of the production rate of the two charginos and of the four neutralinos.

Finally, we consider the Yukawa terms which contribute a term $`[ln\frac{s}{M^2}]`$ in the production of heavy quarks, heavy squarks, Higgs bosons and Goldstones, charginos and neutralinos. The interesting feature here is that the SM parameters $`m_t^2`$ and $`m_b^2`$ are, in the MSSM, replaced by terms that also contain the products $`m_t^2\mathrm{cot}^2\beta `$ and $`m_b^2\mathrm{tan}^2\beta `$. For large values of $`\mathrm{tan}^2\beta `$ this would provide a genuine possibility of measuring this fundamental parameter, as already stressed in .

The general conclusion that can be drawn from our analysis is that the genuine SUSY electroweak Sudakov logarithmic structure differs from the corresponding one met in the SM in very simple and specific ways. This suggest the following strategy.

Through the measurements of the coefficients of the DL and SL, the production of usual particles (leptons, quarks, gauge bosons) should provide global tests of the SM gauge and Higgs structure. The spirit of these tests is similar to the one which motivates the high precision tests with $`g2`$ or $`Z`$ peak measurements. If new particles, candidates for supersymmetry, exist, departures from SM predictions should appear and one should then compare with MSSM predictions for these log coefficients. This can be done both for the production of usual particles and for the production of the new states (sleptons, light or heavy squarks, charged and neutral Higgses, charginos, neutralinos). This should allow to check if the MSSM description is satisfactory or if modifications or extensions are needed (higher gauge bosons, more Higgses, …., or different forms of New Physics).

We would like to conclude by saying that even through the production of usual particles, ”virtual” supersymmetry would have a ”reality” at future accelerators, since it exhibits peculiar ”logarithmic fingerprints” that would be observable . Roughly, one might be tempted to summarize these results via a rough ”thumb rule”, sounding like:

”0,1,2,3,4….. count the logs and check supersymmetry”.