# Crystal-field splittings in rare-earth-based hard magnets: An ab initio
approach

Pascal Delange  Centre de Physique Théorique, École Polytechnique, CNRS,
Université Paris-Saclay, 91128 Palaiseau Cedex, France  Silke Biermann  Centre
de Physique Théorique, École Polytechnique, CNRS, Université Paris-Saclay,
91128 Palaiseau Cedex, France  Collège de France, 11 place Marcelin Berthelot,
75005 Paris, France  Takashi Miyake  CD-FMat, AIST, Tsukuba 305-8568, Japan
Leonid Pourovskii  Centre de Physique Théorique, École Polytechnique, CNRS,
Université Paris-Saclay, 91128 Palaiseau Cedex, France  Collège de France, 11
place Marcelin Berthelot, 75005 Paris, France  Materials Modeling and
Development Laboratory, National University of Science and Technology “MISIS”,
119049 Moscow, Russia

###### Abstract

We apply the first-principles density functional theory + dynamical mean field
theory framework to evaluate the crystal-field splitting on rare-earth sites
in hard magnetic intermetallics. An atomic (Hubbard-I) approximation is
employed for local correlations on the rare-earth 4$f$ shell and self-
consistency in the charge density is implemented. We reduce the density
functional theory self-interaction contribution to the crystal-field splitting
by properly averaging the 4$f$ charge density before recalculating the one-
electron Kohn-Sham potential. Our approach is shown to reproduce the
experimental crystal-field splitting in the prototypical rare-earth hard
magnet SmCo5. Applying it to $R$Fe12 and $R$Fe${}_{12}X$ hard magnets ($R=$Nd,
Sm and $X=$N, Li), we obtain in particular a large positive value of the
crystal-field parameter $A_{2}^{0}\langle r^{2}\rangle$ in NdFe12N resulting
in a strong out-of-plane anisotropy observed experimentally. The sign of
$A_{2}^{0}\langle r^{2}\rangle$ is predicted to be reversed by substituting N
with Li, leading to a strong out-of-plane anisotropy in SmFe12Li. We discuss
the origin of this strong impact of N and Li interstitials on the crystal-
field splitting on rare-earth sites.

##  I Introduction

Permanent magnets are a key component of modern electronic devices, ranging
from electric motors to medical imaging. An important breakthrough in the
quest for high-performance permanent magnets was the discovery of rare-earth
intermetallic magnets, starting with SmCo5 in 1966Strnat (1970). Since its
discovery in 1982, the champion of hard magnets has been Nd2Fe14BHerbst
(1991). More recently, rare-earth iron-based hard magnets $R$Fe${}_{12}X$ with
the ThMn12 structure such as NdFe12N have been under renewed scrutinyKörner et
al. (2016); Suzuki et al. (2016); Harashima et al. (2015a); Hirayama et al.
(2015); Miyake et al. (2014). The underlying reason is the high price and
strategical importance of rare-earths and cobalt, and the ongoing research
effort to find good permanent magnets with reduced rare-earth
concentrationCoey (2011). The ThMn12 structure has a reduced ratio of rare-
earth vs transition metal compared to Nd2Fe14B, but nevertheless conserves
strong hard magnetic properties (large magnetization and Curie temperature,
and strong anisotropy) when doped with light elements such as nitrogenSuzuki
et al. (2016); Harashima et al. (2015a); Hirayama et al. (2015).

The main physical ingredients for a rare-earth hard magnet are the high
magnetic anisotropy energy provided by rare-earth ions combined with the high
magnetization and Curie temperature from the transition metal sublattice,
typically composed of Fe or Co atomsCoey (1996); Gutfleisch et al. (2011);
Coey (2011). The 3$d$ transition metal atoms carry little anisotropy; because
of their rather small spin-orbit coupling, their magnetization direction is
essentially fixed by that of the rare-earth ion through an exchange coupling.
The majority of rare-earth elements, especially heavy rare-earth elements, are
very expensive. Moreover, the magnetic moment of heavy rare-earth is normally
anti-parallel to the transition-metal one reducing the net magnetizationCoey
(2011). Hence, one advantage of new compounds like $R$Fe${}_{12}X$ is a
reduced rare-earth concentration. In turn, a higher Fe concentration is
favorable for achieving a large magnetization, which is another advantage of
$R$Fe${}_{12}X$ compounds. However, this reduced rare-earth concentration
means each rare-earth ion must contribute a strong magnetic anisotropy to keep
the overall magnetic hardness. The preferred magnetization direction (in-plane
or out-of-plane) of a given rare-earth ion is determined by the interplay
between the crystal-field (CF) splitting and spin-orbit (SO) interaction. To
the first order, the crystalline magnetic anisotropy energy reads:

$E_{A}\approx K_{1}\textrm{sin}^{2}\theta$

where $\theta$ is the angle between the magnetization and the easy axis, and

$K_{1}=-3J(J-\frac{1}{2})\alpha_{J}A_{2}^{0}\langle r^{2}\rangle n_{R}$ (1)

where $J$ is the total angular momentum for the rare-earth 4$f$ shell, $n_{R}$
the concentration of rare-earth atoms, $\alpha_{J}$ the corresponding Stevens
factor and $A_{2}^{0}\langle r^{2}\rangle$ is the lowest-order crystal-field
parameter (CFP). Additional small doping of light elements is also found to
strongly modify the anisotropy by affecting the rare-earth CF splitting
Matsumoto et al. (2016); Harashima et al. (2016); Harashima et al. (2015a).
They also modify the structural stability: doping B makes the Nd2Fe14B phase
more stable, while interstitial nitrogen has only a minor effect in structural
stability.

It follows that the CF splitting on the rare-earth 4$f$ shell is a crucial
quantity defining the magnetic hardness of rare-earth intermetallics. The
theoretical search for new rare-earth hard magnets thus requires a reliable
approach to calculating CFP. The importance of crystal-field effects for the
optical, magnetic, and other properties of solids has been recognized long
ago, and semi-empirical models of the CF Hamiltonian, such as the point charge
modelHutchings (1964) and the superposition modelNewman and Ng (1989), have
been developed since the 60’s. While they provide an inexpensive and
physically transparent description of CF parameters, their predictive power is
limited as they require experimental input to determine the actual values.
Experimental information is readily available for large band-gap rare-earth
insulators, where the CFP can be extracted from measurements of dipole-
forbidden optical transitions between $f$-statesBrecher et al. (1967). In the
case of rare-earth intermetallics, where the $f-f$ transitions are hidden by
the optical response of conduction electrons, inelastic neutron spectroscopy
can be used to determine CFPBuschow et al. (1974); Sankar et al. (1975);
Givord et al. (1979); Tie-Song et al. (1991); Tils et al. (1999), but its
results are more ambiguous as one needs to sort out the contributions of
phonons and the effect of inter-site exchange interactions.

Ab initio calculations do not rely on experimental input and can have truly
predictive power. First-principles techniques for computing the CF
parametersCoehoorn (1991); Daalderop et al. (1992); Steinbeck et al. (1994);
Novák and Kuriplach (1994); Steinbeck et al. (1996); Novák (1996); Ning and
Brivio (2007); Hu et al. (2011); Novák et al. (2013); Miyake et al. (2014) can
be separated into two main approaches. The first one Daalderop et al. (1992);
Steinbeck et al. (1994); Novák and Kuriplach (1994); Steinbeck et al. (1996);
Novák (1996); Ning and Brivio (2007); Miyake et al. (2014) consists in
extracting the nonspherical Kohn-Sham potential $V_{lm}$ and the 4$f$ charge
density $\rho_{4f}$ around the rare-earth site and then computing the
corresponding crystal-field parameter. As the density functional theory (DFT)
is not able to fully capture the physics of partially-filled localized 4$f$
shells, one imposes their localization by treating the 4$f$ orbitals as semi-
core states. The non-spherical 4$f$ charge density $\rho_{4f}(\boldsymbol{r})$
of the rare-earth ion includes an unphysical contribution to the CFP stemming
from the local-density-approximation (LDA) self-interaction error. This is
usually corrected by spherically-averaging the 4$f$ charge density, but then
approximations have to be made for the long-range “tails” of
$\rho_{4f}(\boldsymbol{r})$.

The importance of excluding the self-interaction of the non-spherical part of
the partial 4$f$ charge density to obtain proper crystal-field energies was
first recognized by Brooks _et al._ in a publication aimed at calculating the
spin Hamiltonian parameters of rare-earth compoundsBrooks et al. (1997).

In the second, more recent, approach the 4$f$ states are represented by
Wannier functionsHu et al. (2011); Haverkort et al. (2012); Novák et al.
(2013), while the charge density and, correspondingly, the Kohn-Sham potential
are generated by self-consistent DFT calculations with 4$f$ states treated as
semi-core. An additional ad hoc parameter is used to correct the charge
transfer energy between 4$f$ and conduction bands.

One may also mention recent work on determining the CFP of lanthanides and
transition metals using quantum-chemical methods, in particular, in order to
understand the properties of magnetic molecules. Such approaches employ, for
instance, the complete active space self-consistant field methodLiddle and van
Slageren (2015) or multireference second order perturbation theoryUngur and
Chibotaru (2017). Here, however, we choose to focus on perfect crystals rather
than on molecules.

Overall, ab initio calculations of CFP for rare-earth ions are a formidable
theoretical problem, due to generally small values of those CFP and their
extreme sensitivity to computational details. The main weak point of
previously proposed DFT-based approaches is that they are not able to
correctly treat the localized valence 4$f$ states. Hence, the charge density
is derived under the drastic approximation of treating them as fully localized
core states, spherically-averaged inside the atomic sphere. The DFT+U method
provides a more realistic treatment for the 4$f$ density in the limit of
strong ordered magnetism. However, it is usually not able to capture the true
quasi-atomic (multiplet) nature of rare-earth shells in the paramagnetic or
partially-polarized state. The DFT+U calculations can nevertheless be used to
estimate the CFP by converging them to the on-site density matrix
corresponding to a given atomic wave function. The CF splitting can be then
evaluated from the difference in DFT+U total energy between such calculations
for relevant CF states. This method in fact makes use of the (usually
inconvenient) tendency of DFT+U to remain in a local energy minimum instead of
converging to the ground state density. Zhou _et al_ Zhou et al. (2011)
employed this approach together with an orbital-dependent self-interaction
correctionZhou et al. (2009), to obtain total energies for different orbital
occupancies in UO2 and deduce its CFP.

In this work, we propose an approach to ab initio CFP calculations based on
self-consistent DFT+Dynamical Mean-Field Theory (DFT+DMFT)Anisimov et al.
(1997); Lichtenstein and Katsnelson (1998) treating the local many-body
problem for the 4$f$ shell in the quasi-atomic (Hubbard-I) approximation.
While this approach of using DFT+DMFT with the Hubbard I approximation, which
we may call DFT+Hub-I, is rather simple and computationally efficient, it was
shown to capture not only the 4$f$ multiplet structure in the paramagnetic
stateLichtenstein and Katsnelson (1998); Lebègue et al. (2005, 2006);
Pourovskii et al. (2009); Locht et al. (2016) and in the ferromagnetic
stateGrånäs et al. (2012), but also the 4$f$–conduction band exchange
interaction and the resulting exchange splitting of the Fermi-
surfacePourovskii et al. (2009). This scheme also provides a rather natural
way of averaging the 4$f$ partial density to reduce the self-interaction error
from the CF Hamiltonian. We validate it by applying it to the well-known hard
magnet SmCo5, for which the crystal-field splitting has been measured in
multiple experimentsBuschow et al. (1974); Sankar et al. (1975); Givord et al.
(1979); Tie-Song et al. (1991); Tils et al. (1999). We then apply our method
to much less investigated new hard magnets of the $R$Fe${}_{12}X$ family,
computing their CFP for different rare-earth elements (Sm or Nd) and
considering N and Li interstitials. Our calculations predict the hypothetical
SmFe12Li compound to possess a strong axial anisotropy and, possibly,
interesting hard-magnetic properties.

The paper is organized as follows: in Sec. II.1 we introduce basic notions as
well as relevant notations of the CF theory. Our first-principles
computational approach is presented in more details in Sec. II.2. Our results
for the DFT+Hub-I electronic structure and CFP for the $R$Fe${}_{12}(X)$ hard
magnets are presented in Secs. III.1 and III.2, respectively. In Sec. IV we
analyze the shape 4$f$ Wannier functions (WF) in real space and employ a
projective approach to evaluate the WF localization and the contribution of
hybridization effects to CFP.

##  II Method

###  II.1 Crystal-field parameters: notation and symmetry

We start by introducing crystal-field parameter notations. The local
Hamiltonian for a rare-earth ion with a partially-filled 4$f$ shell subject to
the exchange field created by the transition-metal sublattice and to a
crystal-field potential reads

$\hat{H}=\hat{H}_{\mathrm{1el}}+\hat{H}_{U}=\hat{E}_{0}+\lambda\sum_{i}s_{i}l_{i}+2\mu_{B}B_{\mathrm{ex}}\hat{S}_{a}+\hat{H}_{cf}+\hat{H}_{U}$
(2)

where the one-electron part of the Hamiltonian corresponds to the first four
terms on the right-hand side, namely, a uniform shift, spin-orbit, exchange-
field, and crystal-field terms. $\hat{S}_{a}$ is the in-plane or ouf-of-plane
spin operator, corresponding to the case where $B_{\mathrm{ex}}$ is along $x$
or $z$, respectively. $\hat{H}_{U}$ represents the electron-electron Coulomb
repulsion term of the many-body Hamiltonian. The crystal-field term
$\hat{H}_{cf}$ is defined as the non-spherically symmetric part of the one-
electron Hamiltonian. The corresponding non-spherical part
$V_{ns}(\boldsymbol{r})$ of the one-electron potential can be expanded into
spherical harmonics as follows

$V_{ns}(\boldsymbol{r})=\sum_{k=1}^{\infty}\sum_{q=-k}^{k}A_{k}^{q}(r)Y_{kq}(\hat{\boldsymbol{r}}).$
(3)

where $Y_{kq}(\hat{\boldsymbol{r}})$ is the spherical harmonic function with
total angular moment $k$ and projected angular moment $q$. The matrix elements
of $V_{ns}(\boldsymbol{r})$ between 4$f$ orbitals define $\hat{H}_{cf}$. Due
to the properties of the spherical harmonics, only $A_{k}^{q}$ for $k\leq 2l$,
i.e. $k\leq 6$ in the case of an $f$ shell, can contribute to $\hat{H}_{cf}$ .
In the point-charge CF calculations $A_{k}^{q}(r)$ is reduced to
$A_{k}^{q}r^{k}$. While we do not assume this form for $A_{k}^{q}(r)$ in the
present formalism we still employ the now standard notation $\langle
A_{k}^{q}(r)\rangle\equiv A_{k}^{q}\langle r^{k}\rangle$. For historic
reasons, several conventions exist for the parametrization of $\hat{H}_{cf}$,
leading to a rather confusing variety of definitions for the crystal-field
parameters. Using the so-called Stevens operator equivalentsStevens (1952),
$\hat{H}_{cf}$ is decomposed as follows

