Intrinsic long range antiferromagnetic coupling in dilutely V doped CuInTe$_2$

TL;DR

This study reveals an intrinsic long-range antiferromagnetic coupling in V-doped CuInTe$_2$, driven by p-d interactions and electron redistribution, with potential implications for magnetic semiconductor design.

Contribution

It uncovers a novel long-range AFM coupling mechanism in V-doped CuInTe$_2$ using DFT and a Hubbard model, advancing understanding of magnetic interactions in chalcopyrite DMS.

Findings

Discovery of intrinsic long-range AFM in V-doped CuInTe$_2$

Identification of p-d coupling as the main mechanism

Development of a Hubbard model for fundamental understanding

Abstract

Despite the various magnetic orders mediated by superexchange mechanism, the existence of a long range antiferromagnetic (AFM) coupling is unknown. Based on DFT calculations, we discovered an intrinsic long range AFM coupling in V doped CuInTe$_2$. The AFM coupling is mainly due to the $p-d$ coupling and electron redistribution along the interacting chains. The relatively small energy difference between $p$ and $d$ orbitals and the small energy difference between d orbitals of the dopants and that of stepping stone sites can enhance the stability of this AFM configuration. A multi-bands Hubbard model was proposed to provide fundamental understanding to the long range AFM coupling in chalcopyrite diluted magnetic semiconductors(DMS).

Tables

Table 1. Table 1: The positions of the second V atom.

NN	1	2	3	4	5	6
position	(2, 8)	(2, 6)	(5, 1) or (5, 3)	(5, 6)	(3, 8)	(6, 1)
V-V distance(Å)	4.38	6.19	7.58	8.75	9.79	11.59

Table 2. Table 2: The character table of T d subscript 𝑇 𝑑 T_{d} group.

$T_d$	$E$	$8C_3$	$3C_2$	$6{\sigma }_d$	$6S_4$
$A_1$	1	1	1	1	1
$A_2$	1	1	1	-1	-1
$E$	2	-1	2	0	0
$T_1$	3	0	-1	-1	1
$T_2$	3	0	-1	1	-1
${\Gamma }_p$	3	0	-1	1	-1
${\Gamma }_d$	5	-1	1	1	-1

Equations

$$E_f=E(\text{doped})-n_{\text{V}}{\mu }_{\text{V}}-E(\text{undoped})+n_{\text{In}}{\mu }_{\text{In}}.$$

$$\hat{V}={\hat{V}}_0+{\hat{V}}_{sf}+{\hat{V}}_{ph}.$$

$$\begin{array}{rl}{\hat{V}}_0& =u{\hat{n}}^2-v{\hat{m}}_z^2-(u-v)\hat{n}+8v\sum _{\alpha }{\hat{n}}_{\alpha \uparrow }{\hat{n}}_{\alpha \downarrow },\\ {\hat{V}}_{sf}& =-2v\sum _{\alpha \ne \beta ,\sigma }{\hat{c}}_{\alpha \sigma }^{†}{\hat{c}}_{\alpha ,-\sigma }{\hat{c}}_{\beta ,-\sigma }^{†}{\hat{c}}_{\beta \sigma },\\ {\hat{V}}_{ph}& =2v\sum _{\alpha \ne \beta }({\hat{n}}_{\alpha \beta }{)}^2\\ & =-2v\sum _{\alpha \ne \beta ,\sigma }{\hat{c}}_{\alpha \sigma }^{†}{\hat{c}}_{\alpha ,-\sigma }^{†}{\hat{c}}_{\beta \sigma }{\hat{c}}_{\beta ,-\sigma }.\end{array}$$

$$\begin{array}{rl}u& =\frac{1}{2}U-\frac{1}{4}J+\frac{5}{2}\Delta J,\\ v& =\frac{1}{4}J-\frac{3}{2}\Delta J.\end{array}$$

$$\begin{array}{rl}\hat{H}& =\sum _{⟨i,j⟩}\sum _{\alpha \beta \sigma }{t}_{\alpha \beta }^{ij}{\hat{c}}_{i\alpha \sigma }^{†}{\hat{c}}_{j\beta \sigma }+h.c.\\ & +\sum _{i\alpha \sigma }{ϵ}_{i\alpha }{\hat{c}}_{i\alpha \sigma }^{†}{\hat{c}}_{i\alpha \sigma }+\sum _i{\hat{V}}_0^i\end{array}$$

$$\begin{array}{r}{E}_{\uparrow \downarrow }=E_{0\uparrow \downarrow }-\frac{18t^2}{E_{1\uparrow \downarrow }-E_{0\uparrow \downarrow }}-\frac{18t^4}{(E_{1\uparrow \downarrow }-E_{0\uparrow \downarrow }{)}^2(E_{2\uparrow \downarrow }-E_{0\uparrow \downarrow })}\end{array}$$

$$\begin{array}{rl}{E}_{\uparrow \uparrow }& =E_{0\uparrow \uparrow }-\frac{18t^2}{E_{1\uparrow \uparrow }-E_{0\uparrow \uparrow }}\\ & -\frac{9t^4}{(E_{1\uparrow \uparrow }-E_{0\uparrow \uparrow }{)}^2}(\frac{1}{E_{2\uparrow \uparrow }-E_{0\uparrow \uparrow }}+\frac{1}{E_{2\uparrow \uparrow }^{\prime }-E_{0\uparrow \uparrow }}).\end{array}$$

$$E_{\uparrow \uparrow }-E_{\uparrow \downarrow }=24v_A-\frac{72t^4{v}_A}{(\Delta E_{10}{)}^2(\Delta E_{20}{)}^2}$$

