Introduction

Ceria (CeO2) has been among the very first1 and most widely studied2 materials for catalytic and energy applications. It is used in three-way exhaust automotive catalysts2,3,4,5,6, solid-state fuel cells7,8,9,10, two-step thermochemical water splitting cycles (TWSC)11,12,13,14, low-temperature water–gas shift reactions15, and several other industrial catalytic applications16,17,18,19,20. To a large extent, the performance of ceria in these applications depends strongly on its oxygen storage capacity and facile Ce4+/Ce3+ redox reaction.

Ceria is a highly promising redox active material for TWSC14, in which a metal oxide is reduced at high temperature, T H, and subsequently re-oxidized by exposure to H2O at a lower (but still elevated) temperature T L:

$${\rm{Ce}}{{\rm{O}}_{2 - {\delta _L}}}\mathop{\longrightarrow}\limits^{{{T_{\rm{H}}}}}{\rm{Ce}}{{\rm{O}}_{2 - {\delta _{\rm{H}}}}} + \frac{{{\delta _{\rm{H}}} - {\delta _{\rm{L}}}}}{2}{{\rm{O}}_2},$$
(1)
$${\rm{Ce}}{{\rm{O}}_{2 - {\delta _{\rm{H}}}}} + \left( {{\delta _{\rm{H}}} - {\delta _L}} \right){{\rm{H}}_2}{\rm{O}}\mathop{\longrightarrow}\limits^{{{T_{\rm{L}}}}}{\rm{Ce}}{{\rm{O}}_{2 - {\delta _L}}} + \left( {{\delta _{\rm{H}}} - {\delta _{\rm{L}}}} \right){{\rm{H}}_2}$$
(2)

δ L and δ H are, respectively, the low-temperature (T L typically around 800 °C) and high-temperature (T H typically around 1600 °C) oxygen non-stoichiometries14, 21. Meredig and Wolverton22 examined the thermodynamics of these two-step reaction cycles and showed that a key thermodynamic quantity for increased efficiency is a large solid-state entropy of reduction \(\left( {\Delta S_{{\rm{red}}}^{{\rm{solid}}}} \right)\). Stoichiometric oxides typically have \({\rm{\Delta }}S_{{\rm{red}}}^{{\rm{solid}}} < 0\), so an additional source of entropy in the reduced phase is required for a high performance in TWSC. Among all the studied stoichiometric oxide reactions in ref. 22 2CeO2 → Ce2O3 + \(\frac{1}{2}\)O2 is the only one that shows \({\rm{\Delta }}S_{{\rm{red}}}^{{\rm{solid}}} >0\).

The entropy of reduction of Eq. (1) in the TWSC process is conventionally defined as23:

$${\rm{\Delta }}{S_{{\rm{red}}}} = {\rm{\Delta }}{S_{{\rm{conf}}}} + {\rm{\Delta }}{S_{{\rm{vib}}}} + \frac{1}{2}{S_{{{\rm{O}}_2}}}$$
(3)

where ΔS conf includes both electronic and ionic configurational entropy of ceria. The electronic part associated with electron or hole localization contributes to the total entropy in the same fashion as the ionic configurational entropy due to vacancies. Both configurational contributions can be determined accurately by Monte Carlo simulations based on a cluster expansion24. ΔS vib is the vibrational entropy of reduction and \(\Delta {S_{{{\rm{O}}_2}}}\) is the gas phase entropy. The latter does not depend on the choice of materials, and at temperatures higher than 1000 K, which is relevant for TWSC, is ~ 15 k B per oxygen atom23, 25—hereafter, all the reported entropic contributions are per one oxygen vacancy, unless stated otherwise. By excluding the gas phase entropy, we have the solid-state entropy of reduction \(\left( {{\rm{\Delta }}S_{{\rm{red}}}^{{\rm{solid}}}} \right)\), which is the material dependent contribution.

