Rare Earth Monopnictides and Monochalcogenides from First Principles: Towards an Electronic Phase Diagram of Strongly Correlated Materials

We present results of an ab-initio study of the electronic structure of 140 rare earth compounds. Specifically we predict an electronic phase diagram of the entire range of rare earth monopnictides and monochalcogenides, composed of metallic, semiconducting and heavy fermion-like regions, and exhibiting valency transitions brought about by a complex interplay between ligand chemistry and lanthanide contraction. The calculations exploit the combined effect of a first-principles methodology, which can adequately describe the dual character of electrons, itinerant vs. localized, and high throughput computing made possible by the increasing available computational power. Our findings, including the predicted"intermediate valent"compounds SmO and TmSe, are in overall excellent agreement with the available experimental data. The accuracy of the approach, proven e.g. through the lattice parameters calculated to within 1.5% of the experimental values, and its ability to describe localization phenomena in solids, makes it a competitive atomistic simulation approach in the search for and design of new materials with specific physical properties and possible technological applications.


I. INTRODUCTION
Increased computer power, together with reliable electronic structure codes, have opened up the possibility of exploring electronic properties of huge numbers of materials. 1 A systematic first principles based theoretical study of entire families of compounds has the potential of exploring new trends in their physical properties or discovering materials with improved characteristics of possible interest for technological applications.
The ab-initio determination of materials properties requires electronic structure calculations that are accurate and parameter-free. With the advent of density functional theory (DFT), 2 quantum mechanical calculations of the electronic structure of solids were put on a rigorous theoretical foundation. Nevertheless, the lack of knowledge of the exact exchange-correlation energy functional constitutes a substantial obstacle regarding the predictive power of the calculations, as different mechanisms are at work depending on the relative importance of band formation energy on one hand and on-site electron-electron interaction and electron localization on the other. Standard electronic structure calculations, based on the local spin density approximation (LSDA) 3,4 to density functional theory, work very well for conventional materials like simple metals and their alloys, characterized by delocalized electron states and large band formation energies. However, for systems where localized electron states occur, such as the 4f 's in rare earth compounds, the LSDA suffers from a sizeable self-interaction (SI) error, 5 brought about by the local approximation to the exchange-correlation energy functional. Consequently, different assumptions and/or parameters need to be invoked to describe the electronic structure of strongly correlated electron materials. 6 The self-interaction corrected 7 local-spin-density (SIC-LSD) energy functional is free of the SI error and fulfils the requirement of being parameter free and of sufficient accuracy to study trends in properties of entire families of materials containing both localized and itinerant electrons. [8][9][10] Here we apply the SIC-LSD to provide global understanding of the ground state electronic structure of rare earth monopnictides and monochalcogenides and, through that, also demonstrate the predictive power and accuracy of this approach. Specifically, the present study covers 140 distinct compounds, RX, with R referring to the rare earth atoms Ce, Pr, Nd, Pm, Sm, Eu, Gd, Tb, Dy, Ho, Er, Tm, Yb, Lu, and X standing for the pnictide atoms N, P, As, Sb, Bi and the chalcogenide atoms O, S, Se, Te, Po. Most of these compounds are known to exist in nature, crystallizing predominantly in the rocksalt structure. Due to the intricate dual nature of the 4f -electrons, as a function of chemical environment these compounds display an extraordinarily wide range of electronic, magnetic, optical, and magneto-optical properties with, among other things, potential practical applications in the field of spintronics. 11,12 A considerable number of theoretical studies on the RX compounds have resulted in greatly improved insight into their properties. [13][14][15][16] However those calculations have often been guided by experiment, relying on parameters to reproduce the observed properties. The goal of the present systematic study is to establish a paramater-free unified theoretical picture of the RX electronic properties and possibly provide directions for further experimental investigations. With future applications in mind, we show how large scale computational infrastructures can potentially be exploited in the search for the novel rare earth related materials that play an increasingly prominent role in innovative technologies and the quest for renewable energies. 17

