Y2O3:Eu and the Mössbauer isomer shift coefficient of Eu compounds from ab-initio simulations

We report on a full potential density functional theory characterization of Y2O3 upon Eu doping on the two inequivalent crystallographic sites 24d and 8b. We analyze local structural relaxation, electronic properties and the relative stability of the two sites. The simulations are used to extract the contact charge density at the Eu nucleus. Then we construct the experimental isomer shift (IS) versus contact charge density calibration curve, by considering an ample set of Eu compounds: EuF3, EuO, EuF2, EuS, EuSe, EuTe, EuPd3 and the Eu metal. The, expected, linear dependence has a slope of α = 0.054 mm s−1 Å− 3, which corresponds to nuclear expansion parameter ΔR/R = 6.0 × 10−5. α allows to obtain an unbiased and accurate estimation of the IS for any Eu compound. We test this approach on two mixed-valence compounds Eu3S4 and Eu2SiN3, and use it to predict the Y2O3:Eu IS with the result +1.04 mm s−1 at the 24d site and +1.00 mm s−1 at the 8b site.

(Some figures may appear in colour only in the online journal) Yttrium sesquioxide doped with lanthanide ions are technologically relevant materials with a broad set of applications, especially owing to their luminescence properties [1]. In particular, Eu 3+ doped cubic Y 2 O 3 is one of the best red emitting phosphors [2,3]. In such doped system is of extreme importance how the rare earth ions distribute. Apart from a possible clusterization, the Eu ions can go in two symmetrically distinct sites (with 8b or 24d Wyckoff position in the Ia3 space group). A theoretical study based on atomistic * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. simulations, finds a preferential occupation of the 8b site [4]. Magnetic susceptibility experiments observe a preferential occupation of the 24d site [5,6]. On the contrary, x-ray diffraction [7] suggest a random distribution of Eu, confirmed by Mössbauer spectroscopic studies [8].
Mössbauer spectroscopy is a powerful instrument of characterization because it probes directly at the Eu site, therefore giving information on the local chemical environment. However the Mössbauer signal needs to be fitted in order to extract the valuable physical information. When dealing with multiple Eu sites with unknown occupation such fitting procedure is not straightforward and involves some biases.
In this paper we present an approach to use ab initio density functional theory (DFT) simulations in connection and as support to Mössbauer spectroscopy. The theoretical approach is introduced in section 1. In section 2 we provide an accurate theoretical characterization of Eu doping in Y 2 O 3 and the relative energetics. Simulations allow to extract the contact charge density which is the electronic density at the Eu nucleus. As a second step we construct an empirical map between the contact charge density and the isomer shift (IS) as measured in Mössbauer spectroscopy. The map (discussed in section 3) is built by performing DFT simulations on a set of Eu compounds for which the IS is experimentally known. The materials in this set are chosen to have Eu in different coordination and chemical status and such that there is one single fully occupied Eu Wyckoff position, which makes the experimental estimation of the IS more reliable. Last (section 3) we show that the map can be used to obtain the IS of other Eu compounds, like those where Eu has multiple valences. The map is finally used to predict the ISs for Y 2 O 3 : Eu.

Methods and numerical parameters
We perform density functional [9,10] theory simulations within the all-electron full-potential linearised augmentedplanewave (FP-LAPW) method [11], as implemented in the elk [12] code. Exchange correlation (xc) effects are included in the local spin density approximation (LSDA) and the LSDA + U method, where a Hubbard U correction [13] is applied to the Eu f -shell. We adopt U = 7 eV and J = 0.75 eV, as validated [14] for the isostructural Eu 2 O 3 compound; nevertheless we will discuss the impact of this choice for the results obtained in this work. All calculations include spin-orbit interactions and allow for arbitrary non-collinear spin ordering.
The Brillouin zone integration is performed in a 2 × 2 × 2k-point Monkhorst Pack mesh [15]. The maximum length of the G expansion is fixed by the R MT × |G + k| parameter set to 7. R MT is the smallest among the muffin-tin radii, which are 2.2, 2.2 and 1.6 Bohr, respectively for Eu, Y and O.
An Eu nuclear radius of 6.4 fm is used to compute the Mössbauer contact charge density, and was derived [16] from the relation R = 1.2A 1/3 fm, where A, the mass number of the Mössbauer active Europium, is taken to be 151.
Numerical results were obtained by fully minimizing internal forces, while Bravais lattice parameters are linearly interpolated, as a function of Eu content, from those of the parent compounds Y 2 O 3 and Eu 2 O 3 (10.602 and 10.859 Å respectively [17]).
All other systems studied in this work (listed in the text and used to obtain figure 3) have been studied at the experimental lattice positions and using the same computational parameters and methods as above, as well as a similarly dense reciprocal space discretization.