$\hat{H}_{cf}=\sum_{kq}A_{k}^{q}\langle
r^{k}\rangle\Theta_{k}(J)\hat{O}_{k}^{q}$ (4)

where $\hat{O}_{k}^{q}$ is the Stevens operator equivalent, $A_{k}^{q}\langle
r^{k}\rangle$, as explained above, is the standard notation for the crystal-
field parameter for given $k$ and $q$. $\Theta_{k}(J)$ is the Stevens factor
for a given ground state multiplet defined by the quantum number $J$.
$\Theta_{k}(J)$ for $k=$2, 4, and 6 are often designated by $\alpha_{J}$,
$\beta_{J}$, and $\gamma_{J}$, respectively. The Stevens operator equivalents
are more convenient for analytical calculations and somewhat outdated, but
they are still extensively used in the literature. For numerical calculations
it is more convenient to express $\hat{H}_{cf}$ in terms of Wybourne’sWybourne
and Meggers (1965) spherical tensor operators $\hat{C}_{k}^{q}$:

$\hat{H}_{cf}=\sum_{kq}B_{k}^{q}\hat{C}_{k}^{q}$ (5)

where $\hat{C}_{k}^{q}$ are defined by

$C_{k}^{q}(\hat{\boldsymbol{r}})=\sqrt{4\pi/(2k+1)}Y_{kq}(\hat{\boldsymbol{r}})$

Moreover, the CFP can be made real by employing the Hermitian combination of
Wybourne’s operators $\hat{T}_{k}^{q}$ defined by

$\hat{T}_{k}^{0}=\hat{C}_{k}^{0}\textrm{ and }\hat{T}_{k}^{\pm|q|}=\sqrt{\pm
1}\left[\hat{C}_{k}^{-|q|}\pm(-1)^{|q|}\hat{C}_{k}^{|q|}\right]$

$\hat{H}_{cf}$ can then be expressed as

$\hat{H}_{cf}=\sum_{kq}L_{k}^{q}\hat{T}_{k}^{q}$ (6)

with a set of real parameters $L_{k}^{q}$. $L_{k}^{q}$ are linked to the
Stevens CFP $A_{k}^{q}\langle r^{l}\rangle$ by a set of positive prefactors
$\lambda^{q}_{k}=A_{k}^{q}\langle r^{k}\rangle/L_{k}^{q}$. For a more
extensive discussion of CFP conventions see, e.g., Refs. Newman and Ng, 1989;
Mulak and Gajek, 2000; Newman and Ng, 2007.

The number of a priori non-zero CF parameters $A_{k}^{q}\langle r^{k}\rangle$
is constrained by the point-group symmetry of a given rare-earth site. In
particular, in the presence of inversion symmetry,
$V_{ns}(\hat{\boldsymbol{r}})=V_{ns}(-\hat{\boldsymbol{r}})$, only
$A_{k}^{q}\langle r^{k}\rangle$ for even $k$ can be nonzero (cf. Eq. 3). Other
point-group symmetries further reduce the number of relevant $A_{k}^{q}\langle
r^{k}\rangle$. As a consequence, the crystal-field on Sm 4$f$ in SmCo5 can be
fully described with only four CF parameters: $A_{2}^{0}\langle r^{2}\rangle$,
$A_{4}^{0}\langle r^{4}\rangle$, $A_{6}^{0}\langle r^{6}\rangle$, and
$A_{6}^{6}\langle r^{6}\rangle$. In the case of the $R$Fe${}_{12}X$ family,
the relevant parameters are $A_{2}^{0}\langle r^{2}\rangle$, $A_{4}^{0}\langle
r^{4}\rangle$, $A_{4}^{4}\langle r^{4}\rangle$, $A_{6}^{0}\langle
r^{6}\rangle$, and $A_{6}^{4}\langle r^{6}\rangle$ .

In our calculations, we extract the set of parameters $L_{k}^{q}$ (or
$A_{k}^{q}\langle r^{k}\rangle$), as well as $B_{\mathrm{ex}}$ and $\lambda$
by a least-square fit of ab initio $\hat{H}_{\mathrm{1el}}$ (using the usual
Frobenius norm) obtained within DFT+Hub-I (see Eq. 7 in the next section).
Note that one may assign a spin label to the CF parameters in Eqs. 5 and 6,
hence, allowing for different CF potentials for spin up and down electrons. We
found that this improves the fit for spin-polarized $\hat{H}_{\mathrm{1el}}$.

###  II.2 Calculational approach

We employ the DFT+Hub-I approachAichhorn et al. (2016) based on the TRIQS
libraryParcollet et al. (2015) and the full potential linearized augmented
planewave Wien-2kBlaha et al. (2001) band structure code in conjunction with
the projective Wannier-orbitals constructionAmadon et al. (2008); Aichhorn et
al. (2009). The charge-density self-consistencyLechermann et al. (2006);
Pourovskii et al. (2007) is implemented as described in Ref. Aichhorn et al.,
2011. The Hubbard-I impurity solver is provided by the TRIQS library.

The Wannier orbitals representing the rare-earth 4$f$ states are constructed
from the Kohn-Sham bands within the window
$[-\omega_{\mathrm{win}},\omega_{\mathrm{win}}]=\left[-2,2\right]$ eV relative
to the Fermi level. The choice of the half-window size $\omega_{\mathrm{win}}$
is the only significant parameter in our calculations (indeed, the choice of
Hubbard $U$ and Hund’s coupling $J$ has limited impact on the results, as we
demonstrate in Appendix D). In order to construct a complete orthonormal basis
of Wannier orbitals one needs to choose $\omega_{\mathrm{win}}$ large enough
to include at least all 4$f$-like Kohn-Sham bands. Wannier orbitals
constructed with a “small window” leakAichhorn et al. (2016) to neighboring
sites due to hybridization between 4$f$ states and conduction band states. A
larger window results in more localized Wannier orbitals consisting almost
exclusively of the corresponding 4$f$ partial waves inside the rare-earth
atomic sphere Aichhorn et al. (2009, 2016), as discussed in in Sec. IV and
Appendix E below. DFT+Hub-I studies of rare-earth wide-gap insulators show a
rather strong sensitivity of calculated CFP to the window size; less-localized
small window Wannier 4$f$ orbitals result in a better agreement with
experimental CFPPourovskii . In the present case of rare-earth intermetallics
we find a rather weak dependence of CFP to variations of
$\omega_{\mathrm{win}}$ within the reasonable range from 2 to 8 eV, see
Appendix E. Hence, we employ $\omega_{\mathrm{win}}=$2 eV in our calculations
throughout.

In the Hubbard-I approximation the hybridization function is neglected and
solving of the DMFT impurity problem is reduced to the diagonalization of the
atomic Hamiltonian (2). The one-electron part $\hat{H}_{\mathrm{1el}}$ of (2)
is then given byPourovskii et al. (2007)

$\hat{H}_{\mathrm{1el}}=-\mu+\langle H^{ff}\rangle-\Sigma_{\mathrm{DC}}$ (7)

where $\mu$ is the chemical potential, $\langle H^{ff}\rangle$ is the Kohn-
Sham Hamiltonian projected to the basis of 4$f$ Wannier orbitals and summed
over the Brillouin zone, $\Sigma_{\mathrm{DC}}$ is the double counting
correction term for which we employ the fully-localized-limit (FLL) formCzyżyk
and Sawatzky (1994) that is known to work best for localized states such as
4$f$ orbitals. In our calculations, we evaluate the FLL double-counting using
the occupancy of the DMFT local Green’s function, which comes out to be close
to the nominal 4$f$ occupancy of the corresponding 3+ rare-earth ion. If the
nominal occupancy is used in FLL DC instead one obtains almost the same CFP,
with differences no larger than 10 to 20 K. We carry out DFT+Hub-I iterations
until convergence in the total energy with precision $10^{-5}$ Ry is reached
and then extract the CFP from Eq. 7 as described in the previous section.

Self-consistent DFT+Hub-I calculations produce a non-spherical one-electron
Kohn-Sham potential (3), that includes several non-spherical contributions
acting on 4$f$ states: the long-range electrostatic (Madelung) interaction, as
well as the local-density-approximation (LDA) exchange-correlation potential
due to the conduction electrons and 4$f$ states themselves. This last
“intra-4$f$ shell” contribution to the exchange-correlation potential should
be removed within DFT+Hub-I, since the on-site interaction $H_{U}$ between
4$f$ states is already treated explicitly within DMFT. Hence, the “intra-4$f$
shell” contribution in the one-electron part $\hat{H}_{\mathrm{1el}}$ of Eq. 2
due to LDA is counted twice and should be removed by a double-counting
correction. Moreover, this contribution includes the LDA self-interaction
error for localized states directly impacting CFP: for low-lying CF levels,
the self-interaction error will be larger than for less occupied excited CF
states.

In order to reduce the self-interaction error in the CFP we enforce uniform
occupancy of all states within the 4$f$ ground state multiplet in our self-
consistent DFT+Hub-I calculations. To that end, we define the imaginary-
frequency atomic (Hubbard-I) Green’s function at the fermionic Matsubara
frequency $\omega_{n}=(2n+1)\pi T$, where $T$ is the temperature, as follows:

$\displaystyle\begin{split}G^{\mathrm{at}}_{ab}(i\omega_{n})=\frac{1}{M}\sum_{\gamma\in\mathrm{GSM}\atop\delta\notin\mathrm{GSM}}\biggl{(}&\frac{\langle\gamma|f_{a}|\delta\rangle\langle\delta|f_{b}^{\dagger}|\gamma\rangle}{i\omega_{n}-E_{\gamma}+E_{\delta}}\\\
&+\frac{\langle\delta|f_{a}|\gamma\rangle\langle\gamma|f_{b}^{\dagger}|\delta\rangle}{i\omega_{n}+E_{\gamma}-E_{\delta}}\biggr{)}\end{split}$
(8)

where the eigenstates $|\gamma\rangle$ and $|\delta\rangle$ with eigenenergies
$E_{\gamma}$ and $E_{\delta}$ are obtained by diagonalization of Eq. 2 and
belong to the ground-state multiplet (GSM) and excited multiplets
respectively, $a$ and $b$ label 4$f$ orbitals, $M$ is the degeneracy of the
GSM. In other words, to obtain Eq. 8 we substitute the standard Boltzmann
weight $X_{\gamma}=e^{-E_{\gamma}/T}/Z$, where $Z$ is the partition function,
with the uniform weight $\tilde{X}_{\gamma}=1/M$ for the GSM and
$\tilde{X}_{\delta}=0$ for exited multiplets in the spectral representation of
the Green’s function 111We consider the case of $T$ being much lower than the
inter-multiplet splitting, hence, the contribution of excited multiplets into
the partition function $Z$ can be neglected222An equivalent approach would be
to make the one-electron Hamiltonian that serves as an input to the Hubbard I
solver spherically symmetric: both approaches should be equivalent and nearly
equally simple to implement.. In practice, the degeneracy of the ground state
multiplet $M$ is chosen to be the same as for the corresponding free ion,
hence, it is given by Hund’s rules. Therefore, $M=10$ for Nd and $M=6$ for Sm.
The self-energy thus obtained is then plugged back into the self-consistency
cycle. This leads to a spherically-averaged contribution from the 4$f$
orbitals, both inside and outside the rare-earth atomic sphere, while non-
spherical contributions from other states are taken into account. We verified
the validity of this method by calculating the density matrix from the local
Green’s function and transforming it to the relativistic basis of one-electron
$J=\frac{5}{2}$ and $J=\frac{7}{2}$ orbitals. With the averaging, the
resulting density matrix is made of two identity blocks with deviations of the
order of few percent, to be compared with over 50% without the averaging.

Conceptually speaking, our approach amounts to replacing Eq. 7 by

$\displaystyle\hat{H}_{\mathrm{1el}}=-\mu+\langle
H^{ff}\rangle-\Sigma_{\mathrm{DC}}$ $\displaystyle-
v_{\mathrm{KS}}\left[\rho_{spd}+\rho_{4f}\right]$ (9)
$\displaystyle+v_{\mathrm{KS}}\left[\rho_{spd}+\bar{\rho}_{4f}\right]$

where $v_{\mathrm{KS}}\left[\rho\right]$ is the Kohn-Sham potential
($v_{\mathrm{KS}}=v_{\mathrm{Hartree}}+v_{xc}$ ) evaluated from the total
electronic density $\rho$ and then projected to the basis of 4$f$ Wannier
orbitals. $\rho_{4f}$ designates the projected electronic density belonging to
the rare-earth’s 4$f$ orbitals, $\bar{\rho}_{4f}$ is the same density,
spherically averaged, and $\rho_{spd}$ designates all the remaining density,
belonging to all atoms’ $s,p$ and $d$ orbitals.

The same approach is used in the spin-polarized DFT+Hub-I calculations: in
this case the exchange splitting is also removed within the GSM. We found,
however, that this averaging is not sufficient, since the value of the
exchange field within our DFT+Hub-I iterations may become larger than the
inter-multiplet splitting. Hence we also directly remove the 4$f$ spin
polarization from the resulting DFT+Hub-I density matrix. For a given ${\bf
k}$-point the “averaged” density matrix $\tilde{N}^{{\bf k}}$ in the Bloch
basis reads:

$\tilde{N}^{{\bf k}}=N^{{\bf k}}+\frac{1}{2}P^{\dagger}({\bf
k})\left(Tn^{ff}({\bf k})T^{\dagger}-n^{ff}({\bf k})\right)P({\bf k})$ (10)

where $N^{{\bf k}}$ is the density matrix in the Bloch basis calculated as
described in Refs. Aichhorn et al., 2011 and Aichhorn et al., 2016, $P({\bf
k})$ is the projectorAichhorn et al. (2009, 2016) between the Wannier and
Bloch spaces, $n^{ff}({\bf k})$ is the density matrix in the Wannier basis,
$T$ is the time-reversal operator. The averaged density matrix
$\tilde{N}^{{\bf k}}$ is then used to recalculate the electron density at the
next DFT iteration as described in Ref. Aichhorn et al., 2011. The
contribution of 4$f$ states to the spin density and local-spin-density-
approximation (LSDA) exchange field is thus suppressed. The resulting exchange
field is due to the polarization of the transition-metal sublattice, as
expected for hard magnetic rare-earth intermetallics. In contrast, direct
spin-polarized DFT+Hub-I calculations without the averaging would lead to a
large unphysical exchange field on rare-earth sites due to the magnetization
density of 4$f$ electrons themselves.

In appendix A, we benchmark the present method on the prototypical rare-earth
hard magnet SmCo5, for which several measurements of CFP exist, and show good
agreement between calculated and measured CFPs. Moreover, the actual
eigenstates of the Sm 4$f$ shell in SmCo5 obtained within DFT+HubI are also in
very good agreement with previous neutron scattering and mangetic form-factor
measurements, see Appendix B.