II. THEORETICAL BACKGROUND
In the LSD approximation the total energy functional is decomposed as whereT is the kinetic energy operator, and V ext [n] represents the interaction with an external potential. The electron-electron interaction is given as the sum of the Hartree term U[n] and the local exchange-correlation term, E LSD xc [n ↑ , n ↓ ]. In the Hartree-Fock approximation, the expectation value of the electron repulsion operator in some f n multiplet state Ψ n (LSJ) may be written as 18,19 where Ψ n (LSJ) is in general given by a linear combination of n × n Slater determinants.
Since 1 r ij commutes with total L, S, and J, the above matrix elements are equal for M J quantum numbers belonging to one (LSJ) family. E i are the Racah parameters, 20 which are related to the Slater two-electron integrals (defined later). The coefficients α i depend on the particular LSJ-configuration chosen.
The quantity E 0 is the Coulomb repulsion between two f -electrons, and the Coulomb energy term in the LSD total energy functional (1) corresponds to the 1 2 E 0 n 2 part of the first term in Eq. (2). The term, − 1 2 E 0 n, in Eq. (2), cancels an equivalent contribution in the Hartree energy term, representing the interaction of an electron with itself. A similar cancellation does not occur in the LSD functional, which therefore carries an unphysical selfinteraction. 7 This interaction, which tends to be insignificant for extended band state, may lead to uncontrollable errors in the description of atomic-like localized states, for example the f -electrons in the rare earths. Following the suggestion by Perdew and Zunger 7 , in the SIC-LSD methodology, the LSD functional (1) is corrected for this spurious self-interaction by adding an explicit energy contribution for an electron to localize. The resulting, orbital dependent, SIC-LSD total energy functional has the form 8,10 where The self-interaction energy (4) consists of the self-Coulomb and self-exchange-correlation energies of the occupied orbitals ψ α with charge density n α and spin densityn α = (n ↑ α , n ↓ α ). For itinerant states, the self-interaction δ SIC α vanishes identically, while for localized (atomiclike) states δ SIC α may be appreciable. Thus, the self-interaction correction constitutes a negative energy contribution gained by an electron upon localization, which competes with the band formation energy gained by the electron if allowed to delocalize and hybridize with the available conduction states. Different localized/delocalized configurations are realized by assuming different numbers of localized states -here f -states on rare earth atom sites. The SIC-LSD approach is fully ab-initio, as both localized and delocalized states are expanded in the same set of basis functions, and are thus treated on an equal footing. If no localized states are assumed, E SIC−LSD coincides with the conventional LSD functional, i.e., the Kohn-Sham minimum of the E LSD functional is also a local minimum of E SIC−LSD .
The spin-orbit interaction (5) couples the band Hamiltonian for the spin-up and spin-down channels, i.e. a double secular problem must be solved. The spin-orbit parameter, in atomic Rydberg units, is calculated from the self-consistent potential V .
Of the remaining terms in Eq. (2), the E 1 term accounts for Hund's first rule. The E 2 term only contributes to the level spacing of excited multiplets and will not be relevant in a functional for the ground state energy. The respectively. It is an atomic effect which results in an increased stability at 1/4 and 3/4 filling of the 4f -shell, and this second Hunds's rule effect is often referred to as the tetrad effect (TE). 19 Although, both LSD and SIC-LSD take into account Hund's first and third rules, with respectively the exchange interaction and spin-orbit coupling included in the total energy functional, Hund's second rule is not defined for the homogeneous electron gas, underlying LDA. 21 To account for it in our total energy calculations, we add (a posteriori) to the SIC-LSD total energy functional the relevant correction The E 3 parameter is equivalent to Racah's B-parameter, 20 and is given in terms of the reduced Slater integrals F k as: where F k = F k /D k . Here D k are numerical constants (D k =225, 1089, 7361.64 for k=2, 4, 6) and the Slater integral F k is defined through Here, φ 4f is the f-partial wave as calculated in the self-consistent crystal potential, and r < (r > ) denotes the smaller (larger) of the variables r 1 and r 2 . In reality, correlation effects in the solid environment tend to reduce the multiplet energy level splittings. Hence our calculated TE values are on average 15-20% larger than their experimental counterparts. 19