$$\begin{array}{rl}\hat{n}& =\sum _{\alpha \sigma }{\hat{c}}_{\alpha \sigma }^{†}{\hat{c}}_{\alpha \sigma }=\sum _{\sigma }{\hat{n}}_{\alpha },\\ {\hat{n}}_{\alpha \beta }& =\sum _{\sigma }{\hat{c}}_{\alpha \sigma }^{†}{\hat{c}}_{\beta \sigma },\\ \hat{\mathbf{m}}& =\sum _{\alpha \sigma {\sigma }^{\prime }}{\hat{c}}_{\alpha \sigma }^{†}{\mathbf{\tau }}_{\sigma {\sigma }^{\prime }}{\hat{c}}_{\alpha {\sigma }^{\prime }},\\ {\mathbf{\tau }}_{\sigma {\sigma }^{\prime }}& =({\tau }_{\sigma {\sigma }^{\prime }}^x,{\tau }_{\sigma {\sigma }^{\prime }}^y,{\tau }_{\sigma {\sigma }^{\prime }}^z).\end{array}$$

$$\begin{array}{rl}:{\hat{n}}^2:& =\sum _{\alpha \beta \sigma {\sigma }^{\prime }}{\hat{c}}_{\alpha \sigma }^{†}{\hat{c}}_{\beta {\sigma }^{\prime }}^{†}{\hat{c}}_{\beta {\sigma }^{\prime }}{\hat{c}}_{\alpha \sigma }\\ & =-\sum _{\alpha \beta \sigma {\sigma }^{\prime }}{\hat{c}}_{\alpha \sigma }^{†}({\delta }_{\alpha \beta }{\delta }_{\sigma {\sigma }^{\prime }}-{\hat{c}}_{\alpha \sigma }{\hat{c}}_{\beta {\sigma }^{\prime }}^{†}){\hat{c}}_{\beta {\sigma }^{\prime }}\\ & =-\hat{n}+{\hat{n}}^2,\end{array}$$

$$\begin{array}{rl}:{\hat{n}}_{\alpha \beta }^2:& =\sum _{\sigma {\sigma }^{\prime }}{\hat{c}}_{\alpha \sigma }^{†}{\hat{c}}_{\alpha {\sigma }^{\prime }}^{†}{\hat{c}}_{\beta {\sigma }^{\prime }}{\hat{c}}_{\beta \sigma }\\ & =-\sum _{\sigma {\sigma }^{\prime }}{\hat{c}}_{\alpha \sigma }^{†}({\delta }_{\alpha \beta }{\delta }_{\sigma {\sigma }^{\prime }}-{\hat{c}}_{\beta \sigma }{\hat{c}}_{\alpha {\sigma }^{\prime }}^{†}){\hat{c}}_{\beta {\sigma }^{\prime }}\\ & =-{\hat{n}}_{\alpha \beta }{\delta }_{\alpha \beta }+{\hat{n}}_{\alpha \beta }^2.\end{array}$$

$$\begin{array}{rl}{\hat{n}}_{\alpha }^2& =({\hat{n}}_{\alpha \uparrow }+{\hat{n}}_{\alpha \downarrow }{)}^2\\ & ={\hat{n}}_{\alpha \uparrow }^2+{\hat{n}}_{\alpha \downarrow }^2+{\hat{n}}_{\alpha \uparrow }{\hat{n}}_{\alpha \downarrow }+{\hat{n}}_{\alpha \downarrow }{\hat{n}}_{\alpha \uparrow }\\ & ={\hat{n}}_{\alpha \uparrow }+{\hat{n}}_{\alpha \downarrow }+2{\hat{n}}_{\alpha \uparrow }{\hat{n}}_{\alpha \downarrow }.\end{array}$$

$$\begin{array}{rl}\sum _{\alpha \beta }:{\hat{n}}_{\alpha \beta }^2:& =-\sum _{\alpha }{\hat{n}}_{\alpha }+\sum _{\alpha \beta }{\hat{n}}_{\alpha \beta }^2\\ & =-\hat{n}+\sum _{\alpha }{\hat{n}}_{\alpha }^2+\sum _{\alpha \ne \beta }{\hat{n}}_{\alpha \beta }^2\\ & =2\sum _{\alpha }{\hat{n}}_{\alpha \uparrow }{\hat{n}}_{\alpha \downarrow }+\sum _{\alpha \ne \beta }{\hat{n}}_{\alpha \beta }^2.\end{array}$$

$$\begin{array}{rl}{\hat{\mathbf{m}}}^2& =\sum _{\alpha \sigma {\sigma }^{\prime }}\sum _{\beta \zeta {\zeta }^{\prime }}{\mathbf{\tau }}_{\sigma {\sigma }^{\prime }}\cdot {\mathbf{\tau }}_{\zeta {\zeta }^{\prime }}{\hat{c}}_{\alpha \sigma }^{†}{\hat{c}}_{\alpha {\sigma }^{\prime }}{\hat{c}}_{\beta \zeta }^{†}{\hat{c}}_{\beta {\zeta }^{\prime }}\\ & =2\sum _{\alpha \beta {\sigma }^{\prime }\sigma }{\hat{c}}_{\alpha \sigma }^{†}{\hat{c}}_{\alpha {\sigma }^{\prime }}{\hat{c}}_{\beta \sigma }^{†}{\hat{c}}_{\beta {\sigma }^{\prime }}-{\hat{n}}^2.\end{array}$$

$$\begin{array}{rl}{\hat{m}}_z& =\sum _{\alpha }({\hat{n}}_{\alpha \uparrow }-{\hat{n}}_{\alpha \downarrow }),\\ {\hat{n}}^2+{\hat{m}}_z^2& =2\sum _{\alpha \beta }({\hat{n}}_{\alpha \uparrow }{\hat{n}}_{\beta \uparrow }+{\hat{n}}_{\alpha \downarrow }{\hat{n}}_{\beta \downarrow }).\end{array}$$