The first experimental studies on the entropy of reduction of CeO2 by Bevan et al.26 and successively Panlener et al.23 showed that this quantity has a logarithmic dependence on the non-stoichiometry, δ, of CeO2−δ . They associated this behavior with the change in the configurational entropy (ΔS conf) with varying degrees of non-stoichiometry, while vibrational entropy was estimated to be negligible. Grieshammer et al.27 calculated the ΔS vib of defect formation in ceria at a fixed non-stoichiometry and found it to be ~ 2.5 k B, which is not negligible, but smaller than other entropic contributions. Gopal et al.28 did comprehensive calculations using Monte Carlo simulations based on the DFT-derived cluster expansion Hamiltonian. They calculated all configurational and vibrational entropic contributions (see Eq. (3)) for different values of δ at a temperature of 1480 K and found that the actual configurational entropy is much smaller than that of commonly assumed ideal solution model. For δ > 0.12 their calculated total entropy agrees with that of experiment, however, in the δ range of 0.01–0.12 they underestimate the experimental entropy, leading to a 4.5 k B 28 gap between their calculations and the experimental measurement23.

In the present paper, we demonstrate that particularly in the case of lanthanides, a different type of electronic entropy should be considered. This source of electronic entropy, which becomes important for elements with partially filled f shells arises when electrons can be distributed over a large number of multiplet states. This entropy results only from onsite configurational entropy in extremely localized f orbitals, which stems from the possible configurations associated with occupations of the same atomic orbitals, and is hereafter denoted \({\rm{\Delta }}S_{{\rm{elec}}}^{{\rm{onsite}}}\). Our results, show that the \({\rm{\Delta }}S_{{\rm{elec}}}^{{\rm{onsite}}}\) contributes to the high performance of ceria in TWSC and fills the ≈4.7 k B gap between previous theoretical and experimental entropy values. We calculate the entropy of reduction (M n+/M (n−1)+) for several other lanthanide cations (praseodymium, neodymium, europium, and terbium) that are stable in two valence states and for which reliable spectroscopic data are available (see ref. 29). We calculate the crystal field (CF) parameters of Ce3+ in the host fluorite CeO2 structure, where each cerium atom experiences a cubic crystal field from eight oxygen atoms. The results show that the electronic entropy of reduction of Ce4+ is the highest among f elements, even surpassing the calculated configurational entropy28 at large off-stoichiometry values, thereby explaining the unique entropic properties of the reduction of CeO2.

Results

L − S coupling and crystal-field splitting

Electronic configurational entropy \(\left( {S_{{\rm{elec}}}^{{\rm{onsite}}}} \right)\) arises from thermal excitations among orbital microstates created by Russel–Saunders (L–S) coupling. In multi-electron atoms, L–S coupling between orbital and spin angular momenta (L–S) results in the creation of microstates, which are further split by CF interactions. As the interaction of the very contracted f-orbitals of the lanthanide ions with the CF is small30, it is reasonable to first treat the electronic energy levels of the lanthanide oxides by considering only free ion energies and then subsequently apply the crystal field. For f-orbitals of lanthanides, defining the electronic configuration by just number of valence electrons (4f n) is far less descriptive than the term symbol described by the L–S coupling scheme31. In this scheme, coupling of orbital and spin angular momentum results in 2S+1 L J term symbols in which 2S + 1 is the spin-multiplicity, L is the total orbital quantum number and J is the total angular momentum quantum number, ranging from |L + S| to |L − S| by steps of one. Indeed, Hund’s rules imply that the term with the largest value of multiplicity (2S + 1) is the most stable one. If several terms have the same multiplicity, then the term with the largest L is the most stable one. If the f shell is less than half-filled then the state with the lowest J has the lowest energy. On the other hand, if the f shell is more than half-filled then the state with the highest J has the lowest energy. The degeneracy of each J-multiplet is (2J + 1) and the total number of microstates (m) for a given term symbol 2S+1 L is (2S + 1) × (2L + 1).