III. CALCULATIONAL DETAILS
Given the total energy functional E SIC−LSD , the computational procedure is as for the LSD case, i.e. minimization is accomplished by iteration until self-consistency. In the present work, the electron wavefunctions are expanded in the linear-muffin-tin-orbital (LMTO) basis functions. 22 The atomic spheres approximation (ASA) is used, whereby the crystal volume is divided into slightly overlapping atom-centered spheres of a total volume equal to the actual volume. A known shortcoming of the ASA is that different crystal structures have different degrees of overlap of the ASA spheres resulting in substantial relative errors in the evaluation of the total energy. While this inhibits the comparison of energies of different crystal structures, when comparing the energies of different localization/delocalization scenarios within the same crystal structure the ASA error is of minor influence. To improve the packing of the structure empty spheres have been introduced on high symmetry interstitial sites. Two uncoupled energy panels have been considered when constructing the LMTO's to ensure an accurate description of semicore states. The LS-coupling scheme is adopted for the localized f -states by starting the iterations with Wannier states of appropriate symmetry. During iteration to self-consistency the symmetry of the Wannier states may change, however grossly retaining their overall characteristics due to the fact that the energy scale of spin-orbit interaction is smaller than that of exchange, but larger than that of crystal field for the f -states. The rocksalt structure with the magnetic moments arranged ferromagnetically was assumed in all the calculations. To conduct the present large scale systematic study of rare earth compounds, in addition to a reliable theory and computer codes, we have made an extensive use of grid computing and automated data/metadata management tools. 23,24