###  II.3 Calculational details

The $R$Fe${}_{12}X$ family has the space group $I4/mmm$, with a tetragonal
primitive unit cell. The conventional unit cell, with twice the volume and the
atoms, is orthorhombic. It has equivalent $R$ sites in the corner and the
center at Wyckoff position 2$a$, $X$ interstitial sites between two nearest
$R$ sites on Wyckoff position 2$b$, and contains 24 Fe atoms on three
inequivalent sites, denoted below Fe1, Fe2 and Fe3 at Wyckoff positions 8$j$,
8$i$ and 8$f$ respectively, as displayed in Fig. 1. Calculations are done at
the theoretical lattice constants for $R$Fe${}_{12}X$, summarized in table 1
in the conventional unit cell (from Table II of Ref. Harashima et al., 2015b
and from this work).

Figure 1:  Conventional unit cell of $R$Fe${}_{12}X$. The rare-earth $R$ sites
are yellow, the three types of Fe sites are grey, light and dark blue, and
dopant $X$ sites purple.

The calculated lattice constants agree within 2 % with the measured ones in
the more stable NdFe11Ti(N) and SmFe11Ti(N) compounds.Harashima et al. (2015b)

| Lattice constant (Å)
---|---
Compound |  $a$ |  $c$
NdFe12 |  8.533 |  4.681
NdFe12N |  8.521 |  4.883
NdFe12Li |  8.668 |  4.873
SmFe12 |  8.497 |  4.687
SmFe12N |  8.517 |  4.844
SmFe12Li |  8.640 |  4.863

Table 1:  Conventional unit cell lattice constants used in our calculations.
$b$ = $a$, and the angles are $\alpha=\beta=\gamma=90^{\circ}$

The DFT calculations are performed with spin-orbit coupling included within
the second variational approach. We employ throughout the rotationally-
invariant Coulomb vertex specified by Slater integrals $F^{0}=U=$6.0 eV as
well as $F^{2}=$10.13, $F^{4}=$6.77, and $F^{6}=$5.01 eV corresponding to
Hund’s rule coupling $J_{H}=$0.85 eV. These values of $U$ and $J_{H}$ are in
agreement with those in the literaturePourovskii et al. (2007); Nilsson et al.
(2013); Locht et al. (2016). One may notice, that while the values of $U$ and
$J_{H}$ are important to determine the one-electron spectrum of a material,
they are expected to have a rather small effect on the crystal-field
parameters that we consider in this work. 333In rare-earth ions one typically
has $U\gg J_{H}\gg\lambda$, hence, the atomic multiplet structure is set
within the $LS$-couping scheme. With the crystal-field splitting in
lanthanides being typically much smaller than all those energy scales one may
neglect the mixture between different multiplets, in that case CFP exhibits
only very weak dependence to $U$ and $J_{H}$. We discuss this dependence in
Appendix D. DFT+Hub-I calculations are carried out for the temperature of 290
K.

##  III Results

###  III.1 DFT and DFT + Hubbard I electronic structure of $R$Fe${}_{12}X$

(a) NdFe12N DFT density of states

(b) NdFe12N DFT + Hubbard I spectral function

Figure 2: a. Atom- and orbital-resolved density of states of NdFe12N
calculated with the spin-polarized DFT method. b. Atom- and orbital-resolved
spectral function of the same compound obtained within self-consistent spin-
polarized DFT+Hubbard I. For better readability we take the average over the
three types of Fe atoms: the actual total Fe density of states per unit cell
is three times larger.

We first compare the electronic structure of $R$Fe${}_{12}X$ obtained within
DFT (LSDA) and DFT+Hub-I. A typical DFT density of states (DOS) and a
DFT+Hub-I spectral function for ferromagnetic $R$Fe${}_{12}X$, namely, for
NdFe12N, are shown in Fig. 2. The DFT DOS of Fig. 2a features a strong
polarization of the Fe 3$d$ band. N 2$p$ states are dispersive, with the
bottom of the bands contributing to a peak in the DOS around -6 eV. The Nd
4$f$ band is fully spin-polarized and anti-ferromagnetically aligned to Fe
3$d$, with the total spin moment within the Nd atomic sphere equal to -2.77
$\mu_{B}$, i.e. close to the Hund’s rule value of 3 $\mu_{B}$ for the Nd3+
ion. The Nd majority-spin 4$f$ band is pinned at the Fermi level, its double-
peak structure is due to spin-orbit splitting. This picture of 4$f$ bands
pinned at the Fermi level is qualitatively incorrect and illustrates the
difficulties of DFT with local or semi-local exchange-correlation functionals
to correctly treat strongly-interacting localized valence states.

The spin-polarized DFT+Hub-I spectral function shown in Fig. 2 was calculated
using the averaging approach described in Sec. II.2. It features an almost
fully-polarized Fe 3$d$ band as well as occupied and empty 4$f$ states
separated, to first approximation, by $U$, thus forming lower and upper
Hubbard bands, respectively. The Hubbard bands are split due to the Hund’s
rule and spin-orbit couplings into several manifolds with characteristically
sharp peaks corresponding to transitions from the ground state to different
quasi-atomic multiplets upon electron addition or removal. The 4$f$ multiplet
structure in lanthanides is known to be only weakly sensitive to the
crystalline environment. Indeed, the positions of the Hubbard bands in Fig. 2b
as well as the overall shape of the upper Hubbard band split into two
manifolds of multiplet peaks centered at about 2 and 4 eV are in agreement
with photoemission and inverse-photoemission spectra of the Nd metalLang et
al. (1981). One also sees that the Nd 4$f$ states in DFT+Hub-I are not fully
spin-polarized, in contrast to the DFT case. Indeed the Nd spin moment of
-1.61 $\mu_{B}$ obtained within DFT+Hub-I is only about half of the Hund’s
rule value and is also aligned antiferromagnetically with respect to the spin
moment on iron. The calculated Nd orbital moment is 3.40 $\mu_{B}$. It is
precisely the crystal-field splitting of the Nd 4$f$ shell that prevents the
full saturation of the Nd magnetization.

###  III.2 crystal-field parameters and exchange fields in $R$Fe${}_{12}X$

The calculated CF and exchange fields for Nd and Sm $R$Fe12(N,Li) compounds
are listed in Table 2 and 3, together with the magnetic moments on $R$ and in
the full cell. Comparing the different materials, one sees that $R$Fe12 has
the smallest values of $A_{2}^{0}\langle r^{2}\rangle$ (in absolute value),
while N insertion enhances $A_{2}^{0}\langle r^{2}\rangle$ up to positive
values of about 400 to 600 K. Li insertion has the opposite effect, leading to
large negative $A_{2}^{0}\langle r^{2}\rangle$, in particular for $R=$Nd. We
notice some dependence of the CF parameters $A_{k}^{q}\langle r^{k}\rangle$ on
the spin direction in the ferromagnetic phase. It is mostly weak, of the order
of a few tenths of kelvin for the most important CFP $A_{2}^{0}\langle
r^{2}\rangle$, except in NdFe12N. It can be significant, though, for higher-
order CFP. The magnetic state (paramagnetic of ferromagnetic) has a
significant impact on $A_{2}^{0}\langle r^{2}\rangle$ in some compounds: one
may notice larger values of $A_{2}^{0}\langle r^{2}\rangle$ for paramagnetic
SmFe12(N,Li) than for either spin direction in the ferromagnetic phase.

Finally, the total magnetization appears to be slightly reduced in Sm
compounds, compared to Nd compounds: in the former, the spin magnetic moment
on the rare-earth compensates the orbital magnetic moment, leading to
negligible total moment, while Nd presents a total moment dominated by the
orbital component, and in the same direction as the Fe sublattice
magnetization.

The sign and overall magnitude of our calculated $A_{2}^{0}\langle
r^{2}\rangle$ are in agreement with previous calculations for $R$Fe12(N) in
Ref. Harashima et al., 2015a using the 4$f$-in-core approach, though there are
some differences in the precise values. We obtain a similar value for NdFe12,
a somewhat larger one for NdFe12N, a more negative value for SmFe12 and a
smaller (positive) value for SmFe12N. One may notice that the results in Ref.
Harashima et al., 2015a are quite sensitive to different treatments of the
“tails” of 4$f$ core orbitals: there is no such uncertainty in our approach.

Table 2:  Calculated CF parameters in ferromagnetic (FM) and paramagnetic (PM) NdFe12(N,Li) in kelvin. For the FM case we list the CF parameters for each spin direction. The exchange field in the FM phase (in tesla), the spin and orbital magnetic moments of the rare-earth as well as the total magnetic moment per crystal unit cell (in Bohr magneton $\mu_{B}$) are also listed. | NdFe12 | NdFe12N | NdFe12Li
---|---|---|---
| PM | FM | PM | FM | PM | FM
|  | $\uparrow$ | $\downarrow$ |  | $\uparrow$ | $\downarrow$ |  | $\uparrow$ | $\downarrow$
$A_{2}^{0}\langle r^{2}\rangle$ | -57 | -71 | -116 | 486 | 477 | 653 | -656 | -687 | -742
$A_{4}^{0}\langle r^{4}\rangle$ | -29 | -5 | -1 | 107 | 75 | 112 | -182 | -158 | -186
$A_{4}^{4}\langle r^{4}\rangle$ | -129 | -76 | -270 | 7 | -105 | -141 | -118 | -60 | -228
$A_{6}^{0}\langle r^{6}\rangle$ | 52 | 62 | 54 | 51 | 32 | 63 | -24 | -17 | -31
$A_{6}^{4}\langle r^{6}\rangle$ | 70 | -224 | -107 | -160 | -65 | -91 | 37 | -6 | 96
$B_{\mathrm{ex}}$ ($T$) | - | 265 | - | 217 | - | 410
Nd $M_{\mathrm{spin}}$ | - | -1.48 $\mu_{B}$ | - | -1.61 $\mu_{B}$ | - | -1.69 $\mu_{B}$
Nd $M_{\mathrm{orb}}$ | - | 2.96 $\mu_{B}$ | - | 3.40 $\mu_{B}$ | - | 3.28 $\mu_{B}$
$M_{\mathrm{cell}}$ | - | 26.39 $\mu_{B}$ | - | 29.15 $\mu_{B}$ | - | 27.59 $\mu_{B}$
Table 3:  The same quantities as in Table 2 for SmFe12(N,Li). | SmFe12 | SmFe12N | SmFe12Li
---|---|---|---
| PM | FM | PM | FM | PM | FM
|  | $\uparrow$ | $\downarrow$ |  | $\uparrow$ | $\downarrow$ |  | $\uparrow$ | $\downarrow$
$A_{2}^{0}\langle r^{2}\rangle$ | -32 | -184 | -211 | 249 | 195 | 225 | -458 | -297 | -272
$A_{4}^{0}\langle r^{4}\rangle$ | -11 | -21 | -18 | 99 | 78 | 70 | -116 | -68 | -71
$A_{4}^{4}\langle r^{4}\rangle$ | -215 | -41 | -136 | -122 | 22 | -91 | -124 | 61 | -198
$A_{6}^{0}\langle r^{6}\rangle$ | 47 | 45 | 40 | 71 | 47 | 25 | -13 | -2 | -12
$A_{6}^{4}\langle r^{6}\rangle$ | -85 | -95 | -58 | -184 | -97 | -82 | 44 | 30 | 38
$B_{\mathrm{ex}}$ ($T$) | - | 232 | - | 205 | - | 331
Sm $M_{\mathrm{spin}}$ | - | -3.31 $\mu_{B}$ | - | -2.41 $\mu_{B}$ | - | -3.96 $\mu_{B}$
Sm $M_{\mathrm{orb}}$ | - | 3.29 $\mu_{B}$ | - | 2.35 $\mu_{B}$ | - | 3.60 $\mu_{B}$
$M_{\mathrm{cell}}$ | - | 24.54 $\mu_{B}$ | - | 26.83 $\mu_{B}$ | - | 25.77 $\mu_{B}$

The lowest-order CF parameters $A_{2}^{0}\langle r^{2}\rangle$ and the
corresponding single-ion anisotropy energies $K_{1}$ evaluated using Eq. 1 are
displayed in Fig. 3. One may see that, while NdFe12N and NdFe12Li exhibit
larger $|A_{2}^{0}\langle r^{2}\rangle|$ (upper panel) than their Sm
counterparts, this difference is offset by a larger Stevens prefactor of Sm in
Eq. 1, so that the Sm and Nd-based compounds have a magnetic anisotropy
coefficient $K_{1}$ of similar magnitude. An important difference between Nd
and Sm is the different signs of their Stevens factors $\alpha_{J}$
($\alpha_{J}=-7/1089$ for Nd, $\alpha_{J}=13/315$ for Sm). Consequently, N
insertion leads to a large out-of-plane anisotropy for Nd, but in-plane
anisotropy for Sm. Li has the opposite effect: doping Li into SmFe12 leads to
a rather large out-of-plane anisotropy of SmFe12Li, of comparable magnitude to
that of NdFe12N.

Figure 3: crystal-field parameters $A_{2}^{0}\langle r^{2}\rangle$ (average
for up and down spins in the FM phase) and anisotropy energy $K_{1}$ for _R_
Fe12 _X_ , with _R_ =Nd,Sm and _X_ is either empty, N or Li ($K_{1}$ is
obtained from equation 1).

Performing the averaging over the ground state multiplet as described in Eq. 8
is crucial to obtain reasonable CFP: the lowest-order CFP $A_{2}^{0}\langle
r^{2}\rangle$ is most sensitive to this. The corresponding data without
averaging for NdFe12N are given and discussed in Appendix C.

For the sake of comparison with future experiments we list low-energy
eigenvalues and eigenstates of all $R$Fe${}_{12}X$ compounds in Appendix B. It
is interesting to notice that eigenstates of the ground-state $J=5/2$
multiplet of Sm (Table 6) are often found to exhibit a significant admixture
from exited $J=7/2$ states; the Nd $J=9/2$ states (Table 5) contain a
significantly lower admixture from the first exited multiplet.

A last interesting point is that the exchange fields $B_{\mathrm{ex}}$ on the
rare-earth are enhanced by Li and reduced by N. This is useful because the
exchange field, or exchange coupling between Fe and $R$, is essential for
finite temperature magnetocrystalline anisotropy. The rare-earth-originated
anisotropy becomes ineffective at high temperature, and this threshold
temperature is determined by the exchange coupling $B_{\mathrm{ex}}$. In Fig.
4, we show the difference between the 4$f$ shell atomic energies $E_{\perp}$
and $E_{\parallel}$, computed as
$E=\mathrm{Tr}[\hat{H}e^{-\beta\hat{H}}]/\mathrm{Tr}[e^{-\beta\hat{H}}]$ with
$\hat{H}$ defined in Eq. 2 and the exchange field $B_{\mathrm{ex}}$ is along
the $z$ axis (along the $c$ lattice parameter) and $x$ axis (along the $a$
lattice parameter), respectively. We scale the exchange field
$B_{\mathrm{ex}}$ by a coefficient $M_{Fe}(T)/M_{Fe}(0)$ at non-zero
temperatures, using the measured magnetization ratio of NdFe12N from Hirayama
_et al_ Hirayama et al. (2015). The energy difference plotted in Fig. 4 is
more general than the expression of Eq. 1, because it also contains higher
order CFP and non-zero temperature; to compute $E_{\perp}$ and $E_{\parallel}$
we diagonalize the full Hamiltonian $\hat{H}$, without restricting ourselves
to the ground state multiplet. This gives quite a different picture than Fig.
3: the strongly enhanced exchange coupling due to Li doping causes the
magnetocrystalline anisotropy to persist at much higher temperatures than with
N doping.