When a free ion is placed in a crystal, the CF further splits each of the degenerate J states to several subsets and breaks the spherical symmetry of the f-shell charge distribution, depending on the local symmetry of the ionic environment. Here, we used a fully ab-initio method, opposing crystal potential (OCP)32, to calculate the crystal field parameters of Ce3+ in the host fluorite CeO2 structure, where each cation is coordinated to eight oxygen atoms. This method is motivated by the above mentioned lowering of the spherical symmetry of the d/f charge distribution, e.g., to cubic symmetry in CeO2, due to crystal field interactions. OCP iteratively computes and subtracts CF interactions in order to effectively cancel them and recover spherical d/f charge distribution upon convergence. An on-site potential \({\lambda _{mm'}}\) is introduced as a matrix of Lagrange multipliers in constrained DFT calculations, and crystal field parameters are then obtained as linear combinations of \({\lambda _{mm'}}\) (see ref. 32 for details). As previously discussed32, the goal of the constrained DFT calculations was to extract CF parameters rather than to introduce band gap and self-interaction corrections. Consequently the commonly adopted LDA + U approach was not necessary. As a benchmark, the CF parameters of PrO2 from constrained OCP calculations were in good agreement with our previous LDA + U calculations that directly compute the energy difference of different crystal field levels33.

Figure 1 shows the f 1 (Ce3+) energy-level splitting scheme in the presence of spin-orbit coupling (SOC) and calculated CF. Without CF, the f 1 states split into 2 F 5/2 and 2 F 7/2 separated by approximately 0.28 eV29. The CF interaction further splits the six-fold degenerate 2 F 5/2 ground state into a four-fold degenerate Γ8 and two-fold degenerate Γ7 subsets, separated by 0.12 eV. Crystal field—calculated by OCP method for Ce3+—splits the eight-fold degenerate 2 F 7/2 state into states with energies 0.25, 0.32, 0.46 eV. These results are comparable with experimental and theoretical CF splittings of Ce3+ doped in CaF2 34,35,36 and YAG (Y3Al5O12)34, 37 hosts.

Fig. 1
figure 1

Energy levels of the 4f 1 orbital of Ce3+. Ce3+ splits initially by spin-orbit coupling (SOC) and subsequently by cubic crystal field (CF) of the the fluorite structure. The spin-orbit splitting between J = 5/2 and J = 7/2 is about 0.28 eV29, 49. The color gradient (see color bar) indicates the probability distribution at 1900 K, given by exp(−E i /kB T), and numbers in parentheses stand for the degeneracy of the electronic states. The first predicted Γ8 → Γ7 excitation for CeO2 is 0.12 eV. Predictions for the higher CF levels of J = 7/2 are 0.25, 0.32, 0.46 respectively

Full size image

Onsite electronic entropy

By having the energy E i and degeneracy (g i ) of each microstate, \(S_{{\rm{elec}}}^{{\rm{onsite}}}\) of a system with m different microstates can be calculated by

$$S_{{\rm{elec}}}^{{\rm{onsite}}} = - {k_{\rm{B}}}\mathop {\sum}\limits_i^m {{g_i}\,{p_i}\,{\rm{ln}}\,{p_i}} $$
(4)

where the probability of thermal excitation (p i ) to the state with energy E i is proportional to the Boltzmann factor by

$${p_i} = \frac{{{\rm{exp}}\left( { - {E_i}{\rm{/}}{k_{\rm{B}}}T} \right)}}{Z}$$
(5)

and Z is the partition function defined by

$$Z = \mathop {\sum}\limits_i^m {{g_i}\,{\rm{exp}}\left( { - {E_i}{\rm{/}}{k_{\rm{B}}}T} \right)} $$
(6)

Equations (4) and (6) indicate that \(S_{{\rm{elec}}}^{{\rm{onsite}}}\) directly depends on two key factors; the probability (p) and the total number of microstates (m) that varies from 14 in f 1/f 13 to 66 in f 5/f 9 (see Table 1 for the number of degenerate microstates in other occupations). From Ce to Nd, both the total number of microstates and SOC increase. A stronger SOC implies a larger multiplet splitting between the levels which leads to microstates with higher energies. At low temperature, i.e., limited thermal excitations, those high-energy microstates are less probable to be occupied. However, in Ce3+ with a moderate SOC, at temperatures relevant for TWSC (T ≈ 1900 K), a large fraction of microstates are accessible due to the increased thermal excitation, as seen in Fig. 1. Therefore, \(S_{{\rm{elec}}}^{{\rm{onsite}}}\) of Ce3+ approaches the ideal limit of k B ln(m). As seen in Table 1, \(S_{{\rm{elec}}}^{{\rm{onsite}}}\) weakly depends on the number of electrons and degenerate states. Once the f orbitals get occupied by a single electron—generating 14 microstates with ln (14) k B (≈2.6 k B) entropic contribution—the system gains a large entropy. For example \(S_{{\rm{elec}}}^{{\rm{onsite}}}\) of f 1 (Ce3+) is 2.34 k B and f 4 (Nd2+) is 3.52 k B.