IV. TOTAL ENERGIES AND VALENCY TRANSITIONS
In the SIC-LSD methodology, both the itinerant and localized limits are described by the same energy functional, which enables us to predict the ground state of correlated electron materials from total energy considerations. 8,9 Hence, one can investigate localization phenomena in solids 6 by realizing and studying different localization/delocalization scenarios, giving rise to different valency configurations, and through total energy minimization de- obtained by subtracting from the atomic number (Z) the sum of core (N core ), and localized (N SIC ) electrons. In this paper we will be using two interchangeable nomenclatures, f n and R m+ , to describe the configuration of the rare earth ion, implying n = N SIC and m = N val , respectively. The total number of f -electrons may be larger than n, because, in addition to the n localized f -states, the band states also contribute to the total f -count on a given ion. Note that our calculated valencies, refer to the number of rare earth electrons that contribute to bonding, and thus do not necessarily coincide with the nominal (ionic) valency of a compound.
In our study of the rare earth monopnictides and monochalcogenides, for each given compound the self-consistent calculations involve the total energy minimization as a function of with respect to f -electrons, the fully localized limit. Configurations with one or two delocalized (bonding) f -electrons are respectively referred to as trivalent and tetravalent. The LDA scenario, which treats all the ligand and rare earth valence electrons, including the 11 f -electrons, as band electrons, is found energetically unfavourable by ∼ 1 Ry, as compared to any of the other scenarios, and is therefore not shown in Fig. 1a. Altogether at least some forty self-consistent calculations are required to establish the global energy minimum for a compound like HoSe which, as can be seen in Fig. 1a, is obtained in the trivalent scenario.
From the calculated local energy minima we can evaluate the energy differences between respectively the divalent and trivalent configurations, ∆E II−III , and the tetravalent and trivalent configurations, ∆E IV −III , both of which are a measure of stability of the ground state configuration. From Fig. 1a Yb. We notice that the total energy trends, observed for the selenide series in Fig. 1b, are reproduced in all the ligands, but overall, in the chalcogenides the divalent configuration becomes increasingly more favourable as a function of ligand size, whilst in the pnictides the trivalent configuration remains the ground state, even for the "late" rare earth Eu and Yb (with the exception of EuBi). The trivalent to divalent localization transition thus mainly pertains to the rare earth chalcogenides. In Fig. 2b, the same total energy data are depicted  never actually becomes the true ground state. The trivalent-tetravalent energy difference, ∆E IV −III , is the more relevant energy scale for the Ce compounds for example, as the divalent configuration is always very unfavourable, and in most cases could not even be converged. In the case of CeN, and to some degree in CeO, the tetravalent and trivalent scenarios are actually close to energetically degenerate. Since the divalent scenario is so highly unfavourable, the Ce pnictides and chalcogenides, as well as all the other rare earth nitrides, are not shown in Fig. 2. For the LuX compounds ∆E II−III is not adequate, as a divalent scenario would here imply 15 localized f -electrons.
As shown in Fig. 3   (Tm 3+ ≡ Tm(f 12 )). Compared to the LSD scenario the trivalent configuration is found to be energetically more favourable by ∼ 700 mRy for TbP and ∼ 1200 mRy for TmTe.
In TbP the seven majority spin states , as well as one of the minority spin states, are no longer available for band formation and hybridization and have vanished from the band picture in Fig. 4b. Similarly in the DOS for TmTe (Fig. 4b'), hand, the f -electron localization increases due to the lanthanide contraction, as we move through the rare earth series from Ce to Eu, and again from Gd to Yb. In terms of DOS, these two trends manifest themselves respectively via ligand p-bands moving up in energy (towards lower binding energies), and unoccupied f -states moving towards higher binding energies. In Fig. 5, the combined effect of these two trends is depicted in the two sets of DOSs respectively along the diagonals of the pnictide ((A) → (C)) and chalcogenide compounds ((D) → (F)). As the f -states move towards lower energy they get pinned to the Fermi level and start filling up. Ultimately the degree of f -band filling becomes such that the gain in energy associated with localizing the given f -state becomes more important than the corresponding loss in band formation energy, and the divalent scenario becomes the ground state (scenarios (C') and (F')). 30 As seen in Fig. 5, owing to the additional electron, in the chalcogenides the filling of f -states sets in much earlier than in the pnictides, which explains why the trivalent scenario is relatively more dominant in the pnictides.

A. Lattice parameters
For each given compound, the calculated global energy minimum determines the equilibrium lattice parameter as well as the ground state. The calculated lattice parameters for the rare earth sulphides and the ytterbium pnictides and chalcogenides are compared to about 5 %. Therefore, a comparison to the experimental lattice parameters is expected to provide a good indication as to whether the calculations predict the correct ground state and valency configuration. In Fig. 7, the experimental equilibrium lattice parameters for all the rare earth monopnictides and monochalcogenides 31 are plotted along the x-axis, whilst our calculated equilibrium counterparts are plotted along the y-axis. In this representation, the proximity of the data points to the x = y line is a measure of the agreement between theory and experiment. We find that for the large majority of the compounds the calculated equilibrium lattice parameters lie within ∼1.5% of the experimental values, i.e. well below the 5% that would indicate a wrongly predicted ground state configuration. The only exceptions are CeN, DyTe, and three of the polonide compounds. With respect to the latter compounds, polonium is difficult to handle, and it is unclear to what degree the experimental results refer to stoichiometric samples. 32 It is well established from rare earth halide data 33 as well as XPS measurements on the rare earth elements 26 that the divalent configuration becomes more competitive in Dy compounds due to the 3/4 filling of the f -shell. The inclusion of the TE leads to a considerable reduction of E II − E III for DyX compounds, as can be seen on the example of DySe in Fig. 1b, and quite generally gives rise to the negative dip around Dy in Fig. 2b. This trend reproduces overall very well the experimentally observed data, 27 however we wrongly predict a Dy 2+ ground state in DyTe. This seems to indicate that in some cases we tend to overestimate the tetrad effect's influence as a result of approximating the respective Slater physical properties of these materials.

B. Electronic Phase Diagram
Here we concentrate on establishing global trends in the ground state electronic structure which we base on the calculated DOS at the Fermi energy, E F , as providing a guideline for predicting materials with specific properties when we systematically scan the full manifold of rare earth compounds. The corresponding results are summarized in the form of an "electronic phase diagram" presented in Fig. 8, where based on the calculated DOS we distinguish between semiconductors, semi-metals, metals, and heavy fermion-like compounds.
Looking at the phase diagram, we see that our calculations predict a considerable number of compounds to be semiconductors. In particular, with respect to the pnictides, NdN, The dominant part of the electronic phase diagram in Fig. 8 consists of metallic compounds. The indicated differences in conductivities can be straightforwardly related to the trends in electronic structure observed when moving along the diagonals in Fig 5. In the pnictides, the closing of the semiconducting gap due to the onset of pnictogen p -rare earth sd overlap gives at first rise to a vanishingly small DOS, resulting in the semi-metals CeX to PrX and GdX to DyX (light blue areas). As the overlap gradually increases (due to decreasing electronegativity and/or increasing lanthanide contraction), the behaviour becomes fully metallic (dark blue areas), and eventually in the "late" rare earths, with more and more of the f -states getting occupied, heavy fermion-like behaviour sets in (yellow-red areas). The rare earth chalcogenides display less variety in terms of metallic behaviour. Here we start out with a partially occupied conduction band, and semi-metallic behaviour is not observed.
Instead, as can be seen from Fig. 8, the rare earth systems, CeX to PmX, and GdX to ErX (apart from the polonides and a couple of tellurides) are metallic. Compared to the pnictides, in the late rare earth chalcogenides, the f -states start filling even before they have moved far enough towards lower binding energies to achieve significant hybridization.
Instead of forming a heavy metal, a relatively pure f -state pins the Fermi level (scenarios (E) and (F) in Fig. 5) and the divalent semiconducting scenario (scenario (F') in Fig. 5) correspondingly becomes energetically more favourable.
The agreement between observed and predicted properties is overall very good, except maybe for the calculated heavy fermion pnictides. Thus the "late" rare earths, SmX, and YbX, are described as low carrier semi-metallic compounds, seemingly at odds with the large DOS at the Fermi level that is predicted in the SIC-LSD calculations. The experimentally observed large specific heat coefficient in these systems has been interpreted as heavy fermion behaviour by some authors, 42 whilst others have explained it in terms of the Kondo effect. 43,44 Since our calculations are based on the effective one electron picture, they can not adequately describe the many-body physics of these compounds. However the large peak at the Fermi level that is calculated in the trivalent configuration indicates strong pf -mixing, i.e. a prerequiste for the hybridization between the narrow f -peak and the sea of conduction electrons that underpins the Kondo effect. EuSb, YbBi, EuBi (as well as EuAs) do not actually crystallize in the NaCl structure investigated here, 45 i.e., the very same pnictides that we find to be divalent or close to divalent do not occur in nature. As indicated schematically in Fig. 5, for the pnictides, localization of an additional f -electron ((C) → (C') transition) implies the removal of some bonding p-states. This, however, would destabilize the underlying NaCl crystal structure and, as a result of this competition between band formation and correlations, it seems that these materials prefer to distort and crystallize in a different structure altogether.
In Fig. 8  In summary, this effective one-electron SIC-LSD methodology accurately predicts the phase transitions associated with the f -electron localization/delocalization. However, we should note here that because of the many-body effects, beyond those included in the SIC exchange-correlation energy functional, the representation of the calculated densities of states in terms of the described above phase diagram is most likely not the full picture of the underlying electronic structure. This is in particlular true in the 'IV' and heavy fermion areas of the phase diagram, where those 'left out' many-body interactions will introduce additional refinements into the physical trends based on SIC-LSD. Therefore, compounds that are situated in those areas provide the greater challenge and should be further investigated by methodologies such as the dynamical mean field theory (DMFT) 53 or GW. 54

VII. CONCLUSIONS AND OUTLOOK
We have shown that the orbital dependent DFT approach, based on the LSD energy functional, appended by two important corrections, namely the self-interaction correction and the tetrad effect, can provide a global picture of the electronic structure of the rare earth monopnictides and monochalcogenides, expressed through the uncovered trends in their behaviour. We believe that this SIC-LSD based study provides an important understanding of the physical properties of these compounds. Although for a given compound the DMFT+LDA or GW approaches can give added insight into the details of the electronic structure, they can not at the present time be employed for large scale predictive studies given that they either depend on parameters and/or are very time consuming.
Our study has given rise to a unified picture of the ground state valency and valency transitions across the entire range of rare earth monopnictides and monochalcogenides, starting from tetravalency, through trivalency to divalency, as one moves from the early to late rare earths, in correlation with the ligands. We have also discovered the important links between the underlying electronic structure, valency and physical properties of these compounds, ranging from semiconducting to semi-metallic, metallic and even heavy-fermion-like behaviour.
The present study has also established that if the energy differences between the energetically relevant valency scenarios are small, then fluctuations will play an important role, giving rise to intermediate ground states. This effect can be described within the SIC-LSD formalism by mapping the valence (and also spin) fluctuations on to disorder in the spirit of the Hubbard III approach, 55 which we intend to apply in future to the "intermediate" compounds identified in the present study.