Figure 4:  (Color online) Evolution with temperature of the difference in the
4$f$ shell energy $E_{\perp}-E_{\parallel}$ between the moments on $R$ and Fe
aligned perpendicularly and parallel to the $z$ axis, respectively, for
NdFe12N (blue, full line) and SmFe12Li (red, dashed line). Inset:
magnetization fraction of the Fe sublattice in NdFe12N, as a function of
temperature from Hirayama _et al_ Hirayama et al. (2015).

##  IV Discussion: the effect of hybridization with the interstitials

Let us now analyze the mechanisms determining the CFP on the rare-earth site
and, in particular, the impact of the N and Li interstitials on them. We
consider the NdFe${}_{12}X$ ($X=$Ni, Li) compounds as example. The N atom
nominally carries three 2$p$ electrons, but in the $R$Fe12N compounds the N
2$p$ bands are more than half-filled (Fig. 2). To verify this we have also
performed a Bader-charge analysisBader (1990) for NdFe${}_{12}X$ and found 8.3
electrons on N resulting in an ion charge of -1.3. In contrast, the Li atom is
nominally 2$s^{1}$, but it looses its single 2$s$ electron inside the NdFe12
matrix, the corresponding Bader ion charge is +0.7.

In Fig. 5 we display the complex Wannier orbitals constructed for Nd 4$f$
states with window size $\omega_{\textrm{win}}=2$ eV with magnetic quantum
numbers $m=0$ and $m=-3$, in the presence of interstitial N or Li. The
orbitals with $m=\pm 3$ do not point towards the N or Li atom, and leak only
to neighboring Fe atoms. On the other hand, the orbital with $m=0$
(corresponding to $f_{z^{3}}$ cubic orbital) points towards the interstitial
site, and shows strong leakage to the interstitial atom, particularly in the
Li case. The same applies, to a lesser extent, to the orbitals $m=\pm 1$ that
are also pointing towards the interstitials.

The N (Li) insertion has thus two effects on the CFP. The first one is due to
the electrostatic interaction between the 4$f$ electrons and the interstitial
ions. This interaction with the negative N (positive Li) ion pushes the on-
site energies of the $m=-1,0,1$ orbitals, which point towards the
interstitial, to higher (lower) energies.

The second contribution is due to hybridization between the 4$f$ states and
the N 2$p$ (Li 2$s$ and 2$p$) bands, which is expected to mainly affect the
$m=-1,0,1$ orbitals pointing towards the interstitial. Mixing with the empty
Li 2$s$ and 2$p$ bands pushes them to lower energies, while the opposite shift
is induced due to hybridization with mostly filled N 2$p$ located well below
rare-earth 4$f$ states, see Fig. 2. Hence, one sees that both the
electrostatic and hybridization effects act in the same direction, raising the
on-site energies of the $m=-1,0,1$ orbitals in the case of N and lowering them
in the case of Li.

This analysis explains the effect of interstitials on the CFP
$A_{2}^{0}\langle r^{2}\rangle$. Indeed, the contribution due to
$A_{2}^{0}\langle r^{2}\rangle$ into the CF Hamiltonian 6 is $A_{2}^{0}\langle
r^{2}\rangle\hat{T}_{2}^{0}/\lambda_{2}^{0}$, where the matrix of the one-
electron operator $\hat{T}_{2}^{0}/\lambda_{2}^{0}$ reads

$\hat{T}_{2}^{0}/\lambda_{2}^{0}=2\quad\begin{pmatrix}-\frac{1}{3}&&&&\\\
&0&&&&(0)\\\ &&\frac{1}{5}&&\\\ &&&\frac{4}{15}&\\\ &&&&\frac{1}{5}&\\\
&(0)&&&&0&\\\ &&&&&&-\frac{1}{3}\end{pmatrix}\quad$

in the basis of complex 4$f$ orbitals. Hence, the energy level of 4$f$
orbitals $m=\pm 3$ is negatively correlated with $A_{2}^{0}\langle
r^{2}\rangle$, while the energy levels of the orbitals with $m=-1,0,1$ are
positively correlated with $A_{2}^{0}\langle r^{2}\rangle$ (orbitals with
$m=\pm 2$ are unaffected by $\hat{O}_{l}^{m}$). Thus, the effect of N (Li)
insertion is to enhance (reduce) the value of $A_{2}^{0}\langle r^{2}\rangle$.

(a) NdFe12N, $m=-3$ (b) NdFe12N, $m=0$

(c) NdFe12Li, $m=-3$ (d) NdFe12Li, $m=0$

Figure 5: $4f$ Wannier orbitals of NdFe12N and NdFe12Li, for magnetic quantum
number $m=-3$ and $m=0$ and window size $\left[-2,2\right]$ eV. The orbital
with $m=0$ points towards and leak to the N or Li sites, while orbitals with
$m=\pm 3$ do not. All of them leak somewhat to the nearest-neighbor Fe atoms.

One may argue that the Hubbard-I approximation neglects the hybridization
function in solving the quantum impurity problem, hence, hybridization to the
bath is not included explicitly when solving for the self-energy $\Sigma$ in
the DMFT (Hubbard-I) step of our DFT+Hub-I calculations. However, our Wannier
orbitals constructed within the “small” energy window do contain the effect of
hybridization implicitly, which is evidenced by their “leakage” to neighboring
sites due to mixing of rare-earth 4$f$ states with Fe 3$d$, N 2$p$ and Li 2$s$
bands. The real-space Wannier functions of Fig. 5 thus represent a convenient
visualization of hybridization between rare-earth and other orbitals. In order
to quantify the amount of this admixture of the conduction band states we also
expand those extended small-window Wannier orbitals $|w^{\sigma}_{m}({\bf
k})\rangle$ in the basis of localized Wannier functions
$|\tilde{w}^{\alpha\sigma^{\prime}}_{lm^{\prime}}({\bf k})\rangle$ (labeled by
spin $\sigma^{\prime}$, orbital $l$ and magnetic $m^{\prime}$ quantum numbers,
as well as atomic site $\alpha$) constructed within a large energy window for
all relevant bands. In Appendix F we derive the corresponding projection
operators relating $|w^{\sigma}_{m}({\bf k})\rangle$ and
$|\tilde{w}^{\alpha\sigma^{\prime}}_{lm^{\prime}}({\bf k})\rangle$. We employ
it to extract the corresponding contribution $\tilde{\rho}^{m\sigma}_{\alpha
l}(\omega)$ of the shell $l$ on the site $\alpha$ into the spectral function
of the “small-window” 4$f$ orbital $\sigma m$.

The comparison of $\rho^{m\sigma}_{\alpha l}$ for the orbital $m=0$ and $m=3$
are shown in Figs. 7 and 6 for NdFe12Li and NdFe12N, respectively. One may
notice in Fig. 6a that Nd $f_{z^{3}}$ ($m=0$) in NdFe12Li exhibits a strong
hybridization with Li 2$s$ and 2$p$; their contribution is significantly
larger than the admixture of Fe 3$d$ states. We further observe that spin up
states are hybridizing more strongly than spin down states. In contrast, in
the same compound for $m=3$ (Fig. 6b), there is a peak of hybridization with
Fe states but barely any with the Li 2$s$ and 2$p$ ones. The same difference,
but much less pronounced, is noticeable in the case of NdFe12N, see Fig. 7.
Hence, one may conclude, that the effect of the hybridization with the
interstitial on the CF is much larger for Li than for N. In the latter case
the electrostatic shift due to the negative charge on N seems to play the
leading role.

(a)

(b)

Figure 6: a. Projected spectral functions $\rho^{0\sigma}_{\alpha l}(\omega)$
for the 4$f$ orbital with $m=$0 in NdFe12Li, where the atom $\alpha$ and shell
$l$ are given in the legend. The magnitude $\rho^{m\sigma}_{\alpha l}(\omega)$
indicates the amount of admixture of the character $\alpha l$ into a given
4$f$ orbital, for its precise formulation see the text. b. The same for the
4$f$ orbital $m=$3

(a)

(b)

Figure 7: a. Projected spectral functions $\rho^{0\sigma}_{\alpha l}(\omega)$
for the 4$f$ orbital with $m=$0 in NdFe12N. For the notation see caption of
Fig. 6. b. The same for the 4$f$ orbital $m=$3, $\rho^{3\sigma}_{\alpha
l}(\omega)$

##  V Conclusion

In conclusion, we propose a novel first-principles approach for calculating
crystal and exchange fields in rare-earth systems. This approach is formulated
within the DFT+DMFT framework with local correlations on the rare-earth 4$f$
shell treated within the quasi-atomic Hubbard-I approximation. The 4$f$ states
are represented by Wannier functions constructed from a narrow energy range of
Kohn-Sham states of mainly 4$f$ character. We employ a charge-density
averaging that suppresses the contribution due to the self-interaction of the
4$f$ orbitals to the one-electron Kohn-Sham potential. We thus reduce the
effect of this unphysical self-interaction from the crystal-field splitting,
while keeping non-spherical contributions to CFP from other bands. Similarly,
by removing the contribution due to the 4$f$ magnetic density from the
exchange-correlation potential we suppress its unphysical contribution to the
exchange field at the rare-earth site.

The present approach is effectively free from adjustable parameters and can be
applied to evaluate CFP in any localized lanthanide compound. While in the
present work we chose the value for the on-site interaction parameters $U$ and
$J$, they can in principle be evaluated using constrained local-density or
random-phase approximationAryasetiawan et al. (2004). Moreover, we show that
the crystal-field splitting exhibits a rather weak dependence on the value of
$U$ chosen within a reasonable range for lanthanide 4$f$ shells (4 to 8 eV).
Our choice for the local basis representing 4$f$ orbitals, namely, that we
construct it from a narrow range of Kohn-Sham bands with heavy 4$f$ character,
is physically motivated as it allows for the impact of the hybridization on
the CFP being included within DFT+Hub-I.

We apply this approach to evaluate the crystal and exchange-field potentials
as well as the resulting single-ion magnetic anisotropies in several rare-
earth hard-magnetic intermetallics. First, we verify that our ab initio
scheme reproduces the measured crystal-field parameters (CFP) in the well-
known hard magnet SmCo5. We subsequently apply it to prospective rare-earth
hard magnetic intermetallics of the $R$Fe${}_{12}X$ family (where $R=$Nd, Sm
and $X$ can be N, Li or vacancy). Our calculations reproduce the strong out-
of-plane anisotropy of NdFe12N due to a large positive value of the key CFP
$A_{2}^{0}\langle r^{2}\rangle$ induced by insertion of N. Interestingly, we
find that interstitial Li has a strong opposite effect, leading to a large
negative value of $A_{2}^{0}\langle r^{2}\rangle$. We thus predict a strong
out-of-plane anisotropy in the hypothetical compound SmFe12Li. We also find
the anisotropy in SmFe12Li to persist to higher temperatures as compared to
NdFe12N. Hence, Sm-based compounds may represent interesting candidates for
hard-magnetic applications. Of course, the thermodynamic stability of SmFe12Li
and technological feasibility of Li doping still need to be demonstrated by
future studies.

We analyze the effect of N and Li interstitials on $A_{2}^{0}\langle
r^{2}\rangle$ by evaluating the Bader charges as well as by studying the
leakage of 4$f$ Wannier orbitals to interstitial sites and quantifying the
4$f$ hybridization with N 2$p$ and Li 2$s$ states.

Extensions of the present approach beyond the Hubbard-I approximation are
promising for applications to other rare-earth intermetallics. In particular,
a similar DFT+DMFT technique suppressing subtle self-interaction and double-
counting contributions to the Kohn-Sham potential might be necessary to study,
for example, the impact of a spin-polarized transition-metal sublattice on
heavy-fermion behavior in Yb-based intermetallicsMazet et al. (2013, 2015).

###### Acknowledgements.

This work was supported by IDRIS/GENCI Orsay under project t2017091393, the
ECOS-Sud grant A13E04, the French Agence Nationale de la Recherche in the
framework of the international collaborative DFG-ANR project RE-MAP and the
European Research Council under its Consolidator Grant scheme (project
617196). LP acknowledges the financial support of the Ministry of Education
and Science of the Russian Federation (Grant No. 14.Y26.31.0005). TM
acknowledges the financial support by the Ministry of Education, Culture,
Sports, Science and Technology (MEXT) of Japan as a social and scientific
priority issue (Creation of new functional devices and high-performance
materials to support next-generation industries, CDMSI) to be tackled by using
post-K computer, and also by MEXT KAKENHI Grant Number 16H06345. We thank
Tilmann Hickel, Halil Soezen, Antoine Georges and Dominique Givord for useful
discussions.

##  Appendix A crystal-field parameters in SmCo5

SmCo5 has been studied more extensively than other hard magnetic rare-earth
intermetallics, so ample experimental data is available in this case. In
particular, several groups estimated the CF parameters using inelastic neutron
scattering or magnetization measurements. Therefore, this compound is a good
benchmark to test our approach. SmCo5 has already been studied within DFT+Hub
I to evaluate its ground state magnetization and photoemission spectraGrånäs
et al. (2012), but the CF parameters were not calculated in this work.

Figure 8: Atom-resolved spectral function of SmCo5 obtained within self-
consistent spin-polarized DFT+Hubbard I (full lines). The two inequivalent Co
types are summed to give the total Co 3$d$ spectral function. The occupied
part of the experimental spectrum of SmCo5 (light blue from Ref. Cuthill et
al., 1974 and purple from Ref. Laslo et al., 2014) and the full experimental
spectrum of metal Sm (red, from Ref. Lang et al., 1981) are shown for
comparison in dotted lines.

The calculated spectral function of SmCo5 is shown in Fig. 8. We find a total
magnetic moment on the 4$f$ shell of Sm of 0.42 $\mu_{B}$, antiparallel with
the Co moments. This compares well with the measured value of 0.38 $\mu_{B}$
at 4.2 KGivord et al. (1979).

The calculated CFP and exchange fields for SmCo5 are listed in Table 4,
together with experimental data. The calculations on SmCo5 are done at the
experimental lattice constants.

