Theory of defect-induced crystal field perturbations in rare earth magnets

We present a theory describing the single-ion anisotropy of rare earth (RE) magnets in the presence of point defects. Taking the RE-lean 1:12 magnet class as a prototype, we use first-principles calculations to show how the introduction of Ti substitutions into SmFe$_{12}$ perturbs the crystal field, generating new coefficients due to the lower symmetry of the RE environment. We then demonstrate that these perturbations can be described extremely efficiently using a screened point charge model. We provide analytical expressions for the anisotropy energy which can be straightforwardly implemented in atomistic spin dynamics simulations, meaning that such simulations can be carried out for an arbitrary arrangement of point defects. The significant crystal field perturbations calculated here demonstrate that a sample which is single-phase from a structural point of view can nonetheless have a dramatically varying anisotropy profile at the atomistic level if there is compositional disorder, which may influence localized magnetic objects like domain walls or skyrmions.

The unique magnetic properties of lanthanide elements originate from their partially-filled 4f shells [1].Their aspherical charge distribution can yield a highly anisotropic dependence of the total energy on the magnetic moment direction, i.e. a gigantic magnetocrystalline anisotropy (MCA), which has been exploited with huge commercial success in rare earthtransition metal (RE-TM) permanent magnets [2].RE atoms also lie at the center of recent fundamental research into single molecule magnets [3], skyrmionics [4,5] and other forms of topological magnetism [6].A quantitative description of RE MCA is provided by crystal field (CF) theory [7], where the environment of each RE atom is described using a set of CF coefficients.These are inserted via Stevens operators [8] into quantum or classical Hamiltonians acting on the 4f electrons [9], from which one can extract equilibrium properties through statistical mechanics [10], or carry out time-evolving simulations within the framework of atomistic spin dynamics (ASD) [11].Further coarse-graining the atomistic results yields expressions for the MCA energy for use in continuum micromagnetics [12].
A key target in permanent magnet development is a high maximum energy product (BH) max , which in turn requires a high coercivity (resistance to demagnetization by external fields).Although a large MCA is a prerequisite for a large, shape-independent coercivity [13], RE magnets typically only achieve 5-35% of their theoretical limit [14].According to an influential analysis by Kronmüller [15], this performance gap is due to inhomogeneities in the sample (e.g.misaligned or nonmagnetic grains) creating magnetically soft regions susceptible to the nucleation of domains with reversed magnetization.These grow through domain wall motion and eventually lead to destruction of the original magnetic order.
Kronmüller's analysis is a top-down approach, which posits that an inhomogeneity induces a variation in the anisotropy according to an assumed analytical profile.However, computational magnetics is now focused on bottom-up simulations, whereby magnetization reversal and domain wall motion are treated as phenomena emerging from atomistic properties [16].ASD simulations of homogeneous materials are able to reproduce intrinsic quantities in impressive agreement with experiment [17,18].However, to better understand the coercivity gap, it is essential to develop a method of incorporating inhomogeneities into these simulations [19].
As an example, we consider how an ASD simulation based on a classical spin Hamiltonian should be modified to include a point defect, the smallest possible inhomogeneity.In ASD, each magnetic atom contributes to the total energy based on the orientation of its magnetic moment ê = (sin Θ cos ϕ, sin Θ sin ϕ, cos Θ) [11].If a point defect is placed at a position R J with respect to an RE atom, we derive the following additional contribution to the energy from the perturbed crystal field: Both d 0m l (Θ) and A l are known quantities: d 0m l (Θ) is a function related to the associated Legendre polynomials [20,21], and A l is an RE-dependent prefactor formed from Stevens coefficients and the total angular momentum J [8,22].l and m are the usual angular quantum numbers.Therefore, the only unknowns in equation 1 are the ∆B lm quantities, describing the defect-induced change in each CF coefficient.It is the calculation of these quantities, which are also required in a quantum treatment of the RE magnetism, which is the subject of this Letter.Specifically, we initially obtain ∆B lm using first-principles density-functional theory (DFT), then use the calculations to parameterize a highly efficient model suitable for deployment in large-scale ASD simulations.
Point defects, and their effects on the MCA, are particularly important in the "1:12" class of RE-TM alloys [23].These compounds are under intense investigation as "rare earth-lean" magnets, with their RE:TM atom ratio of 1:12 being more economical than 1:7 or 1:8.5 as is not stable in bulk, so it is necessary to introduce a transition metal M, like Ti, to form REFe 12-x M x , where x ∼ 1 [24].Although the M atoms prefer to sit at particular sites in the unit cell, the substitution is probabilistic, so they can be viewed as substitutional point defects embedded in an REFe 12 host.Despite promising intrinsic properties, sufficiently high coercivities have not yet been achieved, leading naturally to the question of how this intrinsic disorder affects 1:12 performance.
We now describe our approach.In CF theory, the potential V (r) experienced by the RE-4f electrons is expanded in terms of atom-centred functions, which here are the complex spherical harmonics Y lm , i.e.V (r) = lm V lm (r)Y lm (r).The radial functions V lm (r) are then combined with the 4f electron charge density n 4f (r) and integrated, to construct the CF coefficients B lm [42].The anisotropy energy of the 4f electrons is affected only by CF coefficients with l equal to 2, 4 or 6.Furthermore, CF coefficients will be nonzero only if they are compatible with the symmetry of the RE site.SmFe Figure 1 shows the effect on the crystal field at the two different sites, plotting the change in each coefficient ∆B lm .We see that, as well as perturbing the existing nonzero CF coefficients like B 20 , the Ti substitution generates new coefficients like B 22 , due to the lower symmetry.
The perturbations show two general features: first, that the largest changes occur at the 2a site closest to the substituted Ti, and second, that the largest changes are in the l = 2 coefficients.The same behavior occurs when the Ti is substituted at an 8j or 8f site [21].
The strong l-and interatomic distance dependence of the CF perturbation brings to mind the Laplace expansion of the Coulomb potential and suggests that the Ti substitution might be efficiently modeled as a screened point charge [7].We use this idea to derive an expression (in Hartree atomic units) for the change in CF coefficient when a dopant is introduced at site J: We have previously calculated ⟨r l ⟩ = r l+2 n 4f (r)dr with self-interaction-corrected DFT [42,47] and found it to be largely independent of the RE environment; for Sm, ⟨r l ⟩ is 1.02, 2.46 and 10.25 for l = 2, 4 and 6 [21].∆Z J is the point charge associated with the dopant, which depends on whether site J is 8i, 8j or 8f .To account for the short-ranged behavior shown in Fig. 1, we impose a cutoff radius R c such that only defects with R J < R c contribute to ∆B lm .Therefore, our model has four parameters, which we determine from a least-squares fit to the DFT calculations of ∆B lm , (Fig. 2): we find ∆Z 8i = 0.238, ∆Z 8j = 0.251, ∆Z 8f = 0.046 and R c = 4 Å.The model captures the most significant changes in CF coefficients, with errors in smaller ∆B lm values of order 50 K.
We note that the use of analytical models is hardly a new idea in CF theory [7], with Li and Cadogan's bonding charge model applied specifically to 1:12 compounds [48].However, our current approach carries two particular advantages.First, we are parameterizing the model using direct DFT calculations of CF coefficients, rather than having to extract them indirectly from experimental measurements [25,49].Second, we are only using the analytical model to describe the perturbation to the crystal field.This both reduces the number of parameters in the model and also increases the number of data points that can be used in the fitting process, since the lower symmetry produces more nonzero CF coefficients.
Having established the model, we now consider its predictions.Since equation 2 is an electrostatic model, the principle of superposition applies, such that introducing multiple substitutions has an additive effect on the CF coefficients.At SmFe 11 Ti stoichiometry, two of the 24 TM sites in the conventional cell are occupied by Ti atoms, on average.It is a simple matter to generate the 276 possible configurations (many of which are not unique), and use equation 2 to evaluate the perturbation to the CF.The predicted effect on the (2, 0) coefficient is shown in Fig. 3 as blue crosses.Here, the CF is quantified using the conventional A 20 ⟨r 2 ⟩ coefficient, which is equal to half of B 20 [7,42].The largest changes are predicted to occur when both Ti atoms are located either on 8i or on 8j sites, close (within R c ) to the RE atom.Furthermore, due to opposing signs of ∆B 20 for 8i or 8j, there is a cancellation effect for mixed 8i-8j substitution.
To verify these predictions, we took a subset of the 276 configurations and calculated the CF coefficients using DFT.The resulting ∆A 20 ⟨r 2 ⟩ values are shown as stars in Fig. 3, and emphasize the excellent qualitative and quantitative agreement between the two methods.
Similarly good agreement is found for other coefficients, e.g.∆A 22 ⟨r 2 ⟩ [21].We stress the independence of the training set in Fig. 2 and the test set in Fig. 3, i.e. the data in Fig. 3 was not used to refine the parameters.for SmFe 12 from a detailed dynamical mean field theory treatment of on-site correlation (-198 K) [38] with our calculated off-site electrostatic effect of Ti doping, averaged over configurations (+72 K).
According to our calculations, if two Ti atoms occupy 8i sites close to Sm, the perturbation to the CF will be so large that it will almost destroy the single-ion anisotropy.
Therefore, it is crucial to know how the dopants are distributed.We calculated the energies of all possible arrangements of Ti atoms within the conventional cell of the analogous material YFe 11 Ti, and like previous calculations and experiments on 1:12 systems [27,29,52,53], found that the Ti atoms overwhelmingly prefer 8i sites.Importantly, we also found the most stable configuration has one Ti atom per (001) plane, which corresponds to both Sm sites having the same +83 K shift in A 20 ⟨r 2 ⟩ shown in Fig. 3.However, there are configurations which are higher in energy by only 6 meV/FU where alternate (001) planes contain two or zero Ti atoms respectively [21].Such a structure would have an intrinsically textured anisotropy, where Sm sites in adjacent (001) planes have alternating A 20 ⟨r 2 ⟩ shifts of +155 and -23 K. Given the closeness in energy of these configurations, it is highly likely that synthesized SmFe 11 Ti samples will have Sm atoms with a distribution of anisotropies, with A 20 ⟨r 2 ⟩ values ranging from ∼ -200 to -50 K. Furthermore, any Ti-rich mesoscopic regions would have a significantly weakened (possibly even planar) anisotropy.
Although we have focused on the (2, 0) coefficient, Fig. 1 clearly shows how the (2, ±2) coefficient is also substantial.It is important to realize that, assuming that the Ti dopants are uniformly distributed across a set of Wyckoff sites, the value of ∆B lm averaged over all Sm sites will be zero unless m = 0 or ±4.Physically, this averaging process is equivalent to treating all of the Sm atoms as a single entity-a magnetic sublattice-with a single, common magnetization direction, and is valid when comparing to measurements of macroscopic phenomena of bulk crystals, e.g.measurements of magnetization like in Ref. [25], which see an averaged value of ∼ −100 K for A 20 ⟨r 2 ⟩.However, any noncollinear and/or localized magnetic texture, such as a domain wall or a skyrmion [4,5], will probe the anisotropy of only a subset of the RE atoms, leading to an incomplete averaging.Therefore, in simulating these phenomena, it is essential to consider these additional CF coefficients, even if they average to zero in the macroscopic limit.We believe that there may be an interesting link between our work and experiments carried out in the 1970s investigating an "intrinsic coercive force' in ternary RETM 5 alloys like SmCo 5-x Ni x [54,55].Intriguingly, it was found that even in samples containing only a single phase, there existed a coercive field whose magnitude depended on the concentration of the substituted element.It was pointed out that the Bloch walls separating domains should only be a few atomic layers thick [56,57], and it was argued that the dopants affected these through modification of the exchange interaction.However, the CF perturbation due to dopants was not considered, and it is important to revisit this phenomenon using our new approach.
We have focused our discussion on magnetocrystalline anisotropy, but it is important also to quantify the effect of Ti substitutions on the RE-TM exchange.Using the DFT formulation of the disordered local moment picture (DFT-DLM [58][59][60][61][62]), we calculated that the exchange fields at either of the two Sm sites in SmFe 11.5 Ti 0.5 are within 3 or 7% of those in SmFe 12 , at 300 K. We found a similarly small variation for SmFe 11 Ti, even with the Ti defects clustered around only one Sm site [21].Based on these calculations, we conclude that the Ti dopants do not appreciably modify the RE-TM coupling to Sm.
We have applied the point charge model to another atomistic defect: the RECo 5 "dumbbell" whereby an RE atom is substituted with a Co 2 pair [63].Carrying out the same procedure as for SmFe 12-x Ti x , we find a negative value of ∆Z Co 2 , -0.465, and that the only RE sites affected by the presence of the Co 2 pair are the nearest neighbors lying directly above or below it.The negative value of ∆Z Co 2 has a detrimental effect on the anisotropy; pristine SmCo 5 has very strong uniaxial anisotropy, with the prolate electron cloud pointing along the c-axis [44].Introducing the negative point charge along the same direction reduces the magnitude of A 20 ⟨r 2 ⟩.This microscopic behavior is reflected in the bulk properties of Sm 2 Co 17 , which can be viewed as an ordered array of Co 2 defects implanted in SmCo 5 [63], and has a substantially reduced A 20 ⟨r 2 ⟩ (50%) compared to SmCo 5 [42].
As a summary, we set out our proposed scheme to carry out ASD simulations of localized magnetic phenomena in the presence of point defects.First, the Hamiltonian of the host material is constructed, taking advantage of the widening literature of CF coefficients and exchange constants for pristine systems.Next, distributions of point defects are generated, which may be totally disordered, partially or fully ordered; the defect distributions might be built based on atom probe tomography data, particularly when such measurements correlate with substantial variations in coercivity [64].Equation 2 is then used to calculate the CF perturbations at each RE atom, which are then fed into equation 1 to produce the Hamiltonian; an example of this step is provided as Supplemental Information [21].Our DLM calculations indicate it is not necessary to modify the exchange part of the Hamiltonian, and we also suggest that only the l = 2 terms need to be included in the CF perturbation.
Ideally, the ∆Z parameters should be obtained by fitting to DFT calculations, but we stress their physical meaning as being the effective charge of the defect.Therefore, an order of magnitude estimate for ∆Z could be obtained from a relatively simple DFT calculation on the defective system, obtaining the charge e.g. through a Bader analysis [65].From here, all of the sophisticated techniques developed for perfect systems [17,18] may be applied.The framework presented here thus opens a new avenue in computational magnetics research, enabling the atomistic study of domain wall propagation and other localized magnetic phenomena in systems with inhomogeneous anisotropy.
Computational Details: CF coefficients were calculated within DFT using the Y-analogue formalism as described in Ref. [42], using the GPAW code [66].A plane-wave basis with cutoff energy of 800 eV and 6×6×6 reciprocal space sampling were used, and exchange and correlation effects were treated within the local-spin-density approximation (LSDA) [67].
The structural parameters of the SmFe 12 cell (a, c = 8.497, 4.687 Å, x 8i = 0.359, x 8j = 0.270) were taken from previous DFT (generalized-gradient approximation) calculations [50] and were used throughout the study, except for the total energy comparison and indicated CF calculations, where the internal co-ordinates were allowed to relax.The DFT CF coefficients were averaged over spin directions before fitting to the point charge model.

