Theory of ferromagnetism driven by superexchange in dilute magnetic semiconductors

Magnetic properties of Ga$_{1-x}$Mn$_x$N are studied theoretically by employing a tight binding approach to determine exchange integrals $J_{ij}$ characterizing the coupling between Mn spin pairs located at distances $R_{ij}$ up to the 16th cation coordination sphere in zinc-blende GaN. It is shown that for a set of experimentally determined input parameters there are no itinerant carriers and the coupling between localized Mn$^{3+}$ spins in GaN proceeds via superexchange that is ferromagnetic for all explored $R_{ij}$ values. Extensive Monte Carlo simulations serve to evaluate the magnitudes of Curie temperature $T_\mathrm{C}$ by the cumulant crossing method. The theoretical values of $T_\mathrm{C}(x)$ are in quantitative agreement with the experimental data that are available for Ga$_{1-x}$Mn$_x$N with randomly distributed Mn$^{3+}$ ions with the concentrations $0.01 \leq x \leq 0.1$.


I. INTRODUCTION
Dilute magnetic insulators constitute an emerging class of magnetic semiconductors in which rather than the p-d Zener mechanism [1], ferromagnetic superexchange accounts for the coupling between diluted transition metal (TM) spins [2]. Due to its compatibility with IIInitrides that have reached the status of the most important semiconductors next to Si, particularly attractive is Ga 1−x Mn x N, in which there are no itinerant holes but nevertheless ferromagnetic spin-spin interactions are observed [2].
Recent progress in epitaxy, contamination-free processing, and (nano)characterization [2][3][4][5] allowed the preparation of Ga 1−x Mn x N films with the randomly distributed Mn 3+ ions up to x = 0.1 showing T C up to about 13 K [6] despite the absence of itinerant carriers. A high degree of crystallinity, a random distribution of the Mn ions, and a weak degree of compensation by residual donors were checked in these samples by a range of electron microscopy, synchrotron radiation, ion beam, optical, and magnetic resonance techniques [2,3,5].
The above experimental results were successfully described by the present authors within a tight binding approach (TBA) and Monte Carlo (MC) simulations [4,6] indicating that ferromagnetic superexchange accounts for ferromagnetism in Ga 1−x Mn x N with randomly distributed localized Mn 3+ spins. Within this model, a change of the Mn charge from 3+ to 2+ results in antiferromagnetic superexchange that dominates in intrinsic II-VI Diluted Magnetic Semiconductors (DMSs) such * Electronic address: csimseri@phys.uoa.gr as Cd 1−x Mn x Te [7,8] and, presumably, in Ga 1−x Mn x N containing a sizable concentration of compensating donor defects. Here we discuss briefly the theoretical approach allowing to understand ferromagnetism in Ga 1−x Mn x N.

II. THEORY
Blinowski, Kacman and Majewski [9] developed a theory of superexchange interactions between substitutional TM ions in zinc-blende semiconductor compounds. The magnetic ions were described in terms of Parmenter's generalization of the Anderson Hamiltonian for the relevant electronic configuration of the TM, taking into account the Jahn-Teller distortion, whereas the host band structure was modeled by the sp 3 s * TBA. They calculated numerically the energy of exchange interaction between Cr 2+ ions for zinc chalcogenides (ZnS, ZnSe, ZnTe) and found that the superexchange is ferromagnetic.
We realized that this theory can be adopted for Ga 1−x Mn x N [4], where holes introduced by the Mn acceptors are tightly bound, so that no insulator-to-metal transition occurs up to at least x = 0.1 [6]. Under these conditions, there are no itinerant carriers and the Mn ions assume the 3+ charge state characterized by the electron configuration identical to Cr 2+ substitutional cations in Cr-doped zinc-blende II-VI zinc chalcogenides [9]. In addition, the lack of mixed valence (all Mn ions are in the same 3+ charge state), precludes the presence of double exchange. In this situation, the superexchange accounts for the spin-spin interactions. Its sign is determined by the Anderson-Goodenough-Kanamori rules.
GaN has a wurtzite (wz) structure with a = 0.3188 nm and c = 0.5185 nm. Hence, to obtain the zinc-blende (zb) analogue with identical density of cation sites, we take the lattice parameter, a 0 = ( √ 3a 2 c) 1/3 = 0.45 nm. Here c/a ≈ 1.626, whereas for the "perfect" wz structure c/a = ( 8 3 ) 1/2 ≈ 1.633. In the equivalent zb structure, in the fcc cation sublattice, we name D "the diameter of hard touching spheres", i.e. D = √ 2 2 a 0 . In the fcc lattice, the distance of the nth nearest neighbors (NNs) up to the 16th ones [10], r n = D √ n, n ≤ 13 and r n = D √ n + 1, 14 ≤ n ≤ 16. There are no neighbors at distance D √ 14 and further away at distance D √ 30. We mention that fcc and hcp (corresponding to the wz structure) have the same number of 1st and 2nd nearest neighbors, i.e. 12 and 6, respectively, as well as at the same distances, i.e. D and D √ 2, respectively. The differences start at greater distances, e.g. fcc has 24 3rd nearest neighbors at distance D √ 3, while hcp has 2 3rd nearest neighbors at a distance D 8 3 [11]. However, overall, fcc and hcp have the same atomic packing factor π/ √ 18 ≈ 0.74. Moreover, hcp and fcc have almost identical bond percolation threshold (p b c ) and site percolation threshold (p s c ) for nearest neighbors [12]. Hence, wz-GaN can be approximated by zb-GaN in a first decent approach.