$${\hat{m}}^2={\hat{m}}_z^2+2\sum _{\alpha \beta \sigma }{\hat{c}}_{\alpha \sigma }^{†}{\hat{c}}_{\alpha ,-\sigma }{\hat{c}}_{\beta ,-\sigma }^{†}{\hat{c}}_{\beta \sigma }.$$

$$\begin{array}{rl}& 2\sum _{\alpha }({\hat{c}}_{\alpha \uparrow }^{†}{\hat{c}}_{\alpha \downarrow }{\hat{c}}_{\alpha \downarrow }^{†}{\hat{c}}_{\alpha \uparrow }+{\hat{c}}_{\alpha \downarrow }^{†}{\hat{c}}_{\alpha \uparrow }{\hat{c}}_{\alpha \uparrow }^{†}{\hat{c}}_{\alpha \downarrow })\\ & =2\sum _{\alpha }[{\hat{n}}_{\alpha \uparrow }(1-{\hat{n}}_{\alpha \downarrow })+{\hat{n}}_{\alpha \downarrow }(1-{\hat{n}}_{\alpha \uparrow })]\\ & =2\hat{n}-4\sum _{\alpha }{\hat{n}}_{\alpha \uparrow }{\hat{n}}_{\alpha \downarrow },\end{array}$$

$${\hat{\mathbf{m}}}^2={\hat{m}}_z^2+2\sum _{\alpha \ne \beta ,\sigma }{\hat{c}}_{\alpha \sigma }^{†}{\hat{c}}_{\alpha ,-\sigma }{\hat{c}}_{\beta ,-\sigma }^{†}{\hat{c}}_{\beta \sigma }+2\hat{n}-\sum _{\alpha }{\hat{n}}_{\alpha \uparrow }{\hat{n}}_{\alpha \downarrow }.$$

$$\begin{array}{rl}:{\hat{\mathbf{m}}}^2:& =2\sum _{\alpha \beta \sigma {\sigma }^{\prime }}{\hat{c}}_{\alpha \sigma }^{†}{\hat{c}}_{\beta {\sigma }^{\prime }}^{†}{\hat{c}}_{\beta \sigma }{\hat{c}}_{\alpha {\sigma }^{\prime }}-:{\hat{n}}^2:\\ & =-2\sum _{\alpha \beta \sigma {\sigma }^{\prime }}{\hat{c}}_{\alpha \sigma }({\delta }_{\alpha \beta }-{\hat{c}}_{\alpha {\sigma }^{\prime }}{\hat{c}}_{\beta {\sigma }^{\prime }}^{†}){\hat{c}}_{\beta \sigma }-{\hat{n}}^2+\hat{n}\\ & =-2\sum _{\alpha \sigma {\sigma }^{\prime }}{\hat{c}}_{\alpha \sigma }^{†}{\hat{c}}_{\alpha \sigma }+{\hat{\mathbf{m}}}^2+\hat{n}\\ & =-3\hat{n}+{\hat{\mathbf{m}}}^2.\end{array}$$

$$\begin{array}{rl}\hat{V}& =\frac{1}{2}[(U-\frac{1}{2}J+5\Delta J):{\hat{n}}^2:-\frac{1}{2}(J-6\Delta J):{\hat{\mathbf{m}}}^2:\\ & +(J-6\Delta J)\sum _{\alpha \beta }:({\hat{n}}_{\alpha \beta }{)}^2:]\end{array}$$

$$\begin{array}{rl}U& =V_{23,23},\\ J& =\frac{1}{2}(V_{15,51}+V_{23,32}),\\ \Delta J& =V_{15,51}-V_{23,32}.\end{array}$$

$$\hat{V}={\hat{V}}_0+{\hat{V}}_{sf}+{\hat{V}}_{ph}.$$

$$\begin{array}{rl}{\hat{V}}_0& =u{\hat{n}}^2-v{\hat{m}}_z^2-(u-v)\hat{n}+8v\sum _{\alpha }{\hat{n}}_{\alpha \uparrow }{\hat{n}}_{\alpha \downarrow },\\ {\hat{V}}_{sf}& =-2v\sum _{\alpha \ne \beta ,\sigma }{\hat{c}}_{\alpha \sigma }^{†}{\hat{c}}_{\alpha ,-\sigma }{\hat{c}}_{\beta ,-\sigma }^{†}{\hat{c}}_{\beta \sigma },\\ {\hat{V}}_{ph}& =2v\sum _{\alpha \ne \beta }({\hat{n}}_{\alpha \beta }{)}^2\\ & =-2v\sum _{\alpha \ne \beta ,\sigma }{\hat{c}}_{\alpha \sigma }^{†}{\hat{c}}_{\alpha ,-\sigma }^{†}{\hat{c}}_{\beta \sigma }{\hat{c}}_{\beta ,-\sigma }.\end{array}$$

$$\begin{array}{rl}u& =\frac{1}{2}U-\frac{1}{4}J+\frac{5}{2}\Delta J,\\ v& =\frac{1}{4}J-\frac{3}{2}\Delta J.\end{array}$$

$$\begin{array}{rl}\hat{H}& ={\hat{H}}_t+{\hat{H}}_e+{\hat{H}}_{int},\\ & =\sum _{⟨i,j⟩}\sum _{\alpha \beta \sigma }{t}_{\alpha \beta }^{ij}{\hat{c}}_{i\alpha \sigma }^{†}{\hat{c}}_{j\beta \sigma }+h.c.\\ & +\sum _{i\alpha \sigma }{ϵ}_{i\alpha }{\hat{c}}_{i\alpha \sigma }^{†}{\hat{c}}_{i\alpha \sigma }+\sum _i{\hat{V}}^i\end{array}$$

$${\Gamma }_p$$

$${\Gamma }_d$$

$$\begin{array}{rl}{E}_{0\uparrow \downarrow }& =4{ϵ}_{A1}+6{ϵ}_{A2}+12{ϵ}_B+4{ϵ}_{C1}\\ & +40u_A-40v_A+60u_B+60v_B+12u_C+20v_C.\end{array}$$

