Control of Dzyaloshinskii-Moriya interaction in Mn1−xFexGe: a first-principles study

Motivated by the recent experiment on the size and helicity control of skyrmions in Mn1−xFexGe, we study how the Dzyaloshinskii-Moriya (DM) interaction changes its size and sign in metallic helimagnets. By means of first-principles calculations, we successfully reproduce the non-trivial sign change of the DM interaction observed in the experiment. While the DM interaction sensitively depends on the carrier density or the detail of the electronic structure such as the size of the exchange splitting, its behavior can be systematically understood in terms of the distribution of anticrossing points in the band structure. By following this guiding principle, we can even induce gigantic anisotropy in the DM interaction by applying a strain to the system. These results pave the new way for skyrmion crystal engineering in metallic helimagnets.

a fascinating possibility of controlling the DM interaction in metallic systems by manipulating the electronic structure. Indeed, the fact that not only the size but also the helicity of skyrmions in Mn 1−x Fe x Ge changes as a function of x indicates that we have a good chance to control the value of D.
In this paper, we show a quantitative analysis of the DM interaction in the metallic helimagnet, Mn 1−x Fe x Ge, based on ab initio density-functional theory (DFT) calculation. From the obtained band structure, we evaluate the off-diagonal spin susceptibility which is a direct measure of the DM interaction. We find that the sign change of D observed in the experiment for Mn 1−x Fe x Ge is successfully reproduced. The carrier-density dependence of D can be systematically understood in terms of the distribution of band anti-crossing points in the electronic structure. We demonstrate that the sign and the size of D can be controlled as a function of the carrier density or the size of the exchange splitting. There is also an interesting possibility to induce gigantic anisotropy in D by applying a strain to the system.