Structural and electronic properties of Y 2 O 3 : Eu
The host compound Y 2 O 3 crystallises in the cubic bixbyte structure with T 7 h (Ia3) symmetry, its primitive cell is body centered and composed by 10 formula units [18]. An important aspect is that there are two symmetry-inequivalent Y sites: the 8b site, has S 6 symmetry and hosts 1/4 of Y atoms; the 24d site has the lower C 2 symmetry and hosts the remaining 3/4 of the Yttrium. The structure is shown in figure 1. The local symmetry of the 8b and 24d sites is highlighted showing the respective oxygen coordination polyhedra. From the electronic point of view Y 2 O 3 is a non-magnetic insulator, with an experimental band gap of 6 eV. The Kohn-Sham DFT approximation, as expected [19], gives a smaller value which is 4.3 eV. An accurate numerical estimation can be usually achieved by many body perturbation theory in the G 0 W 0 approximation for the electronic self-energy [20,21]; however the bixbyite structure is too large and complicate for this kind of calculations. The computed density of states is shown in figure 2(a). As a validation test we have simulated the optical absorption spectrum and compared with available experimental data [22]. For a fair comparison we have included a 1.86 eV scissor correction to the conduction band. The agreement (see figure 2(b)) in the basic RPA approximation (neglecting xc effects on the response function) is fair as it reproduces the main absorption peak between 5 and 15 eV. A second structure is observed at about 30 eV which experimentally seems to be broader and merge with the main peak. The inclusion of xc contribution within the Bootstrap approximation [23][24][25] for f xc improves the shape of the main peak, and especially the form of its near-gap part (see inset in figure 2(b)).
An Y to Eu substitution can be obtained at any level up to the completely substituted system Eu 2 O 3 . However the most relevant regime is the low doping one, with ∼ 5% of Eu, since this corresponds to the maximal optical efficiency [26]. In simulations we consider a low Eu doping limit of the Y 2 O 3 crystal, and assume no interaction among Europium ions. Owing to the large unit cell size this can be simulated by the replacement of one single Y per cell. This Eu → Y substitution causes a local expansion of the O coordination polyhedra, consistently with the larger Eu ionic size. The modified bonding lengths are reported in table 1. Accurate total energy estimation shows that the 24d crystal site is more stable that the 8b by  Δ E 90 meV. This relatively high energy difference (corresponding to ∼1050 K), indicates that the 8b population should be lower than the 24d. A quantitative estimation can be done by using Boltzmann statistic and assuming that occupations to be frozen as the crystal is formed (∼1000 K); in this case then the occupation factor of the 24d sites should be of about 0.73 (against the 0.5 expected for Δ E = 0). It should be possible to observe this difference; however experiments seem to give contrasting results [5,6,8].
The value of U significantly affects the relative stability of the two sites. At U = 3 we compute that Δ E reduces to zero. But in the range of reasonable U values, which for f -electrons is between 6 and 8 eV, Δ E ranges in the interval 77-108 meV. This corresponds to a small thermal occupation error-bar of ±0.04 around the occupation population factor of 0.73, supporting the reliability of our numerical estimation. Such U dependence can be accessed only via the ab-initio calculations and is missing in the classical force fields based scheme [4], which predicted the preferable 8b population.
The magnetic quantization axis of the Eu ions is the same as for the Eu 2 O 3 , being along the (111) cubic direction for the 8b site and along the 100 for the 24d site. In simulations we obtain L, S and J to be 1.63 μ B , 2.78 μ B and 1.15 μ B for Eu in the 8b site; 2.38 μ B , 2.46 μ B , 0.07 μ B in the 24d site, with a total magnetic moment not completely quenched. This aspect in spin density functional theory (SDFT) results has been already observed and discussed in pure Eu 2 O 3 [14].
We will now consider in detail the Mössbauer properties of the Eu ion in the two crystallographic sites.