| PM | FM | Tils et al. | Zhao et al. | Givord et al. | Sankar et al. | Bushow et al. | Richter et al. | Hummler et al. | Novak et al.
---|---|---|---|---|---|---|---|---|---|---
|  | $\uparrow$ | $\downarrow$ |  †Ref. Tils et al.,1999 |  †Ref. Tie-Song et al.,1991 |  †Ref. Givord et al.,1979 |  †Ref. Sankar et al.,1975 |  †Ref. Buschow et al.,1974 | Ref. Richter et al.,1995 | Ref. Hummler and Fähnle,1996 | Ref. Novak and Kuriplach,1994
$A_{2}^{0}\langle r^{2}\rangle$ | -140 | -313 | -262 | -326 | -330 | -200 | -420 | -180 | -755 | -509 | -160
$A_{4}^{0}\langle r^{4}\rangle$ | -40 | -40 | -55 | - | -45 | 0 | 25 | 0 | -37 | -20 | -33
$A_{6}^{0}\langle r^{6}\rangle$ | 33 | 35 | 25 | - | 0 | 50 | 0 | 0 | 11 | 2 | 40
$A_{6}^{6}\langle r^{6}\rangle$ | -684 | -731 | -593 | - | 0 | 0 | 6 | 0 | 290 | -55 | 168
$B_{\mathrm{ex}}$ ($T$) | - | 227 | 260 | 327.5 | 260.5 | 357 | 298 | - | 279 | -

Table 4:  Calculated CF parameters in ferromagnetic (FM) and paramagnetic (PM)
SmCo5 in kelvin. For the FM case we list the CF parameters for each spin
direction. The exchange field in the FM phase (in tesla) is also listed. For
comparison, measured and calculated values from several groups are also given.
References corresponding to an experimental work are denoted by the symbol †.

One may notice that the CF parameter $A_{2}^{0}\langle r^{2}\rangle$ exhibits
a strong dependence on the spin polarization; it is about twice larger in the
FM phase. For other CF parameters this dependence is small.

Our results for $A_{2}^{0}\langle r^{2}\rangle$ are in good agreement with the
experimental (rather wide) range from about -180 to -420 K. The calculated
$B_{\mathrm{ex}}$ also agrees rather well with the experimental range from 260
to 360 T. One may notice that the experimental measurements were performed at
room temperature, hence, in ferromagnetic SmCo5. Also, the most recent
experimental valuesTils et al. (1999) of $A_{2}^{0}\langle r^{2}\rangle$ are
in very good agreement with our results for the FM phase.

The main discrepancy between our theoretical and experimental CFP lies in the
large value that we find for $A_{6}^{6}\langle r^{6}\rangle$. The high-order
CF parameters are usually assumed to be rather small in SmCo5. However, as
noted in Ref. Hummler and Fähnle, 1996, experimental inelastic neutron and
susceptibility data are not particularly sensitive to those high-order
parameters. Hence, they are often assumed to be small from the onset and
neglected in the fitting procedure.

In order to facilitate the reproducibility of our calculations, we provide
below the one-electron Hamiltonian of Eq. 7 for a converged, ferromagnetic
calculation of SmCo5, used to obtain the CFP of Table 4.

$\hat{H}_{\mathrm{1el}}^{\uparrow\uparrow}=\quad\begin{pmatrix}-26.5763&0&0&0&0&0&0.0235\\\
0&-26.5003&0&0&0&0&0\\\ 0&0&-26.4465&0&0&0&0\\\ 0&0&0&-26.3511&0&0&0\\\
0&0&0&0&-26.2796&0&0\\\ 0&0&0&0&0&-26.1675&0\\\
0.0235&0&0&0&0&0&-26.0800\end{pmatrix}\quad$
$\hat{H}_{\mathrm{1el}}^{\downarrow\downarrow}=\quad\begin{pmatrix}-26.1093&0&0&0&0&0&0.0191\\\
0&-26.1926&0&0&0&0&0\\\ 0&0&-26.3023&0&0&0&0\\\ 0&0&0&-26.3802&0&0&0\\\
0&0&0&0&-26.4688&0&0\\\ 0&0&0&0&0&-26.5252&0\\\
0.0191&0&0&0&0&0&-26.6069\end{pmatrix}\quad$
$\hat{H}_{\mathrm{1el}}^{\uparrow\downarrow}=(\hat{H}_{\mathrm{1el}}^{\downarrow\uparrow})^{\dagger}=\quad\begin{pmatrix}0&0.2032&0&0&0&0&0\\\
0&0&0.2632&0&0&0&0\\\ 0&0&0&0.2893&0&0&0\\\ 0&0&0&0&0.2886&0&0\\\
0&0&0&0&0&0.2633&0\\\ -0.0002&0&0&0&0&0&0.2036\\\
0&-0.0003&0&0&0&0&0\end{pmatrix}\quad$

where

$\hat{H}_{\mathrm{1el}}=\begin{pmatrix}\hat{H}_{\mathrm{1el}}^{\uparrow\uparrow}&\hat{H}_{\mathrm{1el}}^{\uparrow\downarrow}\\\
\hat{H}_{\mathrm{1el}}^{\downarrow\uparrow}&\hat{H}_{\mathrm{1el}}^{\downarrow\downarrow}\end{pmatrix}$

##  Appendix B Eigenenergies and eigenstates of the 4$f$ shells in SmCo5 and
$R$Fe${}_{12}X$

In this Appendix we present the actual converged eigenfunctions and
eigenstates of the 4$f$ shell obtained within our DFT+HubI approach. In tables
5 and 6, we list those eigenenergies and the corresponding wavefunctions for
the ground-state multiplet, as well as for the lowest-energy state of the
first exited multiplet, in all the materials considered. The eigenenergies are
given with respect to the ground state. The eigenvalues are expanded in the
basis of total angular momentum $J$ as
$\sum_{J,m_{J}}a(J,m_{J})|J;m_{J}\rangle$; we include all contributions with
$|a(J,m_{J})|>0.03$.

One sees that the eigenstates of Sm belonging to the ground-state multiplet
feature a rather significant admixture of the exited $J=7/2$ multiplet. The
inter-multiplet mixing is markedly lower in the case of Nd.

To our awareness, only 4$f$ eigenstates in SmCo5 have been measured to date.
Our calculated intra- and inter-multiplet splittings are in good agreement
with the results of of Tils _et al._ (Ref. Tils et al., 1999) and Givord _et
al._ (Refs. Givord et al., 1979), see the lowest panel of Table 6. Moreover,
the actual eigenstates and their order are also in very good agreement with
the magnetic form-factor measurementsGivord et al. (1979); Laforest (1981),
especially for the lowest-energy states. 444The opposite sign of the $J=7/2$
contributions in Refs. Givord et al., 1979 and Laforest, 1981 with respect to
ours is due to their choice for the coordination frame, with the $z$ axis
aligned to the spin moment of Sm. We align $z$ with the spin of Co sublattice,
hence opposite to the Sm spin. Flipping the direction of the $z$-axis leads to
the change in the relative sign between the $J=5/2$ and $J=7/2$ contributions.
There is a mismatch between the signs of $m_{J}$ in $|J,m_{J}\rangle$ reported
in Refs. Givord et al., 1979 and Laforest, 1981. The correct signs are
reported in Ref. Laforest, 1981 (D. Givord, private communication).

NdFe12
---
GSM | 0 |  $0.996|9/2;+9/2\rangle+0.081|11/2;+9/2\rangle$
| 8 |  $0.995|9/2;+5/2\rangle+0.094|11/2;+5/2\rangle$
| 18 |  $0.989|9/2;+3/2\rangle+0.145|11/2;+3/2\rangle-0.037|9/2;-5/2\rangle$
| 18 |  $0.999|9/2;+7/2\rangle+0.042|11/2;+7/2\rangle$
| 39 |  $0.989|9/2;+1/2\rangle+0.145|11/2;+1/2\rangle-0.032|9/2;-7/2\rangle$
| 51 |  $0.996|9/2;-3/2\rangle+0.085|11/2;-3/2\rangle$
| 51 |  $0.994|9/2;-1/2\rangle+0.105|11/2;-1/2\rangle$
| 56 |  $0.992|9/2;-5/2\rangle+0.121|11/2;-5/2\rangle+0.038|9/2;+3/2\rangle$
| 82 |  $0.985|9/2;-7/2\rangle+0.166|11/2;-7/2\rangle+0.033|9/2;+1/2\rangle$
| 88 |  $0.998|9/2;-9/2\rangle+0.057|11/2;-9/2\rangle$
FEM | 309 |  $0.987|11/2;+5/2\rangle+0.126|13/2;+5/2\rangle-0.093|9/2;+5/2\rangle$
|  |  $-0.046|11/2;-3/2\rangle$
NdFe12N
GSM | 0 |  $0.994|9/2;+9/2\rangle+0.092|11/2;+9/2\rangle-0.054|9/2;+1/2\rangle$
| 32 |  $0.998|9/2;+5/2\rangle+0.065|11/2;+5/2\rangle$
| 38 |  $0.994|9/2;+3/2\rangle+0.106|11/2;+3/2\rangle$
| 41 |  $0.997|9/2;+7/2\rangle-0.065|9/2;-1/2\rangle$
| 56 |  $0.993|9/2;+1/2\rangle+0.101|11/2;+1/2\rangle+0.055|9/2;+9/2\rangle$
| 66 |  $0.998|9/2;-3/2\rangle+0.052|11/2;-3/2\rangle$
| 66 |  $0.966|9/2;-1/2\rangle-0.242|9/2;-9/2\rangle+0.064|9/2;+7/2\rangle+0.061|11/2;-1/2\rangle$
| 67 |  $0.996|9/2;-5/2\rangle+0.077|11/2;-5/2\rangle$
| 73 |  $0.970|9/2;-9/2\rangle+0.241|9/2;-1/2\rangle$
| 80 |  $0.996|9/2;-7/2\rangle+0.093|11/2;-7/2\rangle$
FEM | 301 |  $0.994|11/2;+11/2\rangle+0.091|13/2;+11/2\rangle-0.051|11/2;+3/2\rangle$
NdFe12Li
GSM | 0 |  $0.993|9/2;+7/2\rangle+0.116|11/2;+7/2\rangle$
| 16 |  $0.990|9/2;+5/2\rangle+0.137|11/2;+5/2\rangle$
| 24 |  $0.994|9/2;+9/2\rangle+0.105|9/2;+1/2\rangle$
| 31 |  $0.989|9/2;+3/2\rangle+0.137|11/2;+3/2\rangle-0.041|9/2;-5/2\rangle$
| 40 |  $0.983|9/2;+1/2\rangle+0.146|11/2;+1/2\rangle-0.107|9/2;+9/2\rangle$
| 51 |  $0.986|9/2;-1/2\rangle+0.163|11/2;-1/2\rangle$
| 65 |  $0.985|9/2;-3/2\rangle+0.166|11/2;-3/2\rangle$
| 80 |  $0.988|9/2;-5/2\rangle+0.148|11/2;-5/2\rangle+0.042|9/2;+3/2\rangle$
| 93 |  $0.990|9/2;-7/2\rangle+0.135|11/2;-7/2\rangle$
| 126 |  $0.984|9/2;-9/2\rangle+0.179|11/2;-9/2\rangle$
FEM | 319 |  $0.983|11/2;+7/2\rangle+0.134|13/2;+7/2\rangle-0.116|9/2;+7/2\rangle$
|  |  $+0.042|11/2;-1/2\rangle$
Table 5:  Energies (in meV) and wavefunctions (expanded in the total angular
momentum $J$ basis) of the atomic eigenstates in the ground state multiplet
(GSM) and of the lowest eigenstate of the first excited multiplet (FEM) for
NdFe${}_{12}X$.  SmFe12
---
GSM | 0 |  $0.986|5/2;+5/2\rangle+0.164|7/2;+5/2\rangle$
| 27 |  $0.983|5/2;+3/2\rangle+0.183|7/2;+3/2\rangle$
| 47 |  $0.968|5/2;+1/2\rangle+0.245|7/2;+1/2\rangle+0.043|9/2;+1/2\rangle$
| 69 |  $0.977|5/2;-1/2\rangle+0.209|7/2;-1/2\rangle+0.041|9/2;-1/2\rangle$
| 85 |  $0.976|5/2;-3/2\rangle+0.218|7/2;-3/2\rangle$
| 99 |  $0.987|5/2;-5/2\rangle+0.159|7/2;-5/2\rangle+0.032|9/2;-5/2\rangle$
FEM | 191 |  $0.987|7/2;+7/2\rangle+0.154|9/2;+7/2\rangle+0.034|7/2;-1/2\rangle$
SmFe12N
GSM | 0 |  $0.997|5/2;+5/2\rangle+0.079|7/2;+5/2\rangle$
| 2 |  $0.990|5/2;+3/2\rangle+0.143|7/2;+3/2\rangle$
| 13 |  $0.976|5/2;+1/2\rangle+0.216|7/2;+1/2\rangle$
| 31 |  $0.982|5/2;-1/2\rangle+0.186|7/2;-1/2\rangle+0.035|9/2;-1/2\rangle$
| 48 |  $0.978|5/2;-3/2\rangle+0.209|7/2;-3/2\rangle$
| 75 |  $0.971|5/2;-5/2\rangle+0.238|7/2;-5/2\rangle$
FEM | 176 |  $0.988|7/2;+5/2\rangle+0.129|9/2;+5/2\rangle-0.078|5/2;+5/2\rangle$
SmFe12Li
GSM | 0 |  $0.977|5/2;+5/2\rangle+0.212|7/2;+5/2\rangle+0.037|9/2;+5/2\rangle$
| 39 |  $0.965|5/2;+3/2\rangle+0.256|7/2;+3/2\rangle+0.047|9/2;+3/2\rangle$
| 67 |  $0.960|5/2;+1/2\rangle+0.278|7/2;+1/2\rangle+0.042|9/2;+1/2\rangle$
| 90 |  $0.946|5/2;-1/2\rangle+0.318|7/2;-1/2\rangle+0.051|9/2;-1/2\rangle$
| 116 |  $0.941|5/2;-3/2\rangle+0.330|7/2;-3/2\rangle+0.068|9/2;-3/2\rangle$
| 137 |  $0.974|5/2;-5/2\rangle+0.221|7/2;-5/2\rangle+0.053|9/2;-5/2\rangle$
FEM | 202 |  $0.977|7/2;+7/2\rangle+0.211|9/2;+7/2\rangle+0.035|11/2;+7/2\rangle$
SmCo5
GSM | 0 (0 / 0) |  $0.984|5/2;+5/2\rangle+0.171|7/2;+5/2\rangle$
| 33 (31 / 28) |  $0.983|5/2;+3/2\rangle+0.181|7/2;+3/2\rangle$
| 52 ( - / 47) |  $0.973|5/2;+1/2\rangle+0.225|7/2;+1/2\rangle+0.033|9/2;+1/2\rangle$
| 71 ( - / 73) |  $0.977|5/2;-1/2\rangle+0.209|7/2;-1/2\rangle+0.035|9/2;-1/2\rangle$
| 86 ( - / 91) |  $0.977|5/2;-3/2\rangle+0.211|7/2;-3/2\rangle-0.032|9/2;+9/2\rangle$
| 95 ( - / 109) |  $0.989|5/2;-5/2\rangle+0.122|7/2;-5/2\rangle+0.073|7/2;+7/2\rangle$
|  |  $+0.033|9/2;-5/2\rangle$
FEM | 188 (166 / - ) |  $0.963|7/2;+7/2\rangle+0.187|7/2;-5/2\rangle+0.164|9/2;+7/2\rangle$
|  |  $-0.093|5/2;-5/2\rangle$
Table 6:  Same as Table 5, for SmFe12X and SmCo5. For SmCo5, between brackets
next to the calculated energies: energies of the atomic eigenstates measured
by Tils _et al._ (left, Ref. Tils et al., 1999) and Givord _et al._ (right,
Ref. Givord et al., 1979). Note that only Tils _et al._ directly measure the
eigenenergies, while Givord _et al._ obtain them from an atomic Hamiltonian
fitted to reproduce the measured magnetic form factor.