$$\begin{array}{rl}{E}_{1\uparrow \downarrow }& =E_{0\uparrow \downarrow }+8u_C-10u_B-10v_B-{\Delta }_{BC}+{\delta }_C,\\ E_{2\uparrow \downarrow }& =E_{0\uparrow \downarrow }+8u_C-8u_A+8v_A-{\Delta }_{AC}+{\delta }_C-{\delta }_A.\end{array}$$

$$E_{\uparrow \downarrow }=E_{0\uparrow \downarrow }-\frac{18t^2}{E_{1\uparrow \downarrow }-E_{0\uparrow \downarrow }}-\frac{18t^4}{(E_{1\uparrow \downarrow }-E_{0\uparrow \downarrow }{)}^2(E_{2\uparrow \downarrow }-E_{0\uparrow \downarrow })}.$$

$$\begin{array}{rl}{E}_{0\uparrow \uparrow }& =4{ϵ}_{A1}+6{ϵ}_{A2}+12{ϵ}_B+4{ϵ}_{C1}\\ & +40u_A-16v_A+60u_B+60v_B+12u_C+20v_C.\end{array}$$

Full text

I Introduction

Anti-ferromagnetic(AFM) order is essential in many spintronic devices such as tunnel junctions, spin valves, Hall devices, AFM/FM bilayers.Wang et al. (2017); Baltz et al. (2018); Fukami et al. (2016). Compared to ferromagnets, AFM materials have high information storage densities, high reading and writing speed, low power consumption, and tetra-hertz spin oscillation frequenciesCheng et al. (2014). AFM order has been observed in various metallic alloysZajac et al. (2001), in which magnetism tuning can be challenging. On the other hand, AFM order in diluted magnetic semiconductors(DMS) might be relatively easy to tune because prototypical semiconductor properties are preserved.

In DMS, a relatively low concentration of transition metal atoms are doped, compared to magnetic semiconductors. This unique property makes tuning of magnetic order in DMS easier than typical magnetic semiconductors. In early theoretical investigations, it has been found that the magnetic moment of magnetic dopants tend to form antiparallel coupling in large magnetic doping concentrations and small simulation cells in conventional semiconductors such as V doped GaAs under carrier free conditionBelhadji et al. (2007). However, such a coupling may not exist in DMS due to the limited dopant concentrations. Whether a long range AFM coupling exists or not in DMS is still an open question.

Generally speaking, superexchangeGoodenough (1955) and coupling between localized spins and carriersRuderman and Kittel (1954); Kasuya (1956); Yosida (1957) are the major mechanisms to induce an AFM coupling in materials. Here, the magnetic configurations are often in the type of ABA chains, where A sites are magnetic elements. In DMS, the short range AFM coupling leads to an intrinsic difficulty to realize long range AFM order. Therefore, to realize the long range AFM order, there must be other mechanisms involved. Recently, topological surface statesZhang et al. (2009) and stepping stone mechanism have been proposed to explain the ferromagnetic coupling in Cr doped (Sb, Bi)2Te3Chan et al. (2018). Stepping stone mechanism is based on spin polarized s lone pair states of cations in C sites of spin chains of ABCBA type. Whether stepping stone mechanism can be extended to spin polarized d states that may stabilize the long range AFM coupling is still unknown, except one report on a possible long range AFM coupling in Mn doped LiZnAsZhang et al. (2018), which is relatively unstable. Also, no quantitative model has been constructed to describe the stepping stone mechanism and fundamental physical understanding is poor.

In this letter we proposed a new kind of long range AFM coupling of magnetic dopants in transition metal(TM) chalcopyrite, CuInTe2, which is an intrinsic narrow band gap semiconductor. The long range is defined as ABCBA type or longer spin chains. An extended super-exchange model was constructed to explain the underlying mechanism which gives rise to the intrinsic AFM coupling in DMS. This microscopic mechanism may lead to the realization of macroscopic AFM order, which can be checked by large scale simulations that is beyond the capability of DFT.

II Methods

All calculations were performed using projected augmented plane wave methodBlöchl (1994) and density functional theory with Perdew-Burke-Ernzerhof generalized gradient approximation(GGA) Perdew et al. (2009) as implemented in VASP codeKresse and Furthmüller (1996). Atoms were relaxed with force tolerance of 0.01 eVÅ*-1*. A plane wave energy cut-off of 300 eV was used in all calculations. The Brillouin zone integrations were performed by using $\Gamma$ centered $5\times 5\times 5$ k-points grid. A strong correlation effect was considered for transition metal, and DFT+$U$ method was usedDudarev et al. (1998). We chose the onsite Coulomb interaction parameter $U=4.70$ eV and onsite exchange interaction $J=0.70$ eV, so that effective parameter $U_{\text{eff}}=4.00$ eV, as suggested in Ref.Zhang et al. (2013). Convergence test about energy cut-off, number of k-points mesh, cell size were performedsm . In order to explore the short range and long range magnetic coupling, we substituted two In atoms with two V atoms with increasing distance. The V doped CuInTe2 was simulated in $3\times 3\times 1$ supercells based on conventional cell with 144 atoms. The doped compound is CuIn1-xVxTe2 with $x$ equals 5.6%. Convergence test was also done in $2\times 2\times 1$ supercell, see supplementary information for detailssm . Stability of V doped CuInTe2 is characterized by its formation energy, which is defined as

[TABLE]

where $E$ are total energies of CuIn1-xVxTe2 and CuInTe2, $\mu$ are chemical potentials of V and In, $n$ is the number of corresponding atoms. In our study, different configurations have the same number of In and V atoms, so difference of formation energies between different configurations cancels the contribution of chemical potentials. And we define this difference between each configuration and that of the most stable one as relative formation energy, which equals the difference of total energies between them.