Isomer shift of Y 2 O 3 :Eu and other Europium compounds
The 151 Eu Mössbauer spectroscopy measures the impact of the electronic configuration on the nuclear transition levels. The leading effects are: IS, a change in the absorption energy [27] due to the presence of electronic charge (called contact charge density: ρ 0 ) inside the nuclear volume; the Zeeman splitting induced by the electronic magnetic moment; and the quadrupole interactions, related to the symmetry breaking effect of the crystal field. The Zeeman splitting is not present in Y 2 O 3 : Eu due to the paramagnetism of the system, while quadrupole effects are negligible for the 8b site and non-negligible for the 24d site [8]. In this work we will focus on the IS, as it carries important information on the site distribution and the bond character of the compounds. For example starting with the configuration of Eu 3+ 4 f 6 : if we add a 6s electron, the IS increases, while if we add a 4 f electron, the IS decreases due to the screening effect [28]. On the other hand, there is a very small increase if a 6p is added and a very small decrease if a 5d is added, always due to a small screening effect [28].
If one assumes an homogeneous nuclear model, the IS (with respect to the source material), is bound to the contact charge density by the relation: where ρ 0 and ρ S 0 are the contact charge density of the compound and of the radioactive source respectively and α is a proportionality coefficient. α is given by: where e is the electron charge, Z is the atomic number, R is the nuclear radius in the ground state and ΔR is the variation between the radius of the excited state (nuclear spin I = 7/2) and the radius of the ground state (nuclear spin I = 5/2). We aim at determining the coefficient α. This will be done Eu reproduced from reference [8]; black dots are the Mössbauer data of a bulk sample. Red is the total fitting curve, green and yellow are the fitted 24d and 8b components. Bottom: experimental isomer shifts [28,29] (IS) versus simulated contact charge density (Δρ 0 ) with respect to the EuF 3 compound. Red squares indicate the Eu compounds used to obtain the linear fit (pink line). Violet and blue points are validation points, which are not used for the fit as refer to compounds with more than one Wyckoff Eu site.
by means of first principle calculations of the contact charge density on several Eu compounds. Therefore, we can obtain the nuclear parameter ΔR/R not by taking a single measure of IS, but by using a fit of IS vs Δρ 0 with several measures of IS. As usual in the Mössbauer spectroscopy of 151 Eu, the value of IS is expressed with respect to a reference material which is EuF 3 . Similarly, we will refer to EuF 3 the contact charge density as well, by defining Δρ 0 := ρ 0 − ρ EuF 3 0 . This will minimize the effect of theoretical errors, computational parameters and chosen numerical approach.
The proportionality constant α between Δρ 0 and the experimental IS can then be obtained from a set of experiments and simulations, and we use the occasion of this work on Y 2 O 3 : Eu to compute it. To this end we have considered, beyond EuF 3 seven other Eu compounds, chosen in a broad valence range: EuO, EuF 2 , EuS, EuSe, EuTe, EuPd 3 and the Eu metal. For all of which the IS is known from the literature [16,28,29], and we can simulate the contact charge. The results are collected in figure 3.
The linear dependence is quite evident and allows for an accurate estimate of α = 0.054 mm s −1 Å −3 . By using this value of α and equation (2), we can calculate the nuclear parameter ΔR/R. In order to do this, it is necessary to consider that the velocity of 1 mm s −1 is actually the energy Doppler shift caused by an object which moves at 1 mm s; this converts to 1.151 × 10 −19 erg, obtained by using the energy of the transition (21.54 keV) and the expression of the relativistic longitudinal Doppler effect [30]; the final results is ΔR/R = 6.0 × 10 −5 . It is worth to point out that while the results shown in this work are based on the LDA + U, the GGA + U approach leads to very similar values for the contact charge densities and to the very same value of α = 0.054 mm s −1 Å −3 . Moreover our results do not depend critically on the chosen U parameter, as soon the value is physically reasonable.
Using the fitted curve the IS of Y 2 O 3 : Eu (or any experimentally unknown Eu compound) can be directly computed from the simulation of ρ 0 . We obtain +1.045 for the 24d site and +1.000 for the 8b site. The IS difference of 0.045 mm s −1 is small and difficult to observe in experiments without relying on an fitting procedure. The experimental data of a bulk sample, fitted with two components with the constraint of equal IS, give IS = 1.14 mm s −1 [8], in good agreement with the calculated value. On the other hand a measure on a nanocrystalline sample fitted allowing for two different ISs at the two sites, gives values of 0.97 mm s −1 for the 24d site and 1.23 mm s −1 for the 8b site [31], compatible with the calculated IS.
As a test on the IS-Δρ 0 relation and our estimation for α we consider two mixed-valence europium compounds: Eu 3 S 4 and Eu 2 SiN 3 . The Eu 3 S 4 sulfide crystallizes in the I42d space group and Eu occupies the 4a and 8d Wyckoff positions. In our SDFT simulations only 4a sites are magnetically ordered while 8d are non-magnetic, in agreement with previous studies [32]. At the two sites we obtain Δρ 0 = −194. 5

Conclusions
We have presented full-potential DFT characterization of Y 2 O 3 upon Eu doping on the two inequivalent crystallographic sites 24d and 8b. Total energy results show that the 24d sites are expected with an occupation probability of about 0.73, therefore ∼ 90% of Eu atoms should be hosted in this crystallographic environment.
We propose an approach to determine the Mössbauer IS from the knowledge of the contact charge density and the universal relation between these. Such relation is constructed by considering a broad set of Eu systems (EuF 3 , EuO, EuF 2 , EuS, EuSe, EuTe, EuPd 3 , Eu), and performing a linear fitting of the experimental IS versus the simulated contact charge for which we obtain a slope α = 0.054 mm s −1 Å −3 . The approach allows for an accurate estimation of the IS of any Eu compound in a way which is unbiased by nuclear modeling, and has been tested on the two multivalence systems Eu 3 S 4 and Eu 2 SiN 3 , for which it gives excellent results. By using the value of α, we also calculated the nuclear parameter ΔR/R = 6 × 10 −5 .
For the present case of Y 2 O 3 : Eu we derive an IS of +1.04 mm s −1 at the 24d site and +1.00 mm s −1 at the 8b site. These values are compatible with experiments, although the experimental estimation of the IS difference is difficult to resolve and depends critically on the fitting procedure.