Table 1 Calculated \({\bf{S}}_{{\bf{elec}}}^{{\bf{onsite}}}\) per oxygen vacancy of selected lanthanide ions before and after reduction at 1900 K
Full size table

Myers et al.38 extracted the electronic entropy contribution of lanthanide ions (Ln3+) in lanthanide trihalides from absolute entropy data. Our calculated electronic entropies per ion at ≈300 K in units of k B compared with Myers et al.38 data (value inside parentheses) are the following: Ce3+, 1.79 (1.77); Pr3+, 2.19 (2.18), Nd3+, 2.30 (2.27); Eu3+, 1.13 (1.10); Tb3+, 2.56 (2.54). The calculated \(S_{{\rm{elec}}}^{{\rm{onsite}}}\) based on L–S coupling shows excellent agreement with previously reported data. Below, we will consider the effect of crystal field splitting on the calculation of \(S_{{\rm{elec}}}^{{\rm{onsite}}}\).

For TWSC applications, the absolute electronic entropy does not matter, only the entropy difference before (f n) and after (f n+1) reduction is relevant, \(\Delta S_{{\rm{elec}}}^{{\rm{onsite}}} = 2\left( {S_{{\rm{elec}}}^{{\rm{n}} - 1} - S_{{\rm{elec}}}^{\rm{n}}} \right)\). The factor two is due to the fact that two Ce4+ ions are reduced per oxygen vacancy. As discussed in the previous paragraph and shown in Table 1, because of the large \(S_{{\rm{elec}}}^{{\rm{onsite}}}\) gain upon occupation of an f orbital, in TWSC systems the electronic entropy of reduction reaches its maximum, when the entropy of the oxidized state is zero, \(S_{{\rm{elec}}}^{{n}} = 0\). As a result, the largest \({\rm{\Delta }}S_{{\rm{elec}}}^{{\rm{onsite}}}\) is found in Ce4+ → Ce3+, which undergoes an f 0 → f 1 redox reaction. Having the oxidized state f 0 (1 S) with zero onsite electronic entropy is a unique feature of ceria, resulting in a large \({\rm{\Delta }}S_{{\rm{elec}}}^{{\rm{onsite}}}\) of 4.68 k B per oxygen vacancy, which is a maximum for the reduction of any rare-earth cation. We assert that this unique entropic characteristic of the Ce4+/Ce3+ redox reaction helps facilitate the TWSC properties of CeO2. The second largest value of \({\rm{\Delta }}S_{{\rm{elec}}}^{{\rm{onsite}}}\) is found in terbium (Tb4+ → Tb3+) with 2.30 k B per oxygen vacancy at 1900 K (Fig. 2). There, the non-reduced Tb4+ has a half-filled shell with only eight spin-degenerate (2S + 1) states but with only one orbital degeneracy (L = 0), while Tb3+, has an orbital degeneracy of 7 (L = 3), providing extra entropy, see Table 1. This extra source of entropy could make Tb4+ based materials promising candidates for TWSC applications, as Tb, like Ce, is stable in two valence states (Tb4+/Tb3+). This prediction agrees with a recent thermodynamic study that also suggested39 TbO2 as a potential candidate for TWSC applications.

Fig. 2
figure 2

Calculated \({\rm{\Delta}}S_{{\rm{elec}}}^{{\rm{onsite}}}\) for lanthanides ions. Predicted electronic entropy of reduction per oxygen vacancy for the lanthanide oxides studied in this work. At high temperature, reduction of CeO2 has the highest \(\Delta S_{{\rm{elec}}}^{{\rm{onsite}}}\) followed by reduction of TbO2 (see Table 1)

Full size image