III Results

The relative formation energies of V doped CuInTe2 as a function of different neighboring configurations were calculated from fully relaxed $3\times 3\times 1$ supercell, shown in Fig.1. The results in $2\times 2\times 1$ supercell and spin-orbit coupling(SOC) effect were given in supplementary materialsm . As shown in the FIG.1, the second NN AFM state is the global minimum, with an AFM-FM energy difference 5.81 meV. The relative formation energy difference between the first NN and second NN are similar in both $2\times 2\times 1$ and $3\times 3\times 1$ supercells. Hence, subsequent calculation results of density of states and spin texture were illustrated in $2\times 2\times 1$ supercell.

To understand the mechanism for such kind of long range carrier free AFM interaction, we calculated the projected density of states(pDOS) for second NN configuration, as shown in FIG.2 The V shows clearly polarized $p-d$ hybridized states right above the Fermi level. And the polarized $p$ states of V atoms mainly come from the hybridization of neighboring Te atoms. These polarized states have similar energy and shape as the polarized $p$ states of Te atoms. The Cu atoms, however, show no spin polarization. This is because Cu atom is at the middle of the chain, and two halves of the chain polarized it with same magnitude but opposite direction.

In order to see the spin polarization more clearly, we calculated the polarized spin density of AFM second NN, as shown in FIG.3. The V atom at the body center and another two at the face center are connected by three V-Te-Cu-Te-V chains and one V-Te-In-Te-V chain. We can see clearly polarized Cu $d$ orbitals in the middle, which serve as ”stepping stones” to pass on the magnetism along the chain. In addition, long range AFM coupling was discovered in V doped AgGaTe2 and AgInTe2, see supplementary information for details.

IV Model

In this section, a multi-bands Hubbard model is proposed to study the mechanism of long range AFM coupling in an A-B-C-B-A like chain which possesses the same geometry as that of the V-Te-Cu-Te-V chain, which is the building block of this AFM structure. For simplicity of theoretical modeling, we will ignore the effect of SOC in the analysis below. Nevertheless, SOC effect can still be treated using perturbation theory to achieve more accurate results, which is out the scope of this paper. We constructed various spin configurations with parallel or anti-parallel spins on $e_{g}$ states at the A(V) sites. Effective Hamiltonians in low energy subspace were derived for these two cases, and configurations with anti-parallel spins show relatively lower energy than that with parallel spins. And this accounts for an effective AFM coupling between V atoms. Finally, an estimation of the effective magnetic coupling was given.

We use the expression of onsite electron electron interaction for p and d orbitals in cubic crystal field following the result from Coury et al. (2016). We neglected the quadrupole moment terms and derived sm :

[TABLE]

where $\hat{V}_{0}$, spin flipping term $\hat{V}_{sf}$ and pair hopping term $\hat{V}_{ph}$ are given by:

[TABLE]

where $\alpha,\beta$ are indices of local orbitals, $\sigma,\sigma^{\prime}$ are indices of spin. $\hat{n}$ is the electron number operator, $\hat{m_{z}}$ is the $z$ component of the magnetic moment operator. $\hat{n}_{\alpha\beta}=\sum_{\sigma}\hat{c}_{\alpha\sigma}^{\dagger}\hat{c}_{\beta\sigma}$ is the onsite hopping between two orbitals. Here site label is omitted for simplicity. The parameters $u,v$ in equation (22) are defined as:

[TABLE]

For $d$ orbitals, $U$ is the Coulomb interaction between $t_{2g}$ orbitals, $J$ is the average exchange splitting of $e_{g}$ and $t_{2g}$ orbitals and $\Delta J$ is the difference of exchange splitting between $e_{g}$ and $t_{2g}$ orbitals following the definitions in Coury et al. (2016). While for $p$ orbitals, $U$ and $J$ are the Coulomb interaction and exchange splitting of three $p$ orbitals, and $\Delta J=0$.

Then the Hamiltonian of multi-orbital Hubbard model under tight binding approximation is:

[TABLE]

where $i,j=A,B,C$ are atomic sites and $t_{\alpha\beta}^{ij}$ are hopping integrals between nearest neighbor atoms. $\sum_{i}\hat{V}^{i}\approx\sum_{i}\hat{V}_{0}^{i}$ is the interaction according to equation (22) on each site, and the spin flip and pair hopping have been ignored by approximation. For detailed discussions, please see the latter part of this letter.

A and C atoms have local Td symmetry, so local $d$ orbital splits to 2 fold degenerated $e_{g}$ states and 3 fold degenerated $t_{2g}$ states. Te $p$ orbitals interact strongly with nearby $t_{2g}$(A, C) states which have the same symmetry and form $\sigma$ bonds. Te $p$ orbitals form only weaker $\pi$ bonds with nearby $e_{g}$ (A, C) states, which play less important role near the Fermi level Jaffe and Zunger (1984); Yoodee et al. (1984), and we ignored these terms in this model. The signs of the hopping integrals depend on the geometry of the chain. Using the Slater-Koster matrix elements Slater and Koster (1954), the signs of hopping integrals can be determinedsm . The magnitudes of hopping were approximated to be the same between A, B and B, C with $t=V_{pd\pi}/\sqrt{3}$.

Starting from the neutral state, different configurations were obtained by nearest neighbor hopping among $p$ orbitals of B and $t_{2g}$ (A, C). The $e_{g}$ (A, C) were fixed in this model. The relative position of each atomic level were taken to be the same as V(A), Te(B), Cu(C), and this relation can as well be applied to compounds with elements one row lower or higher. The total energies of this chain with parallel and anti-parallel initial V spins were calculated using equation (5) sm .