FIG. 1 .
FIG. 1. Calculated perturbations to the CF coefficients when a single Ti atom (blue) is introduced at an 8i site.The Sm atoms at the 2a sites are shown in purple and pink.See Fig. 2 for identification of the individual 8i, 8j, 8f sites.All CF coefficients are real, and ∆B lm = ∆B l−m .
12-x Ti x crystallizes in the ThMn 12 -type structure shown in the inset of Fig 1.The conventional unit cell is tetragonal and contains two formula units, with Sm atoms at the corner and centre positions (Wyckoff label 2a) and transition metals occupying three distinct sublattices, 8i, 8j and 8f .In the pristine SmFe 12 host, these sublattices are all occupied by Fe atoms.This results in the 2a sites having D 4h symmetry, whereby the RE-4f crystal field is specified completely by the (l, m) coefficients (2, 0), (4, 0), (4, ±4), (6, 0) and (6, ±4)[9,43].We have calculated these within DFT using the Y-analogue model[42], which we have previously developed and applied to RE-Co, Nd-Fe-B and REFe 2 -type magnets[44][45][46].Computational Details are provided at the end of this Letter.Now we substitute one Fe atom with Ti at an 8i site, yielding SmFe 11.5 Ti 0.5 .As shown in the inset of Fig.1, this substitution breaks the equivalence of the 2a sites, with the Sm atoms at the cell corners lying further away from the Ti compared to the Sm atom located at the centre position.Both 2a sites have their symmetry lowered to C 2v .