In contrast to TbO2, reduced Eu2+ has a half-filled shell with eight spin- and only one orbital-degeneracy, resulting in a temperature independent entropy of 2 × ln(8)k B = 4.16 k B per oxygen vacancy. On the other hand non-reduced Eu3+ has a large number of total degeneracy of 49 (see Table 1). Since \({\rm{\Delta }}S_{{\rm{elec}}}^{{\rm{onsite}}} = 2\left( {S_{{\rm{elec}}}^{{\rm{E}}{{\rm{u}}^{2 + }}} - S_{{\rm{elec}}}^{{\rm{E}}{{\rm{u}}^{3 + }}}} \right)\), by increasing temperature \(S_{{\rm{elec}}}^{{\rm{E}}{{\rm{u}}^{3 + }}}\) rapidly increases which decreases \({\rm{\Delta }}S_{{\rm{elec}}}^{{\rm{onsite}}}\) making it negative at high temperature.

To examine the effect of crystal field on the reduction entropy, we calculated CF parameters for Ce3+ in the CeO2 host lattice. We should note that a previous study40 using inelastic neutron scattering demonstrated that reducing nearest neighbor oxygen atoms have little effect on the overall energy level. This indicates that formation of vacancy next to Ce3+ can not dramatically change the CF splitting of a perfect cubic field. Using ab initio CF parameters, previous calculations41 of \(S_{{\rm{elec}}}^{{\rm{onsite}}}\) for actinides showed an excellent agreement with the available experimental data. Here, we calculate electronic entropy of Ce3+. As seen in Fig. 2, CF only affects the \({\rm{\Delta }}S_{{\rm{elec}}}^{{\rm{onsite}}}\) at low temperatures while at temperatures higher than 1000 K the CF effect is small; therefore, one can just rely on J-multiplet states resulting from L–S coupling. Table 1 shows that at such elevated temperatures, even pure \(S_{{\rm{elec}}}^{{\rm{onsite}}}\) weakly depends on CF. For example, the effect of CF on \(S_{{\rm{elec}}}^{{\rm{onsite}}}\) of Pr4+ and Ce3+ is about < 3%. This percentage is in fact an upper limit, as the calculated CF is at DFT lattice parameters at zero Kelvin. The temperature dependent measurement of CF splitting by Walsh et al.42 showed that the CF splitting significantly decrease with temperature and lattice thermal expansion. Therefore, at 1900 K, which corresponds to typical temperature required for TWSC applications, CF field plays only a minor role.

Other entropy contributions

Finally, we compare \({\rm{\Delta }}S_{{\rm{elec}}}^{{\rm{onsite}}}\) with the other sources of entropy. For simplicity we consider a fixed composition of δ = 0.03 roughly corresponding to one oxygen vacancy in a 96-atom supercell. For this composition, we were able to find several reported experimental and theoretical data points (Table 2). At this composition the calculated ΔS vib is ~ 2.5 k B 27. The ΔS conf of CeO2−δ , assuming ideal mixing entropy23, 28 (ΔS c = −nk B ln(δ)), where n depends on the defect structure and here n = 3, is ≈10.4 k B. However, we note that a system with extensive ordering of oxygen vacancies43, such as ceria, will have short range order and hence the actual configurational entropy is non-ideal and smaller than the ideal solution model. For instance, the non-ideal ΔS conf + ΔS vib, calculated by Monte Carlo simulation based on a cluster expansion Hamiltonian, is about 5.9 k B 28, less than half of the ideal ΔS conf. Our calculations show that the neglected electronic entropy (ΔS elec) is >4.7 k B, which is comparable to these other widely considered sources of entropy and can explain the ≈5 k B gap between the calculation and experiment. We note that as long as oxygen vacancy is compensated by two polarons Ce′ (i.e., Ce4+) \(\left( {\left[ {{\rm{Ce}}\prime\!_{{\rm{Ce}}} } \right] = 2\left[ {V_{\rm{O}}^{..}} \right]} \right)\), \({\rm{\Delta }}S_{{\rm{elec}}}^{{\rm{onsite}}}\) is not a function of off-stoichiometry (δ). Being independent of off-stoichiometry implies that at large δ the contribution from the electronic entropy surpasses that of the configurational entropy (which decreases with δ) and becomes the major entropy contribution. Using the calculated ΔS conf in ref. 28, we estimate that this crossover occurs at (δ ≈ 0.05).

Table 2 Contribution of different entropic terms for δ = 0.03 and temperature of 1500 K
Full size table