The lowest energy states with parallel($\psi_{0\uparrow\uparrow}$) and anti-parallel($\psi_{0\uparrow\downarrow}$) $e_{g}$ spins at site A are shown in FIG.4. In state $\psi_{0\uparrow\downarrow}$, two V(A) atoms have half-filled $3d$ orbitals with anti-parallel spins, two Te(B) atoms have fully filled $5p$ orbitals and Cu(C) atom has empty $t_{2g}$ orbitals and fully filled $e_{g}$ orbitals. The spin configuration of two A atoms in this state is $3d_{\uparrow}^{5}$, $3d_{\downarrow}^{5}$($S=\frac{5}{2}$) obeying the Hund’s rule. While in state $\psi_{0\uparrow\uparrow}$, electron filling is the same at site B, C, but one V(A) atom has $3d_{\uparrow}^{5}$ filling and another has $e_{g\uparrow}^{2}t_{2g\downarrow}^{3}$ filling. The energy of $\psi_{0\uparrow\uparrow}$ is $24v_{A}$ larger than that of $\psi_{0\uparrow\downarrow}$. This is because $\psi_{0\uparrow\downarrow}$ state has larger local magnetic moment and results in a smaller energy. Although hopping directly between A and C is not allowed in this model, it can happen by virtual hopping. In this process, one spin on the $p$ state at site B hops to $t_{2g}$ state at site C first, then one spin on the $t_{2g}$ state at site A hops to B site, and vise versa. This is a typical process that happens in super-exchange mechanismAnderson (1950). We denote the intermediate states and states after one virtual hopping as $\psi_{1\uparrow\downarrow},\psi_{1\uparrow\uparrow}$ and $\psi_{2\uparrow\downarrow},\psi_{2\uparrow\uparrow}$. In the parallel case, two A atoms are not symmetric, so we denote state after one virtual hopping of an up(down) spin from $t_{2g}$ orbitals at site A as $\psi_{2\uparrow\uparrow}(\psi_{2\uparrow\uparrow}^{\prime})$.

Relatively higher energy states are those with smaller local magnetic moment at site A, the $t_{2g}$ spins at site A can be different. There are huge energy barriers between these states and the two lowest energy states, and we focus only on the lowest energy states and their intermediate states, so these states are not considered in the analysis below. Starting from lowest energy configurations, higher order configurations can be obtained by nearest neighbor hopping. We can apply perturbation theory to calculate the effective ground state energy. Applying the fourth order perturbation theory, we obtain:

[TABLE]

and

[TABLE]

where all $E$ terms refer to the energies of states defined by the subscripts. The coefficients of the fourth order terms come from the degeneracy of $d$ orbitals of C atom. The fourth order terms are different for parallel and anti-parallel cases because two A atoms are inequivalent for the parallel case. This gives us the energy difference between parallel and anti-parallel initial V spins with forth order correction is:

[TABLE]

where $\Delta E_{10}=E_{1\uparrow\downarrow}-E_{0\uparrow\downarrow}=E_{1\uparrow\uparrow}-E_{0\uparrow\uparrow}$, $\Delta E_{20}=E_{2\uparrow\downarrow}-E_{0\uparrow\downarrow}=E_{2\uparrow\uparrow}-E_{0\uparrow\uparrow}$.

We found that anti-parallel initial spins would have lower energy than that of the parallel one. So the system prefers long range AFM interaction. This long range AFM interaction is mediated by low-lying $d$ orbitals of Cu atoms. If we replace the Cu atoms on the chain by Ag atoms, the low energy configurations does not change much, while $\Delta E_{10}$ and $\Delta E_{20}$ become larger due to higher $d$ level of Ag atoms. The energy difference becomes larger according to (8). Remember that two V atoms are connected by six chains, three to the left, three to the right. To further illustrate the stepping stone mechanism mediated by $d$ states, we replaced the six Cu atoms by Na atoms at the stepping stone sites and found that the energy difference between FM and AFM configurations is dramatically reduced to 1.6 meV. This small residue energy is probably due to the magnetism mediated by In atoms. We also replaced the six Cu atoms by Ag atoms at stepping stone sites and found that the energy difference between FM and AFM configurations changes from 7.57 meV for the Cu case to 8.89 meV for the Ag case, which confirms the validity of our model.

The actual energy difference is small due to the following reasons: 1) The two V atoms in this system are actually connected by three V-Te-Cu-Te-V chains. Three Cu atoms in the middle provide more hopping channels, which will enlarge the coefficients of the second term in equation (8); 2) configurations with different local magnetic moments also have low energies. For example, energy of configurations with $3d_{\uparrow}^{4}3d_{\downarrow}^{1}$, $3d_{\downarrow}^{4}3d_{\uparrow}^{1}$ (spin distribution of V atoms is $S=\frac{3}{2}$) is 32$v_{A}$ larger than that of lowest energy state. Although these configurations are ignored due to huge energy barriers, see supplementary information for details, they also help to reduce the magnetic coupling. And the ground state will be the superposition of these states, so that the actual magnetic moments of V atoms are smaller than that of the lowest energy configuration, which has the largest magnetic moments. Also, the effect of spin flipping and pair hopping become more important when $S<\frac{5}{2}$, because in this case there are both up spins and down spins at site A. Spin flipping happens when $d$ orbitals of A are half-filled and pair hopping happens when there are both doubly filled and empty orbitals. And spin flipping effect of ”stepping stones” induce a ferromagnetic superexchange when there is a 90 degree twisted chain structureGoodenough (1955).

Our model gives clues to finding stable long range AFM order in DMS, and several guidelines can be formed: 1) To search for a long range AFM order independent of RKKY interaction, the host cell should be carrier free; 2) There should be low lying states along the path between magnetic dopants which serve as ”stepping stones” to pass on the spin exchange process. Such states are relatively easy to be magnetized due to a large local Coulomb repulsion that may induce charge redistribution along the chain. As a result, long range magnetic order emerges with the assistance of the charge redistribution; 3) Also host cell should consist of large atoms, since they may suppress neighboring superexchange and favor the long range interactions.