Results
DM interaction in the continuum model. Let us first look at the second term in the Hamiltonian (1). This indicates that q-linear term in the spin susceptibility, χ αβ , should be proportional to the DM interaction coefficient, D. Therefore, to estimate D in the continuum limit from the DFT calculation, we compute the long-wave length limit of the spin susceptibility, that is, Here, (α, β, γ) = (x, y, z), (y, z, x), or (z, x, y) and D ∼ β corresponds to the coefficient for M α ∂M γ /∂β in Eq. (1). Since we consider the skyrmions in the x-y plane under the total magnetic moment along the z-axis, hereafter we focus on D x ∼ and D y ∼ . In Eq.
(2), we use the non-interacting spin susceptibility defined as where, σ is the Pauli matrix and G 0 is the non-interacting Green's function in the orbital basis. Using this non-interacting spin susceptibility, we can write as where k n is the eigenvector of the Kohn-Sham Hamiltonian with the eigenvalue of ε nk . Hence, we can discuss the DM interaction in terms of the band structure. Although there is a sophisticated approach to compute D 19 , we employ the current simple approach to explore various parameters and materials. Furthermore, this approach is appropriate to obtain a guiding principle for controlling D as discussed below.
Ab initio band structure. Figure 1(a) shows the DFT band structure of FeGe (black solid lines).
Here, we include the spin-orbit couplings and assume the ferromagnetic moment along the z axis. The calculated local magnetic moment is 1.18 μ B per Fe atom, which is consistent with the experiments 22,23 and previous calculations 24 . The x-and y-component of this local magnetic moment is only about 1% of the z-component indicating that the DM interaction is sufficiently small so that Eq. (2) holds in good accuracy. Using this electronic structure, we construct the tight-binding model made of Fe 3d and Ge 4p Wannier orbitals to reproduce the band structure below the Fermi level as shown in red broken lines. The densities of states for up spin (red line) and down spin (blue line) are also shown in Fig. 1(a). As can be seen, there is a large exchange splitting, Δ . According to the energy difference of up and down spins for the Fe 3d orbitals, we obtain Δ = 1.17 eV. In Fig. 1(b), the obtained tight-binding band structure around the Fermi level is illustrated with colors representing the weight of the up spin. Since we consider ferromagnetic electronic structure, each band can be basically characterized as either up-spin or down-spin band as shown in Fig. 1(b). In addition, due to the spin-orbit couplings, there are several anticrossing points where complex spin texture emerges. In the followings, we use this tight-binding model for all the calculations.
Ab initio evaluation of the DM interaction. Let  We assume that the band dispersion is linear in the k x -k y plane, and dispersionless in the k z direction. We introduce m to open a gap at the band crossing point as shown in Fig. 2(a). The static spin susceptibility can be calculated analytically, and the result is  In real materials, the situation is not so simple as that of this two-band model. The anticrossing points form complex surfaces in a four-dimensional space spanned by the energy and the wave number, and the electronic states around neighboring anticrossing points can hybridize with each other. In Fig. 2(c), as a representative case, we schematically show how D ∼ changes as a function of the chemical potential when two anticrossing points with different energies and spin textures with opposite chiralities reside close to each other. Thus, when many anticrossings are densely clustered around the Fermi level, the sign and the size of D ∼ should change drastically. In fact, the band structure of FeGe has many anticrossing points around the Fermi level (see Fig. 1). To visualize the distribution of the anticrossing points as a function of energy, in Fig. 3, we plot a partial density of states (DOS) where the up-spin weight, w ↑ , satisfies 0.05 < w ↑ < 0.95. We can see that there are two peaks around μ = 0 eV and − 0.5 eV, where D ∼ is expected to change significantly.
To analyze the contribution of each anticrossing point, next, we focus on the anticrossing points around the M (0 1/2 1/2) point. In Fig. 4(a), there are two anticrossing points between Y and M points, and between M and Z points. In k x = 0 plane, such anticrossing points continuously form a closed loop around the M point. In Fig. 4 when the chemical potential is below and above the anticrossing points shown in black lines in Fig. 4(a), respectively. We see that the texture around the M point drastically changes when the chemical potential sweeps across the anticrossing points, which is consistent with the simple two-band calculation 20 around the other Fermi surfaces such as the one between Γ -Y line are not negligible although the spin mixture is not so significant. This is because the effect of spin mixing extends away from anticrossing points; that is, k D( ) ∼ has non-negligible values typically up to 0.2-0.3 eV away from anticrossing points. For comparison, in Fig. 4(d,e), we show the Berry curvature, Ω z (k), which is the origin of intrinsic anomalous Hall conductivity (AHC) 21,25 . Ω z (k) is defined as n v n n v n 2 Im where v x and v y are velocity operators. Since the spin mixing is important for both cases, there are common regions where the contributions to the DM interaction and AHC are large. However, in the Berry curvature, the summation in Eq. (7) is restricted to n ≠ n′ while it is not in Eq. (4). As a result, only the restricted region is important for AHC, which is in sharp contrast to D ∼ . To check how large region around the anticrossing points contributes to D ∼ , we calculate D ∼ with restricted summation in Eq. (4). As shown in Fig. 4(f), we find that the contribution from 10%-20% of the total (k, n) points around the band anticrossing points determines up to 70%-80% of the total D ∼ . Figure 5 shows the resulting D x y ∼ , at T = 300 K as a function of the carrier density together with the AHC, σ xy . In the calculation, we use the rigid band approximation. The relation between the carrier density, n, and the chemical potential, μ, is shown in the inset of  hump structure in D μ ( ) ∼ around μ ~ 0.4 (0.2) eV originates from the peak structure around μ ~ 0.5 (0.0) eV in Fig. 3, respectively. If we assume that contributions from μ ~ 0(− 0.5) eV is negative (positive), we can understand why D ∼ changes its sign around μ ~ 0.3 eV. As for σ xy , the calculated value of σ xy in MnGe is larger than that in FeGe and there is no sign change. This behavior including size and sign agrees well with the experimental anomalous Hall contribution to σ xy for MnGe 26 and for Mn 1−x Fe x Ge 27 .