A. Tight binding approach
The calculations of the exchange energies J ij are performed for zinc-blende GaN, i.e. the cationic sublattice is fcc. The host band structure is modeled by the sp 3 s * TBA, employing the established parametrization for GaN in the cubic approximation [13]. The integrals over the Brillouin zone are performed using 2048 k-points. This guarantees that J ij are computed with an accuracy of 0.0002 K.
The magnetic ions are described in terms of the Parmenter Hamiltonian taking into account the Jahn-Teller distortion [14,15]. Since the effect of spin-orbit splitting is small in the valence band of GaN, the spin-dependent interaction between two Ga-substitutional Mn spins is described by a scalar Heisenberg coupling γ and δ denote the t 2 orbital (xy, xz, or yz) which is empty at the Mn 3+ ions i and j, respectively. k B is the Boltzmann constant. The magnitudes of J γδ ij are evaluated within the fourth-order perturbation theory in V pd for all possible orbital configurations γ and δ. Similarly to the case of Cr 2+ ions in II-VI compounds [9] the main contribution originates from quantum hopping involving occupied t 2 orbitals at the one Mn 3+ ion and the empty orbital at the other Mn 3+ ion. For the orbital configurations in question, we find that the averaged over the shell interaction is ferromagnetic at all distances.
To compare quantitatively the theoretical and experimental results, we assume a statistical distribution of directions corresponding to tetragonal Jahn-Teller distortions and determine the average value of the exchange energy J ij characterizing the coupling of Mn 3+ pairs at a given distance R ij in the fcc cation sublattice. This way we obtain the values of J ij shown in Fig. 1 together with the number of cations corresponding to particular coordination shells. The 1st set of exchange energies corresponds to 10 NNs ( Fig. 1(a) [4]). The 2nd set ( Fig. 1(b)) is obtained [6] for the 16th NNs and reducing the magnitude of e 2 from 4.8 to 4.4 eV, i.e. within its expected experimental uncertainty. We find a better agreement between theoretical and experimental T C (x) values with the 2nd set of exchange energies (see, Sec. III).