V Conclusion

To conclude, we have discovered an intrinsic long range AFM structure in V doped CuInTe2. This AFM coupling can not be explained by RKKY interaction since there are no carriers in this system. In addition, this AFM coupling can not be explained by van Vleck paramagnetism based on band inversion in topological insulators because system is topologically trivial. It can rather be explained by an extended superexchange mechanism in a A-B-C-B-A like chain structure, which we named as ”stepping stone” mechanism. Energy difference between each two of the three atoms, local Coulomb interaction, local exchange splitting, local orbital degeneracy and asymmetry between two halves of the chain help to induce an effective AFM coupling between two magnetic dopants. A multi-orbitals Hubbard model was proposed and confirmed by further DFT calculations. Our model can also be applied to other long range chain structures with different types of localized electronic states locating at the stepping stone sites.

VI Acknowledgement

We are grateful for the financial support of Chinese University of Hong Kong (CUHK) (Grant No.4053084), University Grants Committee of Hong Kong (Grant No. 24300814), and start-up funding of CUHK.

VII Convergence test and calculation in $\text{AgGaTe}_{2}$ and $\text{AgInTe}_{2}$

All calculations were performed using projected augmented plane wave methodBlöchl (1994) and density functional theory with Perdew-Burke-Ernzerhof generalized gradient approximation(GGA) Perdew et al. (2009) as implemented in VASP codeKresse and Furthmüller (1996). The V doped CuInTe2 was simulated in $2\times 2\times 1$ supercell and $3\times 3\times 1$ supercell. Atoms were relaxed with force tolerance of 0.01 eVÅ*-1*. A plane wave energy cut-off of 300 eV was used in all calculations. The Brillouin zone integrations were performed by using $\Gamma$ centered $5\times 5\times 5$ k-points grid. The convergence test of energy of the system using different number of k points was shown in figure 5.

Similar calculations of the relative formation energy of different configurations were done in AgGaTe2 and AgInTe2 in $2\times 2\times 1$ supercells, as shown in figure 6. Similar as in CuInTe2, V atoms in AgInTe2 have a lowest energy in the second nearest neighbor configuration, while lowest energy configuration in AgGaTe2 is the forth nearest neighbor.

VIII Calculations in $2\times 2\times 1$ supercell

The calculation setup in $2\times 2\times 1$ supercell was the same as stated in the main text. An In atom at body center of the supercell was substituted with a V atom. Another V atom replaces an In atom at various neighboring sites to the center V atom, from the first nearest neighbor(NN) to the forth NN, as shown in Fig.7. Cu, V and In atoms have approximately local $T_{d}$ symmetry in the host cell, which is a property of chalcopyrite. The point group of this supercell is $\bar{4}2d$ ($D_{2d}$). Under this symmetry, the third NN has two nonequivalent configurations, while other three have only one configuration for each.

Due the limitation of calculation resources, SOC calculation was only done in $2\times 2\times 1$ supercell. The relative formation energy without and with SOC is shown in FIG.8.

Although the effect of spin orbit coupling (SOC) is ignored in the theoretical model in the main text, SOC may induce further splitting of energy levels of heavy atoms like Te. Due to the SOC effect, the $p$ level of Te will become lower, while this effect at Cu and In atoms is not so obvious. Hence, according to the perturbation theory, smaller energy difference will cause a larger magnetic coupling. The results of DFT calculation with SOC effect considered confirmed this point. At the second nearest neighbor, AFM has energy 7.57 meV lower than that of FM without SOC effect, while this value becomes 10.74 meV with SOC considered. Although the magnetic coupling of the first NN changes from FM to AFM once SOC is considered, the general trends of formation energy as functions of neighboring sites are similar.

IX Calculations in $3\times 3\times 1$ supercell

The setup is the same as above, but there are more NN configurations. In FIG.10, In atoms are labeled by ($n$, $m$), where $n$ is the label of cell and $m$ is the label of In atom. In the magnetic doping process, two In atoms are replaced by two V atoms. Due to the symmetry in chalcopyrite, we fix one V atom at position (1, 6). The position of the second V atom is shown in TAB.1.

The setup of calculation in $3\times 3\times 1$ supercell was the same as that of $2\times 2\times 1$ supercell. Static calculation with SOC was performed to qualitatively check the effect of SOC. The results without and with static SOC calculation were shown in FIG.10. The general trends of the formation energy and energy difference between FM and AFM are similar. We see that the second NN is still the most stable configuration in larger supercell without SOC, while the stability of the sixth NN is comparable with that of the second NN. The formation energy is determined both by magnetic coupling strength and the local stress that may slightly vary with or without SOC. Nevertheless, the stabilization of long range configurations versus the first nearest neighboring configuration is valid for both setup.

X Normal ordering of local operators

The electron number operator, onsite hopping and local magnetic moment operator at each site are defined as:

[TABLE]

where $\alpha,\beta$ are indices of local orbitals, $\sigma,\sigma^{\prime}$ are indices for spins and $\tau^{x},\tau^{y},\tau^{z}$ are three Pauli matrices. The normal ordering of $\hat{n}^{2}$ and $\hat{n}_{\alpha\beta}^{2}$ are defined as:

[TABLE]

and

[TABLE]

Note that $[\hat{n}_{\alpha\sigma},\hat{n}_{\alpha\sigma^{\prime}}]=0$ and $\hat{n}_{\alpha\sigma}^{2}=\hat{n}_{\alpha\sigma}$, so we have:

[TABLE]

Combining (11) and (12) gives us:

[TABLE]