##  Appendix C Importance of the charge averaging

In this Appendix we explicitly demonstrate the effect of averaging of 4$f$
charge density (eq. 8) by comparing the CFP calculated with and without this
averaging (but in both cases the 4$f$ magnetic density is suppressed following
Eq. 10) in two materials, NdFe12N and SmCo5, that are known to have an out-of-
plane magnetic anisotropy.

The corresponding values are displayed in Table 7. One sees that the
difference is largest for the lowest-order CFP $A_{2}^{0}\langle
r^{2}\rangle$, where calculations without averaging lead to the wrong sign
with respect to experiment (suggesting in-plane anisotropy in both cases).
Hence, the proper averaging of 4$f$ charge density is crucial for a correct
description of the single-ion anisotropy. For the higher order terms the
difference between two approaches is smaller. This suggests that the self-
interaction contribution in the CFP has predominantly $l=2$ symmetry.

Table 7:  crystal-field parameters and exchange field in NdFe12N and SmCo5 in the ferromagnetic phase, calculated with and without averaging over the ground state multiplet.  | NdFe12N | SmCo5
---|---|---
| with | without | with | without
|  $\uparrow$ |  $\downarrow$ |  $\uparrow$ |  $\downarrow$ |  $\uparrow$ |  $\downarrow$ |  $\uparrow$ |  $\downarrow$
$A_{2}^{0}\langle r^{2}\rangle$ |  477 |  653 |  -190 |  26 |  -313 |  -262 |  278 |  331
$A_{4}^{0}\langle r^{4}\rangle$ |  75 |  112 |  30 |  82 |  -40 |  -55 |  -30 |  -37
$A_{4}^{4}\langle r^{4}\rangle$ |  -105 |  -141 |  -65 |  -124 |  0 |  -0 |  0 |  0
$A_{6}^{0}\langle r^{6}\rangle$ |  32 |  63 |  27 |  64 |  35 |  25 |  38 |  25
$A_{6}^{4}\langle r^{6}\rangle$ |  -65 |  -91 |  -61 |  -112 |  0 |  0 |  0 |  0
$A_{6}^{6}\langle r^{6}\rangle$ |  0 |  -0 |  0 |  0 |  -731 |  -593 |  -945 |  -806
$B_{\mathrm{ex}}$ ($T$) | 217 | 206 | 227 | 235

##  Appendix D Dependence of results on Coulomb $U$ and Hund’s $J_{H}$

To perform DFT+DMFT calculations, we have to choose a value for the on-site
screened Coulomb interaction parameter $U$ and for the Hund’s coupling
parameter $J_{H}$. Several methods have been developed in order to compute
those parameters from first principles, most notably the constrained local
density approximationDederichs et al. (1984) and, more recently, the
constrained random phase approximationAryasetiawan et al. (2004).

In the present work, however, we do not attempt a first principles
determination. We use $U=6$ eV and $J_{H}=0.85$ eV because these values have
given satisfactory results in other calculations on rare-earth
materialsPourovskii et al. (2007); Locht et al. (2016). They are also in line
with reported values calculated from first principlesNilsson et al. (2013).
Nevertheless, it is preferable that results obtained by our calculation scheme
do not depend too strongly on the value of $U$ and $J_{H}$. In Fig. 9, we show
that the dependence of the CFP $A_{2}^{0}\langle r^{2}\rangle$ in NdFe12N is
very moderate, as long as the values of U and J are chosen within a reasonable
ranges for rare-earth ions.

Furthermore, we observe that smaller values of $U$ lead to slightly larger
values of $A_{2}^{0}\langle r^{2}\rangle$: this is not surprising if we keep
in mind that a large $U$ is favorable to a strong localization of the 4$f$
electrons, hence to a weaker coupling to the crystal-field.

Figure 9: CFP $A_{2}^{0}\langle r^{2}\rangle$ in NdFe12N as a function of $U$
for $J_{H}=0.85$ eV (left-hand panel), and as a function of $J_{H}$ for $U=6$
eV (right-hand panel). Our reference values are ($U=6$ eV, $J_{H}=0.85$ eV).

##  Appendix E Dependence of results on window size

Another important parameter of our calculations is the size of the window
around the Fermi level that we use to construct the $4f$ Wannier functions. In
Fig. 10 we compare the Wannier orbitals constructed for the same orbital $m=0$
in NdFe12Li for two different window sizes: a small window with
$\omega_{\textrm{win}}=2$ eV, and a large one with $\omega_{\textrm{win}}=20$
eV. For the large window, the Wannier orbital (WO) takes essentially pure Nd
4$f$ orbital character, while the small-window WO leaks significantly to
neighboring sites, in particular, to Li.

(a) $\left[-20,20\right]$ eV window

(b) $\left[-2,2\right]$ eV window

Figure 10: NdFe12Li $4f$ Wannier orbital $m=0$ constructed with a large window
$\left[-20,20\right]$ eV (left) and a small window $\left[-2,2\right]$ eV
(right). The use of a large window essentially removes all hybridization
between the rare-earth and neighboring atoms.

The effect of the window size on the CF parameters is shown more
quantitatively in Fig. 11, which displays those parameters computed for
several window choices $[-\omega_{\mathrm{win}},\omega_{\mathrm{win}}]$ for
different materials. The smallest window size of $\omega_{\mathrm{win}}=$2 eV
is required to enclose all the 4$f$-like bands, increasing it to 4 eV includes
most of the Fe states and part of the N or Li states inside the window. The
largest size of 20 eV gives Wannier functions with essentially pure orbital
character. One may notice a relatively mild dependence of the CFP on the
choice of the window up to $\omega_{\mathrm{win}}=$8 eV.

Figure 11: Absolute value of the CFP $A_{2}^{0}\langle r^{2}\rangle$ as a
function of window size $\omega_{\textrm{win}}$ for NdFe12, NdFe12N, NdFe12Li
and SmCo5.

##  Appendix F Projection of extended Wannier orbitals to localized Wannier
basis

In this Appendix we derive the projection operator between localized and
extended Wannier spaces. A set of Wannier-like functions
$|\tilde{w}^{\alpha\sigma}_{lm}({\bf k})\rangle$ is constructed for an atom
$\alpha$ of the unit cell and quantum numbers $(lm\sigma)$ as a combinations
of Kohn-Sham Bloch waves for a range of bands within the chosen energy window
$\tilde{\cal{W}}$:

$|\tilde{w}^{\alpha\sigma}_{lm}({\bf
k})\rangle=\sum_{\nu\in\cal{\tilde{W}}}\tilde{P}^{\alpha\sigma}_{lm\nu}({\bf
k})|\phi^{k}_{\nu}\rangle$ (11)

where $\phi^{k}_{\nu}$ are the Bloch functions and
$\tilde{P}^{\alpha\sigma}_{lm\nu}({\bf k})$ is the corresponding matrix
element of the projector constructed as described in Refs. Amadon et al., 2008
and Aichhorn et al., 2009. The corresponding real-space Wannier functions are
then obtained by a Fourier transformation

$\tilde{w}^{\alpha\sigma}_{lm}({\bf r})=\sum_{{\bf k}}e^{-i{\bf k}.{\bf
r}}|\tilde{w}^{\alpha\sigma}_{lm}({\bf k})\rangle.$ (12)

We assume that the window $\tilde{\cal{W}}$ in Eq. 11 is large, i.e. that it
includes both rare-earth 4$f$ states and all relevant valence bands that are
expected to hybridize with them. In result, with such a large-window
construction one obtains a set of mutually-orthogonal and rather well
localized Wannier orbtials (WO). In particular, the large-window 4$f$ WOs
$\tilde{w}^{\alpha\sigma}_{lm}({\bf r})$ almost do not leak onto neighboring
sites, as discussed in the previous section, see Fig. 10a. If one constructs
as many WOs as the number of Kohn-Sham bands within $\tilde{\cal{W}}$ then the
projection operator $\tilde{P}({\bf k})$ is just a unitary transformation,
hence, Eq. 11 can be inverted

$|\phi^{k}_{\nu}\rangle\approx\sum_{\alpha\sigma
lm}\left[\tilde{P}^{\alpha\sigma}_{lm\nu}({\bf
k})\right]^{*}|\tilde{w}^{\alpha\sigma}_{lm}({\bf k})\rangle,$ (13)

where the equality is approximate because high-energy empty bands usually
cross and, hence, one cannot generally chose such a window as to have the same
number of bands for all ${\bf k}$-points. However, those high-energy states
are far from the relevant region close to the Fermi level, and if one applies
Eq. 13 to the bands within a small window $\cal{W}$ around the the Fermi
energy the resulting small non-unitarity of $P({\bf k})$ can be neglected.

Alternatively, one may construct 4$f$ Wannier orbitals from the bands within
that small window $\cal{W}$ enclosing mainly 4$f$-like Kohn-Sham bands:

$|w^{\sigma}_{m}({\bf k})\rangle=\sum_{\nu\in\cal{W}}P^{\sigma}_{m\nu}({\bf
k})|\phi^{k}_{\nu}\rangle,$ (14)

where 4$f$ WOs are constructed for the single rare-earth site in the unit cell
for the compounds under consideration. Hence, the site and $l$ labels are
suppressed in $|w^{\sigma}_{m}({\bf k})\rangle$. The resulting small-window
WOs are rather extended in real space, as one sees in Figs. 5 and 10b.

Inserting the expansion Eq. 13 of the KS states $|\phi^{k}_{\nu}\rangle$ into
Eq. 14 one obtains

$\displaystyle|w^{\sigma}_{m}({\bf k})\rangle$
$\displaystyle=\sum_{\nu\in\textrm{\cal{W}}}\sum_{lm^{\prime}\sigma^{\prime}}P^{\sigma}_{m\nu}({\bf
k})\left[\tilde{P}^{\alpha\sigma^{\prime}}_{lm^{\prime}\nu}({\bf
k})\right]^{*}|\tilde{w}^{\alpha\sigma^{\prime}}_{lm^{\prime}}({\bf
k})\rangle$
$\displaystyle=\sum_{\alpha\sigma^{\prime}lm^{\prime}}U^{\sigma,\alpha\sigma^{\prime}}_{m,lm^{\prime}}({\bf
k})|\tilde{w}^{\alpha\sigma^{\prime}}_{lm^{\prime}}({\bf k})\rangle$ (15)

where

$U^{\sigma,\alpha\sigma^{\prime}}_{m,l^{\prime}m^{\prime}}({\bf
k})=\sum_{\nu\in\cal{W}}P^{\sigma}_{m\nu}({\bf
k})\left[\tilde{P}^{\alpha\sigma^{\prime}}_{lm^{\prime}\nu}({\bf
k})\right]^{*}.$ (16)

We use these projectors
$U^{\sigma,\alpha\sigma^{\prime}}_{m,l^{\prime}m^{\prime}}({\bf k})$ to
project the 4$f$ spectral function computed in the small-window WO basis on
large-window localized WOs representing other states (Fe 3$d$, N 2$p$, Li 2$s$
and so on). Namely, having obtained the real-axis lattice Green’s function in
the small-window Wannier basis for the orbital ($\sigma m$) of the 4$f$ shell,
$G_{m\sigma}({\bf k},\omega+i\delta)$, as well as the corresponding partial
spectral function
$\rho_{m\sigma}(\omega)=-\frac{1}{\pi}\textrm{Im}G_{m\sigma}({\bf
k},\omega+i\delta)$, we compute the different orbital contributions into it as
follows:

$\displaystyle\tilde{\rho}^{m\sigma}_{\alpha
l}(\omega)=-\frac{1}{\pi}\textrm{Im}\sum_{k}\sum_{m^{\prime}\sigma^{\prime}}$
$\displaystyle\left[U^{\sigma,\alpha\sigma^{\prime}}_{m,lm^{\prime}}({\bf
k})\right]^{*}$ (17) $\displaystyle\times G_{m\sigma}({\bf
k},\omega+i\delta)U^{\sigma,\alpha\sigma^{\prime}}_{m,lm^{\prime}}({\bf k})$

where $\tilde{\rho}^{m\sigma}_{\alpha l}(\omega)$ is the fraction of the 4$f$
spectral function of orbital index ($\sigma m$) with the character $(\alpha
l)$. Using Eq. 13 and the orthonormality of small-window WOs

$\langle w^{\sigma}_{m}({\bf k})|w^{\sigma^{\prime}}_{m^{\prime}}({\bf
k})\rangle=\delta_{mm^{\prime}}\delta_{\sigma\sigma^{\prime}}=\sum_{\nu}\left[P^{\sigma}_{m\nu}({\bf
k})\right]^{*}P^{\sigma^{\prime}}_{m^{\prime}\nu}({\bf k})$

one may easily show the completeness of the expansion (17)

$\sum_{\alpha l}\tilde{\rho}^{m\sigma}_{\alpha
l}(\omega)=\rho_{m\sigma}(\omega).$