B. Monte Carlo simulations
We use the Heisenberg Hamiltonian We treat the Mn +3 spins S i , S j as classical vectors with norm S = 2. Mean values (i.e. per Mn ion) are denoted by · · ·, statistical averages [19] by · · · . At each MC sweep, we calculate the mean spin projections (l = x, y, z,) and the mean spin norm S l = N i=1 S il /N and S = N i=1 S i /N . N = N Mn is the number of Mn ions. S l = nt n=1 S l /n t and S = nt n=1 S/n t . n (n t ) denotes successive (the total number of) MC sweeps used for the statistical average. Similarly, for any integer p one could define S p = nt n=1 S p /n t . The spin susceptibility components per spin (χ S l ) and the spin susceptibility per spin (χ S ) [19] are: We use the Metropolis algorithm. In one MC sweep all Mn spins are rotated. We usually keep 2000 initial MC sweeps to thermalize the system. The typical total number of MC sweeps is 120000. Our simulation cubes (L × L × L) have typically L = 40a 0 = 18 nm, L = 50a 0 = 22.5 nm, L = 60a 0 = 27 nm. We use the Mersenne Twister (pseudo)random number generator due to its huge period and its very high order of dimensional equidistribution [20]. Usually T C is identified with the the peak of the susceptibility χ S . However, a subtle point is that in reality the results of any MC simulation depend on the size of the system, especially, for small system sizes. One avenue out of this complexity is the cumulant crossing method [21]. The fourth-order cumulant for a lattice of linear size L (3) S 2 L and S 4 L are the statistical averages of the squares and of the fourth powers of S, respectively, averages taken over systems at equilibrium at a constant temperature T . U L has a size-independent catholic fixed point, i.e. all U L (T ) curves for different L cross at T C .
In other words, one may plot U L (T ) for various L's and estimate T C from the common intersection point. We expect that, according to Eq. 3, for very low temperatures U L = 2 3 , and at very high temperatures U L = 4 9 . Although there may be some scatter at the intersection points for different pairs (L,L ′ ) especially for very small sizes, the accuracy of the fourth-order cumulant method is very satisfactory. Suppose that we plot U L (T ) and U L ′ (T ) with L ′ > L. For T < T C ⇒ U L ′ > U L and for T > T C ⇒ U L ′ < U L , i.e. the plots of U L (T ) for a number of sizes should intersect at a point or at least their pairwise intersections should be fairly close: the intersection point is good estimation of T C . A example of the fourth-order cumulant method [21] applied to our system for x = 0.03 and using the 1st set of J ij is shown in Fig. 2.

III. RESULTS
In Fig. 3 we show magnitudes of T C obtained by MC simulations using the sets of J ij values shown in Figs. 1(a) and 1(b). The computation results for the 2nd set are in remarkable agreement with the available experimental data [4,6,22]. The 2nd set of J ij values, obtained for a smaller value of e 2 , leads to improve agreement with the experimental data, and allows to compute T C (x) down to x = 0.01. As shown, the use of 10 NNs within the 2nd set (up triangles) results already in accurate magnitudes of T C . Assuming no insulator-to-metal transition, i.e. the absence of delocalized holes, we predict roomtemperature ferromagnetism for x 0.5 in Ga 1−x Mn x N with randomly distributed Mn 3+ ions.

IV. CONCLUSION
Our theoretical results on T C (x) in Ga 1−x Mn x N agree quantitatively with the measured values in the experimentally explored range 0.01 ≤ x ≤ 0.1 [4,6,22]. This agreement supports the view that ferromagnetic superexchange is the dominant coupling mechanism between Ga-substitutional Mn 3+ ions in Ga 1−x Mn x N, leading to T C ≈ 13 K at x = 0.1. The ferromagnetic character of the coupling is in accord with the Anderson-Goodenough-Kanamori rules for partly filled t 2 states of tetrahedrally  Fig. 1(a)), the 2nd set up to 16th NNs ( Fig. 1(b)), and the 2nd set up to 10th NNs, respectively. Empty symbols represent experiments results (star [22], hexagons [4], and pentagons [6]).
coordinated TM ions. Furthermore, the identical value of the exponent m = 2.2±0.2 in the dependence T c (x) ∝ x m both for ferromagnetic ordering as found in Ga 1−x Mn x N, and spin-glass freezing observed in Mn-and Co-based II-VI DMSs (in which t 2 states are entirely occupied for the majority spin direction) [7,23] verifies the scaling law [24] implying that independently of the sign of J ij , m = λ/d, where J ij ∝ R −λ ij and d is the space dimensionality. According to our theoretical model for randomly distributed Mn 3+ ions over cation sites, room-temperature ferromagnetism will appear for x 0.5, if the high-T C regime will not be shifted to even lower x values by the insulator-to-metal transition and the associated delocalization of holes supplied by Mn ions [25]. A future growth effort will show whether it is possible to obtain Ga 1−x Mn x N with merely randomly distributed Mn 3+ in a concentration x sizably exceeding 0.1.
Comparing to the ab initio data [26], our J ij magnitudes are significantly smaller and result in a stronger dependence of T C on Mn content x. This suggests that the current first principles methods overestimate coupling between TM levels and band states, which-in turn-is adequately taken into account within the Parmenter's generalization of the Anderson Hamiltonian [9] employed here. Furthermore, since within the present formalism it is possible to compute the magnitudes of J ij for virtually any distance R ij , we have been able to evaluate T C (x) down to the experimentally relevant range of x = 0.01.