Using the relation $\bm{\tau}_{\sigma\sigma^{\prime}}\cdot\bm{\tau}_{\zeta\zeta^{\prime}}=2\delta_{\sigma\zeta^{\prime}}\delta_{\sigma^{\prime}\zeta}-\delta_{\sigma\sigma^{\prime}}\delta_{\zeta\zeta^{\prime}}$, we have:

[TABLE]

Note that the $z$ component of magnetic moment operator satisfies:

[TABLE]

Combining (14) and (15) gives us:

[TABLE]

The $\alpha=\beta$ terms in the above sum is :

[TABLE]

so

[TABLE]

And the normal ordering of $\hat{\bm{m}}^{2}$ is:

[TABLE]

Finally, all the normal ordering of local operators have been expressed using local operators in equation (10), (13), (18) and (19).

XI Interacting Hamiltonian for p and d orbitals

The electron-electron interaction for d orbitals can be expressed asCoury et al. (2016):

[TABLE]

where the quadrupole moment operator has been ignored since $\Delta J$ is approximately of magnitude one order lower than that of $J$Oles and Stollhoff (1984). The interaction for p orbitals can be obtained be setting $\Delta J=0$.

$U$ is the Coulomb interaction between $t_{2g}$ orbitals, $J$ is the average exchange splitting of $e_{g}$ and $t_{2g}$ orbitals and $\Delta J$ is the difference of exchange splitting between $e_{g}$ and $t_{2g}$ orbitals. Using the conventional indices (1,2,3,4,5) to represent d orbitals ($3z^{2}-r^{2},zx,yz,xy,x^{2}-y^{2}$), respectively, the parameters can be written as:

[TABLE]

Plug equations (10), (19) and (13) into (20), we have:

[TABLE]

where $\hat{V}_{0}$, spin flipping terms $\hat{V}_{sf}$ and pair hopping terms $\hat{V}_{ph}$ are given by:

[TABLE]

And where site label is omitted for simplicity. The parameters $u,v$ in equation (23) are defined as:

[TABLE]

Now the Hamiltonian of multi-orbital Hubbard model under tight binding approximation is:

[TABLE]

where the hopping integrals are only nonzero between nearest neighbor atomic levels and their magnitudes and signs are discussed in the next section XII. $\hat{V}^{i}$ is the interacting Hamiltonian of each site according to equation (22).

XII The hopping integrals

The $p$ and $d$ orbitals are irreducible representation of spherical symmetry with angular momentum $l=1$ and $l=2$. However, in a local environment with $T_{d}$ symmetry, which is the local symmetry of V atoms and Cu atoms in our case, some of atomic orbitals become reducible, and will split according to the group theory. As shown in table 2, there are 5 irreducible representations of $T_{d}$ group, and $p(l=1)$ and $d(l=2)$ become reducible and split into linear combination of the irreducible representations.

From the character table we get the splitting relations:

[TABLE]

So the $p$ orbitals in $T_{d}$ symmetry only interact with $t_{2g}$ states of $d$ orbitals, and interact only weakly with $e_{g}$ states Jaffe and Zunger (1984). The magnitude and signs of hopping matrix elements for each pair of $p$ orbitals and $t_{2g}$ orbitals can be determined by the Slater Koster matrix Slater and Koster (1954). The magnitude is $t=V_{pd\sigma}/\sqrt{3}$, while the signs are shown in FIG. 11.

XIII Low energy configurations

For anti-parallel case, the lowest energy configuration has energy:

[TABLE]

where $A1$, $C1$ and $A2$, $C2$ are corresponding $e_{g}$ and $t_{2g}$ states. In this configuration, two V(A) atoms have half-filled $3d$ orbitals, two Te(B) atoms have fully filled $5p$ orbitals and Cu(C) atom has empty $t_{2g}$ orbitals and fully filled $e_{g}$ orbitals. The spin configuration of two A atoms in this state is $3d_{\uparrow}^{5}$, $3d_{\downarrow}^{5}$ obeying the Hund’s rule. Electrons on $t_{2g}$ orbitals of A can effectively hop to $t_{2g}$ orbitals of C mediated by $p$ electrons on B. In this process, one B electron hop to C first, and one A electron hop to B next. The energy of intermediate states and the states after one effective hopping are:

[TABLE]

where $\Delta_{AC}=\epsilon_{A1}-\epsilon_{C1}$ and $\Delta_{BC}=\epsilon_{B}-\epsilon_{C1}$ are energy difference between A, C and B $p$, C $e_{g}$ states. And $\delta_{A}$ and $\delta_{C}$ are the crystal field splitting of A and C atoms. Applying the fourth order perturbation theory gives us:

[TABLE]

where the coefficients in front of each terms come from the degeneracy of $p$ and $d$ orbitals.

Similar procedure can be applied to parallel initial spins of V atoms. The lowest energy is:

[TABLE]

We find that it is $24v_{A}$ lower than $E_{\uparrow\downarrow}$. This is because $\psi_{0\uparrow\downarrow}$ state has larger local magnetic moment and results in a smaller energy. Electrons at sites A can hop to site C via a similar effective hopping process. In the parallel case, two A atoms are asymmetric, so we denote the state after one effective hopping of an up(down) spin from $t_{2g}$ orbitals at site A as $\psi_{2\uparrow\uparrow}(\psi_{2\uparrow\uparrow}^{\prime})$. Starting from the lowest energy configuration, we listed the energies of the states with the second and forth order corrections:

[TABLE]

Applying the forth order perturbation theory gives us:

[TABLE]

So the energy difference between parallel and anti-parallel $e_{g}$ spins of V atoms is:

[TABLE]

where $\Delta E_{10}=E_{1\uparrow\downarrow}-E_{0\uparrow\downarrow}=E_{1\uparrow\uparrow}-E_{0\uparrow\uparrow}$, $\Delta E_{20}=E_{2\uparrow\downarrow}-E_{0\uparrow\downarrow}=E_{2\uparrow\uparrow}-E_{0\uparrow\uparrow}$.