## References

  * Strnat (1970) K. Strnat, IEEE Transactions on Magnetics 6, 182 (1970).
  * Herbst (1991) J. Herbst, Reviews of Modern Physics 63, 819 (1991).
  * Körner et al. (2016) W. Körner, G. Krugel, and C. Elsässer, Scientific reports 6, 24686 (2016).
  * Suzuki et al. (2016) S. Suzuki, T. Kuno, K. Urushibata, K. Kobayashi, N. Sakuma, K. Washio, M. Yano, A. Kato, and A. Manabe, Journal of Magnetism and Magnetic Materials 401, 259 (2016).
  * Harashima et al. (2015a) Y. Harashima, K. Terakura, H. Kino, S. Ishibashi, and T. Miyake, Physical Review B 92, 184426 (2015a).
  * Hirayama et al. (2015) Y. Hirayama, Y. Takahashi, S. Hirosawa, and K. Hono, Scripta Materialia 95, 70 (2015).
  * Miyake et al. (2014) T. Miyake, K. Terakura, Y. Harashima, H. Kino, and S. Ishibashi, Journal of the Physical Society of Japan 83, 043702 (2014).
  * Coey (2011) J. Coey, IEEE Transactions on Magnetics 47, 4671 (2011).
  * Coey (1996) J. M. D. Coey, _Rare-earth iron permanent magnets_ , 54 (Oxford University Press, 1996).
  * Gutfleisch et al. (2011) O. Gutfleisch, M. A. Willard, E. Brück, C. H. Chen, S. Sankar, and J. P. Liu, Advanced materials 23, 821 (2011).
  * Matsumoto et al. (2016) M. Matsumoto, H. Akai, Y. Harashima, S. Doi, and T. Miyake, Journal of Applied Physics 119, 213901 (2016).
  * Harashima et al. (2016) Y. Harashima, K. Terakura, H. Kino, S. Ishibashi, and T. Miyake, Journal of Applied Physics 120, 203904 (2016).
  * Hutchings (1964) M. T. Hutchings, Solid state physics 16, 227 (1964).
  * Newman and Ng (1989) D. Newman and B. Ng, Reports on Progress in Physics 52, 699 (1989).
  * Brecher et al. (1967) C. Brecher, H. Samelson, A. Lempicki, R. Riley, and T. Peters, Physical Review 155, 178 (1967).
  * Buschow et al. (1974) K. H. J. Buschow, A. M. Van Diepen, and H. W. De Wijn, Solid State Communications 15, 903 (1974).
  * Sankar et al. (1975) S. Sankar, V. Rao, E. Segal, W. Wallace, W. Frederick, and H. Garrett, Physical Review B 11, 435 (1975).
  * Givord et al. (1979) D. Givord, J. Laforest, J. Schweizer, and F. Tasset, Journal of Applied Physics 50, 2008 (1979).
  * Tie-Song et al. (1991) Z. Tie-Song, J. Han-Min, G. Guang-Hua, H. Xiu-Feng, and C. Hong, Physical Review B 43, 8593 (1991).
  * Tils et al. (1999) P. Tils, M. Loewenhaupt, K. H. J. Buschow, and R. Eccleston, Journal of alloys and compounds 289, 28 (1999).
  * Coehoorn (1991) R. Coehoorn, Journal of Magnetism and Magnetic Materials 99, 55 (1991).
  * Daalderop et al. (1992) G. Daalderop, P. Kelly, and M. Schuurmans, Journal of Magnetism and Magnetic Materials 104, 737 (1992).
  * Steinbeck et al. (1994) L. Steinbeck, M. Richter, H. Eschrig, and U. Nitzsche, Physical Review B 49, 16289 (1994).
  * Novák and Kuriplach (1994) P. Novák and J. Kuriplach, Physical Review B 50, 2085 (1994).
  * Steinbeck et al. (1996) L. Steinbeck, M. Richter, U. Nitzsche, and H. Eschrig, Physical Review B 53, 7111 (1996).
  * Novák (1996) P. Novák, physica status solidi (b) 198, 729 (1996).
  * Ning and Brivio (2007) L. Ning and G. P. Brivio, Phys. Rev. B 75, 235126 (2007).
  * Hu et al. (2011) L. Hu, M. F. Reid, C.-K. Duan, S. Xia, and M. Yin, Journal of Physics: Condensed Matter 23, 045501 (2011).
  * Novák et al. (2013) P. Novák, K. Knížek, and J. Kuneš, Physical Review B 87, 205139 (2013).
  * Brooks et al. (1997) M. Brooks, O. Eriksson, J. Wills, and B. Johansson, Physical review letters 79, 2546 (1997).
  * Haverkort et al. (2012) M. W. Haverkort, M. Zwierzycki, and O. K. Andersen, Phys. Rev. B 85, 165113 (2012).
  * Liddle and van Slageren (2015) S. T. Liddle and J. van Slageren, Chemical Society Reviews 44, 6655 (2015).
  * Ungur and Chibotaru (2017) L. Ungur and L. F. Chibotaru, Chemistry-A European Journal 23, 3708 (2017).
  * Zhou et al. (2011) F. Zhou, V. Ozoliņš, et al., Physical Review B 83, 085106 (2011).
  * Zhou et al. (2009) F. Zhou, V. Ozoliņš, et al., Physical Review B 80, 125127 (2009).
  * Anisimov et al. (1997) V. I. Anisimov, A. I. Poteryaev, M. A. Korotin, A. O. Anokhin, and G. Kotliar, J. Phys.: Condens. Matter 9, 7359 (1997).
  * Lichtenstein and Katsnelson (1998) A. Lichtenstein and M. Katsnelson, Physical Review B 57, 6884 (1998).
  * Lebègue et al. (2005) S. Lebègue, G. Santi, A. Svane, O. Bengone, M. Katsnelson, A. Lichtenstein, and O. Eriksson, Physical Review B 72, 245102 (2005).
  * Lebègue et al. (2006) S. Lebègue, A. Svane, M. Katsnelson, A. Lichtenstein, and O. Eriksson, Journal of Physics: Condensed Matter 18, 6329 (2006).
  * Pourovskii et al. (2009) L. V. Pourovskii, K. T. Delaney, C. G. Van de Walle, N. A. Spaldin, and A. Georges, Phys. Rev. Lett. 102, 096401 (2009).
  * Locht et al. (2016) I. L. M. Locht, Y. O. Kvashnin, D. C. M. Rodrigues, M. Pereiro, A. Bergman, L. Bergqvist, A. I. Lichtenstein, M. I. Katsnelson, A. Delin, A. B. Klautau, et al., Phys. Rev. B 94, 085137 (2016).
  * Grånäs et al. (2012) O. Grånäs, I. Di Marco, P. Thunström, L. Nordström, O. Eriksson, T. Björkman, and J. Wills, Computational materials science 55, 295 (2012).
  * Stevens (1952) K. Stevens, Proceedings of the Physical Society. Section A 65, 209 (1952).
  * Wybourne and Meggers (1965) B. G. Wybourne and W. F. Meggers, _Spectroscopic properties of rare earths_ (Interscience Publishers, New York, 1965).
  * Mulak and Gajek (2000) J. Mulak and Z. Gajek, _The effective crystal field potential_ (Elsevier, 2000).
  * Newman and Ng (2007) D. J. Newman and B. Ng, _Crystal field handbook_ (Cambridge University Press, 2007).
  * Aichhorn et al. (2016) M. Aichhorn, L. Pourovskii, P. Seth, V. Vildosola, M. Zingl, O. E. Peil, X. Deng, J. Mravlje, G. J. Kraberger, C. Martins, et al., Computer Physics Communications 204, 200 (2016).
  * Parcollet et al. (2015) O. Parcollet, M. Ferrero, T. Ayral, H. Hafermann, I. Krivenko, L. Messio, and P. Seth, Computer Physics Communications 196, 398 (2015).
  * Blaha et al. (2001) P. Blaha, K. Schwarz, G. Madsen, D. Kvasnicka, and J. Luitz, _WIEN2k, An augmented Plane Wave + Local Orbitals Program for Calculating Crystal Properties_ (Techn. Universitat Wien, Austria, ISBN 3-9501031-1-2., 2001).
  * Amadon et al. (2008) B. Amadon, F. Lechermann, A. Georges, F. Jollet, T. O. Wehling, and A. I. Lichtenstein, Phys. Rev. B 77, 205112 (2008).
  * Aichhorn et al. (2009) M. Aichhorn, L. Pourovskii, V. Vildosola, M. Ferrero, O. Parcollet, T. Miyake, A. Georges, and S. Biermann, Phys. Rev. B 80, 085101 (2009).
  * Lechermann et al. (2006) F. Lechermann, A. Georges, A. Poteryaev, S. Biermann, M. Posternak, A. Yamasaki, and O. Andersen, Physical Review B 74, 125120 (2006).
  * Pourovskii et al. (2007) L. V. Pourovskii, B. Amadon, S. Biermann, and A. Georges, Phys. Rev. B 76, 235101 (2007).
  * Aichhorn et al. (2011) M. Aichhorn, L. Pourovskii, and A. Georges, Phys. Rev. B 84, 054529 (2011).
  * (55) L. Pourovskii, private communication.
  * Czyżyk and Sawatzky (1994) M. Czyżyk and G. Sawatzky, Physical Review B 49, 14211 (1994).
  * Harashima et al. (2015b) Y. Harashima, K. Terakura, H. Kino, S. Ishibashi, and T. Miyake, Proceedings of Computational Science Workshop 2014 (CSW2014) p. 011021 (2015b).
  * Nilsson et al. (2013) F. Nilsson, R. Sakuma, and F. Aryasetiawan, Phys. Rev. B 88, 125123 (2013).
  * Lang et al. (1981) J. Lang, Y. Baer, and P. Cox, Journal of Physics F: Metal Physics 11, 121 (1981).
  * Bader (1990) R. F. W. Bader, _Atoms in Molecules: a quantum theory_ (Clarendon: Oxford University Press, 1990).
  * Aryasetiawan et al. (2004) F. Aryasetiawan, M. Imada, A. Georges, G. Kotliar, S. Biermann, and A. I. Lichtenstein, Phys. Rev. B 70, 195104 (2004).
  * Mazet et al. (2013) T. Mazet, D. Malterre, M. François, C. Dallera, M. Grioni, and G. Monaco, Phys. Rev. Lett. 111, 096402 (2013).
  * Mazet et al. (2015) T. Mazet, D. Malterre, M. François, L. Eichenberger, M. Grioni, C. Dallera, and G. Monaco, Phys. Rev. B 92, 075105 (2015).
  * Cuthill et al. (1974) J. Cuthill, A. McAlister, N. Erickson, and R. Watson, AIP Conference Proceedings 18, 1039 (1974).
  * Laslo et al. (2014) A. Laslo, R. Dudric, M. Neumann, O. Isnard, M. Coldea, and V. Pop, Solid State Communications 199, 43 (2014).
  * Richter et al. (1995) M. Richter, L. Steinbeck, U. Nitzsche, P. M. Oppeneer, and H. Eschrig, Journal of Alloys and Compounds 225, 469 (1995).
  * Hummler and Fähnle (1996) K. Hummler and M. Fähnle, Phys. Rev. B 53, 3272 (1996).
  * Novak and Kuriplach (1994) P. Novak and J. Kuriplach, IEEE Transactions on Magnetics 30, 1036 (1994), ISSN 0018-9464.
  * Laforest (1981) J. Laforest, Ph.D. thesis, USM-INP, Grenoble (1981).
  * Dederichs et al. (1984) P. Dederichs, S. Blügel, R. Zeller, and H. Akai, Physical review letters 53, 2512 (1984).


# Jacob’s ladders and the $\tilde{Z}^{2}$-transformation of the orthogonal
system of trigonometric functions

Jan Moser  Department of Mathematical Analysis and Numerical Mathematics,
Comenius University, Mlynska Dolina M105, 842 48 Bratislava, SLOVAKIA
<EMAIL_ADDRESS>
###### Abstract.

It is shown in this paper that there is a continuum set of orthogonal systems
relative to the weight function $\tilde{Z}^{2}(t)$. The corresponding
integrals cannot be obtained in known theories of Balasubramanian, Heath-Brown
and Ivic.

###### Key words and phrases:

Riemann zeta-function

##  1\. The first result

### 1.1.

In this paper we obtain some new properties of the signal

(1.1) $Z(t)=e^{i\vartheta(t)}\zeta\left(\frac{1}{2}+it\right)$

that is generated by the Riemann zeta-function, where

(1.2)
$\vartheta(t)=-\frac{t}{2}\ln\pi+\text{Im}\ln\Gamma\left(\frac{1}{4}+i\frac{t}{2}\right)=\frac{t}{2}\ln\frac{t}{2\pi}-\frac{t}{2}-\frac{\pi}{8}+\mathcal{O}\left(\frac{1}{t}\right).$

Let us remind that

(1.3) $\tilde{Z}^{2}(t)=\frac{{\rm d}\varphi_{1}(t)}{{\rm d}t},\
\varphi_{1}(t)=\frac{1}{2}\varphi(t)$

where

(1.4)
$\tilde{Z}^{2}(t)=\frac{Z^{2}(t)}{2\Phi^{\prime}_{\varphi}[\varphi(t)]}=\frac{Z^{2}(t)}{\left\\{1+\mathcal{O}\left(\frac{\ln\ln
t}{\ln t}\right)\right\\}\ln t}$

(see [12], (5.1)-(5.3)) and $\varphi_{1}(T),\ T\geq T_{0}[\varphi_{1}]$ is the
Jacob’s ladder.

### 1.2.

It is known that the system of trigonometric functions

(1.5)
$\left\\{1,\cos\left(\frac{\pi}{l}t\right),\sin\left(\frac{\pi}{l}t\right),\dots,\cos\left(\frac{\pi}{l}nt\right),\sin\left(\frac{\pi}{l}nt\right),\dots\right\\}$

is the orthogonal system on the segment $[0,2l]$. In this direction the
following theorem holds true.

######  Theorem 1.

Let $\mathcal{J}(2l)=\varphi_{1}\\{\mathring{\mathcal{J}}(2l)\\}$, where

Then the system of functions

(1.7)
$\left\\{1,\cos\left(\frac{\pi}{l}\varphi_{1}(t)\right),\sin\left(\frac{\pi}{l}\varphi_{1}(t)\right),\dots,\cos\left(\frac{\pi}{l}n\varphi_{1}(t)\right),\sin\left(\frac{\pi}{l}n\varphi_{1}(t)\right),\dots\right\\}$

is the orthogonal system on $\mathring{\mathcal{J}}(2l)$ with respect to the
weight function $\tilde{Z}^{2}(t)$, i.e. the following new system of integrals