Our results show that the electronic contribution to the entropy of reduction explains the gap between the results of the currently most detailed theoretical calculations of ref. 28 and the experimental data of Panlener et al.23 for small deviations from stoichiometry corresponding to δ < 0.03. At larger deviations, adding a constant onsite electronic entropy to the vibrational and configurational entropies from ref. 28 would overestimate the experimental data. There could be several reasons for this apparent discrepancy. For instance, at higher δ values most of the polarons become bound to oxygen vacancies forming singly charged V\(_{\rm{O}}^{2 - }\)–Ce3+ or neutral V\(_{\rm{O}}^{2 - }\)–2Ce3+ complexes44; the proximity of Ce3+ to an oxygen vacancy could slightly modify the electronic structure and hence reduce the electronic entropy associated with Ce3+, but as already discussed the overall effect of oxygen vacancy on the energy levels40 and electronic entropy is expected to be small. Furthermore, the experimental measurements of Panlener et al.23 found that the enthalpy of reduction is composition dependent even at very small δ; however, this finding has been challenged due to the large experimental uncertainty23, 26. As the entropy is obtained from TΔS = ΔH − ΔG, the entropy values of Panlener et al. may be contaminated by contributions from the composition dependent contribution to ΔH. Indeed, the results of ref. 28 suggest that the entropy stays approximately constant for δ in the 0.05 to 0.15 range, while the data of Panlener et al.23 show a pronounced decrease in this range.

Measurements of the Seebeck coefficient provide another means of estimating the electronic entropy contribution in the dilute limit where all polarons are unbound6, 45. Unfortunately, the experimental data here are also contradictory. The data of Tuller and Nowick45 suggest that for small δ the spin degeneracy factor is one, which contradicts the Kramers theorem requiring that the ground state must be at least doubly degenerate. However, a later study by the same authors6 concluded that the agreement between the polaron model with spin degeneracy one and the experimental data for the Seebeck coefficient was poor, especially at low δ where impurities were thought to play an important role. On the theory side, the vibrational entropy of an isolated Ce3+ polaron has not been established accurately. Grieshammer et al.27 have calculated a very large value of about 7 k B for the entropy of polaron formation at zero pressure, but the largest contribution to this value is due to a volume contribution from the CeO2 host, which was treated in an approximate fashion. Such a large positive entropy is inconsistent with the available data on the Seebeck coefficients in the dilute limit6, 45. Hence, thermoelectric measurements on pure, well-equilibrated samples of CeO2 and more accurate calculations of the vibrational entropy associated with free polarons are highly desirable.

Discussion

We calculated electronic entropies of different lanthanides in the presence of SOC and CF. We calculated CF splittings for Ce3+ and Pr4+ and found that at temperatures above 1000 K, CF interactions affect the S elec by >3%. The results show that, in ceria, the magnitude of the entropy of reduction due to the commonly neglected onsite electronic entropy (ΔS elec) reaches a maximum of 4.7 k B per oxygen vacancy, which is twice as large as the vibrational entropy contribution and can be larger than the configurational entropy. This surprisingly large entropy is the result of the very unique electronic structure of cerium in ceria, where redox reactions change its electronic state from f 0 to f 1. These entropic properties, together with the excellent chemical stability and tolerance for large non-stoichiometry, put ceria in a unique position for two-step solar thermochemical CO2/H2O splitting cycles. In addition, we find that Tb (IV) based materials have the next highest electronic entropy, for Tb4+ → Tb3+ redox reactions. We therefore propose compounds containing Tb4+ should be experimentally investigated as promising candidates for TWSC applications.

Methods

Constrained DFT calculations, required for the calculation of crystal field parameters using OCP method, were performed in the VASP package46 using the Perdew–Becke–Ernzerhof47 functional, projector augmented-wave potentials48, a 6 × 6 × 6 k-point mesh an energy cut-off of 520 eV and a convergence criterion of 10−10 eV/atom for the electronic structure to ensure satisfactory convergence of CF parameters. One electron was added to the neutral CeO2 cell in order to model Ce3+. Note that as previously discussed32, the goal of the constrained DFT calculations was only to extract CF parameters rather than to introduce band gap and self-interaction corrections.

Data availability

The data that support the findings of this study are available within the paper or from the corresponding author on request.