FIG. 2 .
FIG. 2. Perturbation to CF coefficients ∆B lm calculated within DFT or using equation 2, when one Ti atom is placed in the conventional cell at an 8i, 8j or 8f site.All (l, m) coefficients are shown for both 2a sites.+ symbols are imaginary parts.The inset shows the different sites in pristine SmFe 12 .

FIG. 3 .
FIG. 3. Perturbation to the A 20 ⟨r 2 ⟩ CF coefficient when two Ti defects are introduced into the conventional unit cell at the indicated crystal sites, demonstrating excellent agreement between point charge (blue cross) or DFT (stars) calculations.

Figure 3 shows how the A 20 ⟨r 2 ⟩
Figure 3 shows how the A 20 ⟨r 2 ⟩ CF coefficient is affected by several hundreds of K when one or two Ti atoms are substituted close to the RE site.The three DFT-calculated values of ∆A 20 ⟨r 2 ⟩ shown in Fig. 3 for two 8i substitutions are +155, +83 and -23 K, corresponding to two, one and zero Ti atoms lying close to the RE site, respectively.To place these numbers into context, calculations of A 20 ⟨r 2 ⟩ for the REFe 12 host find values between -50 and -200 K [38, 50, 51], with our own calculations giving -88 K or -136 K depending on whether experimental or DFT-optimized internal co-ordinates are used.Meanwhile, the latest experimental analysis of magnetization curves deduced a value of A 20 ⟨r 2 ⟩ of -110 K for SmFe 11 Ti [25].This agrees very well with the value obtained by combining A 20 ⟨r 2 ⟩ We thank A. Gabay and G. Hadjipanayis for useful discussions.This work was supported in part by a Royal Society Research Grant No. RGS\R1\201151 and by the U.S. Department of Energy, Office of Basic Energy Sciences under Award Number DE SC0022168.Y. H. acknowledges funding from the Ironmongers' Foundation.We acknowledge the University of Oxford Advanced Research Computing (ARC) facility in carrying out this work [68].Figures were rendered using VMD [69].