Strain-induced huge anisotropy of the DM interaction.
According to the above discussion, there are several ways to change the size and sign of D. As we have seen above, D can be efficiently controlled by carrier doping. We can also exploit the temperature dependence of the exchange splitting. When the relative position of up-and down-spin band changes, the distribution of anticrossing points in the band structure also changes, which will have a direct impact on D. As an example, we will show later the change in D x y ∼ , of FeGe as a function of the moment per Fe atom. This mechanism can be related to the temperature dependence of the magnetic moment and skyrmion size in MnGe 26 . Another interesting possibility is to make use of the strain effect. If we apply a strain to the system, the symmetry of the electronic structure can be lowered, and the distribution of the anticrossing points will change drastically. This effect is expected to be prominent especially when D changes its sign. Figure 6(a) shows the calculated D x y ∼ , at n = − 0.45 for which we apply the uniaxial strain along the y direction. We find that the difference between D x ∼ and D y ∼ actually enhances particularly by the elongation along the y axis; D x ∼ is about 40% (400%) larger than D y ∼ at + 2% (+ 5%) strain. The carrier density dependence of the anisotropy, D D y x / ∼ ∼ , for fixed strain of + 5% is shown in Fig. 6(b). We can see that the anisotropy becomes large particularly around the region of sign change.

Discussion
In the present calculation, we employed the rigid band approximation. To examine its validity, we have performed a calculation for the other end material MnGe and doped negative carriers by the rigid band approximation. As shown in Fig. 7, we have obtained a qualitatively similar result in that D ∼ is positive (negative) for the end material MnGe (FeGe). Thus the result that D ∼ for Mn 1−x Fe x Ge changes its sign between x = 0 and 1 should be robust, even when we go beyond the rigid band approximation.
Regarding the crystal structure, it has been known that the magnetic moment for the optimized structure is much smaller than the experimental value within the local density approximation 28 . On the other hand, for the experimental structure, the size of magnetic moment is similar to that in the experiment (~1μ B ). In fact, within our calculations, the magnetic moment in FeGe decreases with decreasing the lattice constant as shown in Fig. 8(a), which is consistent with previous studies 24,29 . As a result, D ∼ exhibits non-trivial magnetic moment dependence as shown in Fig. 8(b). Since the magnetic moment dependence of D is also an important issue, let us next discuss how the change in the lattice structure or the magnetic moment affect D in Mn 1−x Fe x Ge. In Fig. 8(c), we plot the lattice constant dependence of D ∼ , together with the distribution of band anti-crossings in Fig. 8(d,e) for lattice constants a = 0.99a 0 and 1.01a 0 , where a 0 is the experimental lattice constant. As can be seen, for larger a and magnetic moment, the energy difference of two peaks in Fig. 3 becomes larger (Fig. 8(e)). Consequently, the density at which D ∼ changes its sign (n c ) becomes larger (Fig. 8(c)). However, the qualitative feature of D ∼ is robust against the change in a. Therefore, our calculation successfully explains why D ∼ is positive (negative) for MnGe  Table I).

DFT calculations.
To evaluate 0 χ αβ in FeGe, we perform the electronic structure calculation within the generalized-gradient approximation (GGA) 32 based on the density functional theory 33 . We use ultrasoft pseudopotentials 34 and plain-wave basis sets to describe the charge densities and wave functions with cutoff energies of 40Ry and 500Ry, respectively. We use 8 × 8 × 8 k-point mesh. With including the spin-orbit couplings and assuming the ferromagnetic moment along the z axis, we obtain non-collinear magnetic structure. Using this electronic structure, we calculate Wannier functions for Fe 3d and Ge 4p    orbitals using wannier90 code [35][36][37] . Based on the Wannier functions, we construct a tight-binding model on the restricted Hilbert space and calculate χ αβ using 64 × 64 × 64 k-point mesh at T = 300 K. The AHC is also calculated using the Wannier interpolation technique with 200 × 200 × 200 k-point mesh 38 .