(1.8)
$\begin{split}&\int_{\mathring{\mathcal{J}}(2l)}\cos\left(\frac{\pi}{l}m\varphi_{1}(t)\right)\cos\left(\frac{\pi}{l}n\varphi_{1}(t)\right)\tilde{Z}^{2}(t){\rm
d}t=\left\\{\begin{array}[]{rcl}0&,&m\not=n,\\\ l&,&m=n,\end{array}\right.\\\
&\int_{\mathring{\mathcal{J}}(2l)}\sin\left(\frac{\pi}{l}m\varphi_{1}(t)\right)\sin\left(\frac{\pi}{l}n\varphi_{1}(t)\right)\tilde{Z}^{2}(t){\rm
d}t=\left\\{\begin{array}[]{rcl}0&,&m\not=n,\\\ l&,&m=n,\end{array}\right.\\\
&\int_{\mathring{\mathcal{J}}(2l)}\sin\left(\frac{\pi}{l}m\varphi_{1}(t)\right)\cos\left(\frac{\pi}{l}n\varphi_{1}(t)\right)\tilde{Z}^{2}(t){\rm
d}t=0,\\\
&\int_{\mathring{\mathcal{J}}(2l)}\cos\left(\frac{\pi}{l}n\varphi_{1}(t)\right)\tilde{Z}^{2}(t){\rm
d}t=0,\\\
&\int_{\mathring{\mathcal{J}}(2l)}\sin\left(\frac{\pi}{l}n\varphi_{1}(t)\right)\tilde{Z}^{2}(t){\rm
d}t=0\end{split}$

for all $m,n\in\mathbb{N}$ is obtained, where

(A) $t-\varphi_{1}(t)\sim(1-c)\pi(t),$ (B) $2l(K+1)<\widering{2lK},$ (C)
$\rho\\{\mathcal{J}(2l);\mathring{\mathcal{J}}(2l)\\}\sim(1-c)\pi(t)\to\infty,$

as $K\to\infty$, and $\rho$ denotes the distance of the corresponding
segments, $c$ is the Euler constant and $\pi(t)$ is the prime-counting
function.

######  Remark 1.

Theorem 1 gives the contact point between the functions
$\zeta\left(\frac{1}{2}+it\right),\ \pi(t),\ \varphi_{1}(t)$ and the
orthogonal system of trigonometric functions.

######  Remark 2.

It is clear that the formulae (1.8) - for the modulated function
$\tilde{Z}^{2}(t)$ \- cannot be obtained in the known theories of
Balasubramanian, Heath-Brown and Ivic (comp. [1]).

This paper is a continuation of the series [2]-[15].

##  2\. New method of the quantization of the Hardy-Littlewood integral (a
special case)

### 2.1.

We obtain from the first two formulae in (1.8)

(2.1)
$\begin{split}&\int_{\mathring{\mathcal{J}}(2l)}\cos^{2}\left(\frac{\pi}{l}m\varphi_{1}(t)\right)\tilde{Z}^{2}(t){\rm
d}t=\frac{1}{2}|\mathcal{J}(2l)|,\\\
&\int_{\mathring{\mathcal{J}}(2l)}\sin^{2}\left(\frac{\pi}{l}m\varphi_{1}(t)\right)\tilde{Z}^{2}(t){\rm
d}t=\frac{1}{2}|\mathcal{J}(2l)|\end{split}$

for all $m\in\mathbb{N}$. Next, from (2.1) we obtain

######  Corollary 1.

(2.2) $\int_{\mathring{\mathcal{J}}(2l)}\tilde{Z}^{2}(t){\rm
d}t=|\mathcal{J}(2l)|;\ |\mathcal{J}(2l)|=2l.$

### 2.2.

Let us consider now the problem concerning the solid of revolution
corresponding to the graph of the function (comp. [5])

$\tilde{Z}(t),\ t\in[\widering{2lK},+\infty),\ 2lK>T_{0}[\varphi_{1}].$

######  Problem.

To divide this solid of revolution on parts of equal volumes.

From (2.2) we obtain the resolution of this problem.

######  Corollary 2.

Since

(a)
$[\widering{2lK},+\infty)=\bigcup_{r=1}^{\infty}\mathring{\mathcal{J}}(2l,r),\
\mathring{\mathcal{J}}(2l,r)=[\widering{2l(K+r-1)},\widering{2l(K+r)}],$ (b)
$\pi\int_{\mathring{\mathcal{J}}(2l,r)}\tilde{Z}^{2}(t){\rm d}t=2\pi l,\
r=1,2,3,\dots\ ,$

it follows that the sequence of points

$\\{\widering{2l(K+r-1)}\\}_{r=2}^{+\infty}$

is the resolution to the Problem for arbitrary fixed $2l\in(0,T/\ln T]$.

##  3\. Generalization of the formula (2.2)

### 3.1.

The following theorem holds true.

######  Theorem 2.

Let

$\mathcal{J}(T,U)=[T,T+U],\
J(T,U)=\varphi_{1}\\{\mathring{\mathcal{J}}(T,U)\\};\
\mathring{\mathcal{J}}(T,U)=[\mathring{T},\widering{T+U}].$

Then

(3.1) $\int_{\mathring{\mathcal{J}}(T,U)}\tilde{Z}^{2}(t){\rm
d}t=|\mathcal{J}(T,U)|=U,$

for every $T\geq T_{0}[\varphi_{1}],\ U\in(0,T/\ln T]$.

######  Remark 3.

From (3.1) the general method for quantization of the Hardy-Littlewood
integral follows (comp. Corollary 2: $2lK\to\forall\ T\geq
T_{0}[\varphi_{1}],\ \mathcal{J}(2l)\to\mathcal{J}(T,U)$).

Next, we obtain, using the mean-value theorem in (3.1)

######  Corollary 3.

(3.2)
$\tilde{Z}^{2}(\xi)=\frac{|\mathcal{J}(T,U)|}{|\mathring{\mathcal{J}}(T,U)|},\
\xi\in\xi(\mathring{T},\widering{T+U}),\ \tilde{Z}(\xi)\not=0,$

i.e.

$\tilde{Z}^{2}(\xi):1=|\mathcal{J}(T,U)|:|\mathring{\mathcal{J}}(T,U)|.$

### 3.2.

Let $\\{[T^{\prime},T^{\prime}+1]\\}$ stands for the continuum set of segments
$[T^{\prime},T^{\prime}+1]\subset[T,T+T/\ln T]$. Since

$\frac{1}{|\mathring{\mathcal{J}}(T^{\prime},1)|}=\tilde{Z}^{2}(\xi),\
\xi=\xi(T^{\prime})\in(\mathring{T}^{\prime},\widering{T^{\prime}+1})$

then by the Riemann-Siegel formula

$Z(t)=2\sum_{n\leq\sqrt{\frac{t}{2\pi}}}\frac{1}{\sqrt{n}}\cos\\{\vartheta(t)-t\ln
n\\}+\mathcal{O}(t^{-1/4})$

we obtain (see (1.4))

######  Corollary 4.

(3.3)
$\frac{1}{\sqrt{|\mathring{\mathcal{J}}(T^{\prime},1)|}}\sim\frac{2}{\sqrt{\ln\xi}}\left|\sum_{n\leq\sqrt{\frac{\xi}{2\pi}}}\frac{1}{\sqrt{n}}\cos\\{\vartheta(\xi)-\xi\ln
n\\}+\mathcal{O}(\xi^{-1/4})\right|$

where $\xi=\xi(T^{\prime})$.

######  Remark 4.

The formula (3.3) describes the complicated oscillations of the value
$|\mathring{\mathcal{J}}(T^{\prime},1)|$ generated by the nonlinear
transformation
$\mathcal{J}(T^{\prime},1)=\varphi_{1}\\{\mathring{\mathcal{J}}(T^{\prime},1)\\}$.

##  4\. Proof of Theorems 1 and 2

### 4.1.

Let us remind that the following lemma is true (see [12], (5.1)-(5.3))

######  Lemma.

For every integrable function (in the Lebesgue sense) $f(x),\
x\in[\varphi_{1}(T),\varphi_{1}(T+U)]$ the following is true

(4.1) $\int_{T}^{T+U}f[\varphi_{1}(t)]\tilde{Z}^{2}(t){\rm
d}t=\int_{\varphi_{1}(T)}^{\varphi_{1}(T+U)}f(x){\rm d}x,\ U\in(0,T/\ln T],$

where $t-\varphi_{1}(t)\sim(1-c)\pi(t)$.

######  Remark 5.

The formula (4.1) is true also in the case when the integral on the right-hand
side of eq. (4.1) is convergent but not absolutely (in the Riemann sense).

### 4.2.

If $\varphi_{1}\\{[\mathring{T},\widering{T+U}]\\}=[T,T+U]$ then we obtain
from (4.1) the following formula

(4.2)
$\int_{\mathring{T}}^{\widering{T+U}}f[\varphi_{1}(t)]\tilde{Z}^{2}(t){\rm
d}t=\int_{T}^{T+U}f(x){\rm d}x,\ U\in(0,T/\ln T].$

Next, in the case $[T,T+U]=[2lK,2lK+2l]=\mathcal{J}(2l)$, we have

(4.3) $\int_{\mathcal{J}(2l)}F(t){\rm d}t=\int_{0}^{2l}F(t){\rm d}t$

for every (integrable) $2l$-periodic function $F(t)$. Then from the known
formulae

$\begin{split}&\int_{\mathcal{J}(2l)}\cos\left(\frac{\pi}{l}m\varphi_{1}(t)\right)\cos\left(\frac{\pi}{l}n\varphi_{1}(t)\right){\rm
d}t=\left\\{\begin{array}[]{rcl}0&,&m\not=n,\\\ l&,&m=n,\end{array}\right.\\\
&\int_{\mathcal{J}(2l)}\sin\left(\frac{\pi}{l}m\varphi_{1}(t)\right)\sin\left(\frac{\pi}{l}n\varphi_{1}(t)\right){\rm
d}t=\left\\{\begin{array}[]{rcl}0&,&m\not=n,\\\ l&,&m=n,\end{array}\right.\\\
&\int_{\mathcal{J}(2l)}\sin\left(\frac{\pi}{l}m\varphi_{1}(t)\right)\cos\left(\frac{\pi}{l}n\varphi_{1}(t)\right){\rm
d}t=0,\\\
&\int_{\mathcal{J}(2l)}\cos\left(\frac{\pi}{l}n\varphi_{1}(t)\right){\rm
d}t=0,\quad\int_{\mathcal{J}(2l)}\sin\left(\frac{\pi}{l}n\varphi_{1}(t)\right){\rm
d}t=0,\quad m,n\in\mathbb{N},\end{split}$

by the $\tilde{Z}^{2}$-transformation (see (4.2), (4.3);
$[\mathring{T},\widering{T+U}]=\mathring{\mathcal{J}}(2l)$) the formulae (1.8)
follow. The properties (B), (C) in Theorem 1 are identical with [13], (A1),
(B1).

### 4.3.

The formula (3.1) follows from (4.2) in the case $f(x)\equiv 1$.

##  5\. Another type of the orthogonal systems

It follows from (4.2) that the continuum set $\mathcal{S}(T,2l)$ of the
systems

$\begin{split}&\left\\{|\tilde{Z}(t)|,|\tilde{Z}(t)|\cos\left(\frac{\pi}{l}(\varphi_{1}(t)-T)\right),|\tilde{Z}(t)|\sin\left(\frac{\pi}{l}(\varphi_{1}(t)-T)\right),\dots,\right.\\\
&\left.|\tilde{Z}(t)|\cos\left(\frac{\pi}{l}n(\varphi_{1}(t)-T)\right),|\tilde{Z}(t)|\sin\left(\frac{\pi}{l}n(\varphi_{1}(t)-T)\right),\dots\right\\},\\\
&t\in[\mathring{T},\widering{T+2l}]\end{split}$

for all

$T\geq T_{0}[\varphi_{1}],\ 2l\in(0,T/\ln T]$

is the set of orthogonal systems on $[\mathring{T},\widering{T+2l}]$.

######  Remark 6.

Let us call the elements of the system $\mathcal{S}(T,2l)$ for fixed
$2l\in(0,T/\ln T]$ and for all $T\geq T_{0}[\varphi_{1}]$ as _the clones_ of
the known orthogonal trigonometric system

$\left\\{1,\cos\left(\frac{\pi}{l}t\right),\sin\left(\frac{\pi}{l}t\right),\dots,\cos\left(\frac{\pi}{l}nt\right),\sin\left(\frac{\pi}{l}nt\right),\dots\right\\},\
t\in[0,2l].$

I would like to thank Michal Demetrian for helping me with the electronic
version of this work.

## References

  * [1] A. Ivic, ‘The Riemann zeta-function‘, A Willey-Interscience Pub., New York, 1985.
  * [2] J. Moser, ‘Jacob’s ladders and the almost exact asymptotic representation of the Hardy-Littlewood integral’, (2008), arXiv:0901.3973.
  * [3] J. Moser, ‘Jacob’s ladders and the tangent law for short parts of the Hardy-Littlewood integral’, (2009), arXiv:0906.0659.
  * [4] J. Moser, ‘Jacob’s ladders and the multiplicative asymptotic formula for short and microscopic parts of the Hardy-Littlewood integral’, (2009), arXiv:0907.0301.
  * [5] J. Moser, ‘Jacob’s ladders and the quantization of the Hardy-Littlewood integral’, (2009), arXiv:0909.3928.
  * [6] J. Moser, ‘Jacob’s ladders and the first asymptotic formula for the expression of the sixth order $|\zeta(1/2+i\varphi(t)/2)|^{4}|\zeta(1/2+it)|^{2}$’, (2009), arXiv:0911.1246.
  * [7] J. Moser, ‘Jacob’s ladders and the first asymptotic formula for the expression of the fifth order $Z[\varphi(t)/2+\rho_{1}]Z[\varphi(t)/2+\rho_{2}]Z[\varphi(t)/2+\rho_{3}]\hat{Z}^{2}(t)$ for the collection of disconnected sets‘, (2009), arXiv:0912.0130.
  * [8] J. Moser, ‘Jacob’s ladders, the iterations of Jacob’s ladder $\varphi_{1}^{k}(t)$ and asymptotic formulae for the integrals of the products $Z^{2}[\varphi^{n}_{1}(t)]Z^{2}[\varphi^{n-1}(t)]\cdots Z^{2}[\varphi^{0}_{1}(t)]$ for arbitrary fixed $n\in\mathbb{N}$‘ (2010), arXiv:1001.1632.
  * [9] J. Moser, ‘Jacob’s ladders and the asymptotic formula for the integral of the eight order expression $|\zeta(1/2+i\varphi_{2}(t))|^{4}|\zeta(1/2+it)|^{4}$‘, (2010), arXiv:1001.2114.
  * [10] J. Moser, ‘Jacob’s ladders and the asymptotically approximate solutions of a nonlinear diophantine equation‘, (2010), arXiv: 1001.3019.
  * [11] J. Moser, ‘Jacob’s ladders and the asymptotic formula for short and microscopic parts of the Hardy-Littlewood integral of the function $|\zeta(1/2+it)|^{4}$‘, (2010), arXiv:1001.4007.
  * [12] J. Moser, ‘Jacob’s ladders and the nonlocal interaction of the function $|\zeta(1/2+it)|$ with $\arg\zeta(1/2+it)$ on the distance $\sim(1-c)\pi(t)$‘, (2010), arXiv: 1004.0169.
  * [13] J. Moser, ‘Jacob’s ladders and the $\tilde{Z}^{2}$ \- transformation of polynomials in $\ln\varphi_{1}(t)$‘, (2010), arXiv: 1005.2052.
  * [14] J. Moser, ‘Jacob’s ladders and the oscillations of the function $|\zeta\left(\frac{1}{2}+it\right)|^{2}$ around the main part of its mean-value; law of the almost exact equality of the corresponding areas‘, (2010), arXiv: 1006.4316
  * [15] J. Moser, ‘Jacob’s ladders and the nonlocal interaction of the function $Z(t)$ with the function $\tilde{Z}^{2}(t)$ on the distance $\sim(1-c)\pi(t)$ for a collection of disconneted sets‘, (2010), arXiv: 1006.5158