Theoretical study of mixing energetics in homovalent fluorite-structured oxide solid solutions

Journal of Nuclear Materials 444 (2014) 292–297

Contents lists available at ScienceDirect

Journal of Nuclear Materials journal homepage: www.elsevier.com/locate/jnucmat

Theoretical study of mixing energetics in homovalent fluorite-structured oxide solid solutions Vitaly Alexandrov a,⇑, Niels Grønbech-Jensen b,c, Alexandra Navrotsky e,b, Mark Asta d,e,b a

Physical Sciences Division, Pacific Northwest National Laboratory, Richland, Washington 99352, USA Department of Chemical Engineering and Materials Science, University of California, Davis, California 95616, USA c Department of Mechanical and Aerospace Engineering, University of California, Davis, California 95616, USA d Department of Materials Science and Engineering, University of California, Berkeley, California 94720, USA e Peter A. Rock Thermochemistry Laboratory and NEAT ORU, University of California, Davis, California 95616, USA b

a r t i c l e

i n f o

Article history: Received 3 July 2013 Accepted 1 October 2013 Available online 12 October 2013

a b s t r a c t Mixing energies (DHmix) for fluorite-structured (Zr1xCex)O2 and (Th1xCex)O2 solid solutions are computed from density functional theory (DFT), employing cluster-expansion (CE), special-quasirandomstructure (SQS), and continuum-elasticity approaches. These systems are of interest as models for actinide-dioxide mixtures, due to the availability of calorimetric data which allows a direct assessment of the accuracy of the different computational methods for calculating DHmix in such fluorite-structured solid solutions. The DFT-based SQS and CE results for solid solutions with random configurational disorder are in very good agreement, and are used along with the calorimetry data to test the accuracy of a linear-elasticity model which allows predictions of the DHmix under the assumption that the dominant contribution in these homovalent solid solutions arises from elastic strain energy. The linear-elasticity models describe the mixing energies to within an accuracy of approximately 2 and 0.1 kJ/mol for the Zr and Th based systems, respectively. The excellent accuracy for the ThO2-based system is interpreted to result from the smaller size mismatch, and corresponding high accuracy of the linear elasticity approximation. We thus apply elasticity theory to estimate the magnitudes of DHmix for (Th1xMx)O2 and (U1xMx)O2 actinide-dioxide solid solutions, with M = U, Th, Ce, Np, Pu and Am, for which the degree of size mismatch is comparable to that in (Th1xCex)O2; the results yield elastic contributions to DHmix with a maximum magnitude of 3 kJ/mol. Ó 2013 Elsevier B.V. All rights reserved.

1. Introduction The actinide dioxide UO2 and ThO2 systems form the basis for the majority of fuels used in nuclear energy production worldwide. During burnup in a nuclear reactor, a variety of fission and transmutation products form within the fuel which are soluble in the fluorite structure of these compounds [1,2]. As a consequence, spent nuclear fuels generally feature a matrix phase that is a complex multicomponent solid solution of UO2 or ThO2, with a significant concentration of different types of minor actinide and lanthanide solute elements. Although the energetics associated with the formation of these solid solutions is an important factor governing thermochemical properties, experimental characterization of the mixing properties for the actinide dioxide systems relevant to nuclear fuels remains incomplete. The practical challenges associated with direct measurements of such properties, particularly for systems containing transuranic elements, have led to increasing interest in the application of first-principles elec⇑ Corresponding author. Tel.: +1 5093716462. E-mail address: [email protected] (V. Alexandrov). 0022-3115/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jnucmat.2013.10.001

tronic-structure methods as a framework for computing the energetic and thermodynamic properties of actinide-dioxides and their solid solutions (e.g., [1,3–6]). In the present study we consider mixing energetics in the simplest class of fluorite-structured actinide-dioxide solid solutions, namely (Th1xMx)O2 and (U1xMx)O2 systems containing homovalent cation species (Th, U and M, with 4+ oxidation states) and stoichiometric oxygen concentrations. We consider first (Zr1xCex)O2 and (Th1xCex)O2 as relevant model systems for which previously published calorimetric measurements of the mixing enthalpies (DHmix) [7,4] can be used to assess the accuracy of the current and previously published [5,4,6] first-principles computational results. We note that in addition to serving as useful model systems, (Zr1xCex)O2 and (Th1xCex)O2 are also directly relevant in the context of nuclear fuels research, as both ThO2 and ZrO2 have been previously considered as candidate inert matrix fuels (IMF), and ZrO2 is an oxidation product of zirconium alloys used as fuel cladding material [8,9,2]. Recently it was suggested that the mixing energetics of (Th1xCex)O2 solid solutions are dominated by an elastic energy contribution originating from the mismatch in cation size (ionic

V. Alexandrov et al. / Journal of Nuclear Materials 444 (2014) 292–297

radius) [4]. Further, it was shown that accurate results for DHmix could be obtained from continuum linear-elasticity theories of the energetics of solid solutions [10], based solely on a knowledge of the lattice parameters and single crystal elastic constants of the constituent oxides. Here we investigate the applicability of this model in further detail, by comparing DFT and experimental values of DHmix for both (Zr1xCex)O2 and (Th1xCex)O2 with results obtained from the linear elasticity formulas. The comparisons show that the elasticity theory approach is more accurate for (Th1xCex)O2, for which the cation size mismatch is Dr=r ¼ 8% (where Dr and r denote the difference and average, respectively, of the cation ionic radii), than it is for (Zr1xCex)O2, for which Dr=r is much larger at 14%. These results likely reflect a larger magnitude of anharmonic contributions to the elastic energy for the (Zr1xCex)O2 system, which displays a larger size mismatch and an intrinsic instability of the fluorite structure for one of the constituent oxides (ZrO2). The elasticity model is used to estimate the magnitudes of DHmix for (Th1xMx)O2 and (U1xMx)O2 actinide-dioxide solid solutions, with M = U, Th, Ce, Np, Pu and Am, for which the constituent oxides are mechanically stable in the fluorite structure and the ionic size mismatch is comparable to that in (Th1xCex)O2, such that the elasticity-theory predictions are expected to be reasonably accurate. The remainder of the paper is organized as follows. In the next section the details of the computational approaches employed in the current study are reviewed, along with the formulas describing the theory for the elastic energy of cubic solid solutions. Results for (Zr1xCex)O2 and (Th1xCex)O2 are presented in Section 3, along with the continuum elasticity predictions for (Th1xMx)O2 and (U1xMx)O2 actinide-dioxide solid solutions. The results are summarized in Section 4.

2. Methodology 2.1. Ab initio total energy calculations Ab initio total energy calculations were performed within the DFT plane-wave formalism using the generalized gradient approximation of Perdew, Burke and Ernzerhof (GGA-PBE) for the exchange–correlation functional [11]. For the ThxCe1xO2 and ZrxCe1xO2 systems we also examined the effect of exchange–correlation potential by performing additional mixing-energy calculations for SQS structures at x = 0.5, using the GGA+U formalism of Dudarev et al. [12] with U = 4 eV for the Ce and Th f orbitals. We obtained results that agreed to within 1 kJ/mol with those obtained with standard GGA functional. The close agreement between the GGA and GGA+U results are expected due to the fact that in the 4+ oxidation state Th and Ce cations have 5f0 and 4f0 electronic configurations, respectively. All calculations were performed within the projector augmented wave (PAW) formalism [13] as implemented in the Vienna Ab initio Simulation Package (VASP) [14–16]. In the calculations we employed a plane-wave cutoff energy of 500 eV and Monkhorst– Pack [17] grids for Brillouin zone integration corresponding to a 8  8  8 k-point mesh for primitive cells of fluorite MO2. Structural relaxations were performed using a conjugate-gradient algorithm until atomic forces were converged to magnitudes within 1 meV/Å. The structural relaxations were followed by an additional static calculation of the energy for each optimized geometry using tetrahedron integration with Blöchl corrections [18]. To calculate single-crystal elastic constants, for use in the elasticity models, we started by determining the equilibrium lattice constant a, bulk modulus B, and pressure derivative of the bulk modulus B0 for each fluorite-structured compound from a least-squares fit of the ab initio energy-volume data using the

293

third-order Birch equation of state. In the calculations of energyvolume data the total energy was converged numerically to less than 108 eV. The bulk modulus is related to elastic constants for cubic lattices as B = (C11 + 2C12)/3. To derive independent values for each of the cubic single-crystal elastic constants (C11, C12 and C44), we follow the approach described in Ref. [19] in which (C11–C12) and C44 are computed by applying volume-conserving orthorhombic and monoclinic strains, respectively, and fitting the total energy as an even polynomial function of the imposed strains. 2.2. Cluster expansion and special quasirandom structures To compute mixing energetics for the (AxB1x)O2 fluorite-structured solid solutions, we employ the DFT-based cluster expansion (CE) formalism, as implemented in the alloy theoretic automated toolkit (ATAT) [20,21]. According to the CE formalism [22,23], the total energy for an arbitrary arrangement of cations on the fcc sublattice of the fluorite crystal structure can be parametrized as

EðrÞ ¼

X ma J a h/a ð~ rÞi

ð1Þ

a

where ma is the multiplicity of cluster a, Ja is the effective cluster interaction coefficient, and h/a ð~ rÞi is the lattice average of the cluster functions taken over all symmetrically-equivalent clusters of type a. In ATAT the coefficients Ja are obtained from a least-squares fit to a set of DFT-calculated energies for hypothetical ordered cation arrangements, and the optimal set of clusters is chosen to maximize the predictive power of the fit, as estimated through the cross-validation score (for further details see [21]). In this work, the CE coefficients were fit to a database of 137 DFT energies for structures obtained by generating all possible symmetry-distinct atomic arrangements in binary solid solutions containing up to 18 atoms total. The CV scores for the CE fits were 0.027 kJ/mol-cation (0.28 meV/cation) for (Th1xCex)O2, based on a cluster expansion that contained 29 pair and 33 triplet clusters. For (Zr1xCex)O2 a cluster expansion containing 29 pair, 33 triplet and 1 quadruplet clusters produced a CV score of 0.22 kJ/mol-cation (2.3 meV/cation). The special-quasirandom structure (SQS) approach [24] is an alternative method for calculating the energetics of randomly disordered solid solutions. For modeling such systems, the best SQS for a given composition corresponds to a structure that matches the correlation functions for a random alloy for a given target set of clusters as closely as possible. In this study, we use 96-atom SQSs (containing 32 cations and 64 oxygen anions) [25] to obtain values of DHmix for (AxB1x) O2 solid solutions sampling x from 0 to 1 with a step of 0.125, for comparison with the predictions of the CE approach. In addition to providing an independent check of the CE results, the SQS structures are also used in the current work to compute the magnitudes of the ionic displacements induced by the cation size mismatch, as described in Section 3.1 below. 2.3. Elasticity theory Microscopic elasticity theory [10] provides a framework for computing the elastic energy contribution to the mixing energetics of a homogeneous solid solution. The approach relies on the knowledge of the force constants in the solid and the Kanzaki forces characterizing the strains induced by insertion of solute atoms into the matrix. In the limit of a random solution, a longwavelength approximation for the force constants and Kanzaki forces can be used to compute an estimate of the elastic contribution to the mixing enthalpy (Eelast) based solely on the knowledge of the lattice constants and elastic constants. For the case of a cubic

294

V. Alexandrov et al. / Journal of Nuclear Materials 444 (2014) 292–297

binary substitutional solid solution, the resulting expression takes the form [10]:

Eelast ¼

  1 C 11 þ 2C 12 ðC 11 þ 2C 12 Þu2 3  huðnÞin xð1  xÞ 2 C 11

ð2Þ

where Cij denote the single-crystal elastic constants of the solid solution, u = dln(a)/dx is the concentration coefficient of the crystal lattice expansion produced by substitutional atoms, and n = k/k denotes the unit vector defining the direction of a concentration wave in k-space. The average hu(n)in can be approximated as [10]:

huðnÞin  1  

4 C 11  C 12  2C 44 5 C 11 þ C 12 þ 2C 44

18 ðC 11  C 12  2C 44 Þ2 35 ðC 11 þ 2C 12 þ 4C 44 ÞðC 11 þ C 12 þ 2C 44 Þ

ð3Þ

Since the atomic volume of the solid solutions was found to be well approximated by a linear function of composition (Vegard’s law), the parameter u in the above equations can be simply estimated from the difference in the lattice constants of the constituent fluorite-structured oxide compounds. For a given composition (x) of the solid solution, we use concentration-weighted values of elastic constants of constituent oxides as the values of Cij in Eqs. (2) and (3). For the systems with nonzero f orbital occupancy (UO2, NpO2, PuO2, AmO2), we used lattice and elastic parameters calculated previously using the DFT+U approach (see Table 2). 3. Results and discussion 3.1. First-principles results for CeO2-based binary solutions Fig. 1 plots CE and SQS results for the mixing energies of (Th1xCex)O2 and (Zr1xCex)O2 solid solutions in the left and right panels, respectively. In these plots, the mixing energies are defined as:

DHmix ¼ E½ðM1x Cex ÞO2   ð1  xÞEðMO2 Þ  xEðCeO2 Þ

ð4Þ

3.5

Formation enthalpy, kJ/mol-cation

Formation enthalpy, kJ/mol-cation

where E[(M1xCex)O2] denotes the energy of a relaxed fluoritestructured ordered compound or disordered solid solution, while E(MO2) and E(CeO2) denote the energies of the constituent oxides in the cubic fluorite crystal structure, optimized with respect to lattice parameter. Note that this definition of the mixing energy differs from the mixing enthalpy measured experimentally by a PDV contribution, where P denotes pressure and DV is the mixing volume. For the values of DV characteristic of the systems studied in this work, the PDV contribution is negligible for a pressure of P = 1 atm. The cross symbols in Fig. 1 denote values of DHmix computed for the hypothetical ordered compounds used to parameterize the CE for each system, while the solid lines represent the CE predictions for solid solutions with random configurational disorder. Also

CE: ordered structures CE: random solid solutions SQS: random solid solutions

3 2.5 2 1.5 1 0.5 0

0 ThO2

0.2

0.4

0.6

Ce mol fraction

0.8

plotted in Fig. 1, with open diamond symbols, are the SQS predictions for the mixing energies of random solid solutions. All of the calculated mixing energies are positive and nearly symmetric about x = 0.5, suggesting the presence of a miscibility-gap phase diagram for fluorite-structured solid solutions, consistent with previous calculations [26,5,4,6]. Very good agreement is obtained between the CE and SQS results, with the maximum discrepancy between the two approaches found to be on the order of 2 kJ/mol for (Zr1xCex)O2 with x = 1/3. The magnitude of the mixing energy for the (Zr1xCex)O2 system is roughly five times larger in magnitude than that in (Th1xCex)O2. Table 1 compares the calculated mixing energies with experimental calorimetry results [7,4]. The experimental data for DHmix in both systems is well described by a regular solution model, with the functional form DHmix = Xx(1  x), which is symmetric about x = 0.5, in good agreement with the present calculations. Consequently, we compare experimental measurements and calculated results in Table 1 only for x = 0.5. Note that for (Zr1xCex)O2 the values of DHmix are taken with respect to the energy of ZrO2 in the cubic fluorite structure, despite the fact that this compound is stable in a monoclinic structure at low temperatures. Thus, the calorimetry measurements have been corrected by the enthalpy difference between the monoclinic and cubic structures of ZrO2, equal to 8.8 kJ/mol [7]. In comparing experimental and theoretical results listed in Table 1, it is important to keep in mind that the measured data refer to room temperature, with the calculations correspond to zero temperature. Thus, the results differ by a vibrational contribution to DHmix, which is however expected to be negligibly small at room temperature. Another potential difference between the calculations and measurements relates to the degree of configurational disorder present in the solid solutions. The calculations have assumed random configurational disorder, and thus neglect any contributions from (clustering-type) short-range order (SRO) which may be present in the experimental solid solutions that are equilibrated at high (but finite) temperature. Using Monte-Carlo simulations based on the CE Hamiltonian for the (Th1xCex)O2 system we have investigated the contribution of SRO to the mixing enthalpies for temperatures from 400 to 600 K, finding that they lead to at most a 5% lowering of the magnitude of the calculated mixing energies over this temperature range. While these results thus suggest that SRO contributions are relatively small (at least for the (Th1xCex)O2 system), it is important to keep in mind that due to the neglect of SRO the calculated results should be viewed as an upper bound to DHmix values for systems equilibrated at finite temperatures. The CE, SQS and calorimetry results for DHmix for (Th0.5Ce0.5)O2 and (Zr0.5Ce0.5)O2 are listed in the second, fourth and seventh columns of Table 1, respectively. Although the calculated results predict the trends in the experimental data very well, namely the

1 CeO2

16

CE: ordered structures CE: random solid solutions SQS: random solid solutions

14 12 10 8 6 4 2 0 0 ZrO2

0.2

0.4

0.6

Ce mol fraction

0.8

1 CeO2

Fig. 1. Cluster expansion (CE) and special quasirandom structure (SQS) predictions of formation enthalpy for (ThxCe1x) O2 (left) and (ZrxCe1x) O2 (right) ordered structures and random solid solutions from constituent fluorite-structured oxides.

295

V. Alexandrov et al. / Journal of Nuclear Materials 444 (2014) 292–297

Table 1 Mixing enthalpies (kJ/mol-cation) predicted by cluster expansion (CE), special quasirandom structure (SQS) approach, elasticity theory and derived by modeling of hightemperature oxide melt solution calorimetric data for the limit of the random solid solution at 50% oxide mixing level. SQS results are obtained using a 96-atom supercell, while elasticity theory predictions were made using both calculated and experimental lattice (a) and elastic (Cij) constants. For the SQS results, the ‘‘V-relaxed’’ column lists calculated results in which the ionic positions are held fixed at ideal fluorite positions and the energies are relaxed with respect to volume only; the ‘‘Fully-relaxed’’ column lists results for which the energy is optimized with respect to volume, cell shape and ionic positions. CE

(Zr0.5Ce0.5) O2 (Th0.5Ce0.5) O2

9.2 1.8

SQS-96

Elasticity theory

Calorimetry

V-relaxed

Fully-relaxed

Calc. a, Cij

Expt. a, Cij

22.0 3.0

9.6 1.7

13.6 2.1

10.7 3.8

positive and symmetric nature of the mixing energies, as well as the larger magnitudes for DHmix in the (Zr1xCex)O2 system relative to (Th1xCex)O2, there are quantitative discrepancies in the results. Specifically, the CE and SQS underestimate the measured magnitudes of DHmix for (Zr0.5Ce0.5)O2 by 3.6 and 3.2 kJ/mol, respectively. For (Th0.5Ce0.5)O2 the discrepancy between measurements and calculations is 1.9 kJ/mol and 2.0 kJ/mol for CE and SQS results, respectively. While these discrepancies are comparable to the error bars reported for the fits of the calorimetry data, the trend for underestimation of the magnitudes of DHmix in the calculations relative to the measurements is consistent across the entire composition range in both systems. In Ref. [4] it was argued that the discrepancy between experiment and calculations for (Th1xCex)O2 stems from an underestimation of magnitude of the elastic constants by GGA for the constituent oxides. This analysis is discussed in further detail in the next section, where a similar analysis will be presented for the (Zr1xCex)O2 system. It is also interesting to compare the results obtained from the present work with previously published DFT–GGA calculations [5,6] based on an approach for modeling configurationally disordered solid solutions that is similar to the CE described above, but with a Hamiltonian that contains shorter-ranged interactions and is based on a more limited set of calculated input structural energies. The previous calculations are in good agreement with our modeling results exhibiting DHmix values that are positive and symmetric about x = 0.5 for both (Th1xCex)O2 and (Zr1xCex)O2 solid solutions. Quantitatively, DHmix was found to be slightly higher than our CE values for both solid solutions yielding the maximum DHmix of 2.7 kJ/mol-cation for (Th1xCex)O2 and 11.8 kJ/molcation for (Zr1xCex)O2.

3.2. Elasticity theory predictions for mixing energies For the homovalent ionic solid solutions considered in this work, we expect that a dominant contribution to DHmix arises from the elastic strain energy originating from the difference in cation ionic radii. The importance of elastic-strain-energy in governing the mixing energetics of solid solutions has been discussed in detail by several authors previously (e.g., Refs. [10,27–29]). For the systems considered in the present work, this contribution arises due to the fact that the constituent cations prefer different metal–oxygen bond lengths, and thus the formation of a crystalline solid solution inevitably involves a positive contribution to the energy arising from bond strains. This strain energy can be thought of as originating from two contributions: (i) the strain energy required to compress or expand each metal–oxygen bond to the length dictated by the lattice parameter of an undistorted fluorite crystal structure, and (ii) the relaxation energy associated with atomic displacements away from the sites of the perfect fluorite structure, which serve to partially relieve the strain energy arising in (i). The relative magnitudes of contributions (i and ii) can be estimated from the difference in energy between the volume-relaxed and fully-relaxed SQS mixing energies given in the third

12.8 ± 2.0 Ref. [7] 3.7 ± 2.5 Ref. [4]

and fourth columns of Table 1; these results show that the relaxation energy reduces the magnitude of DHmix by roughly one-half. To estimate the strain-energy contribution (Eelast) to DHmix for (Th1xCex)O2 and (Zr1xCex)O2 we make use of the anisotropic elasticity theory described in Section 2.3. The values of Eelast obtained using DFT-GGA calculated lattice constants and elastic constants (see Table 2) are listed in the fifth column of Table 1. For (Th0.5Ce0.5)O2 the magnitude of Eelast shows excellent agreement (to within 0.4 kJ/mol) with the CE and (fully-relaxed) SQS results. Similarly, the value of Eelast computed using experimentally measured elastic constants and lattice parameters (see column six in Table 1) shows excellent agreement with the calorimetry measurements. A comparison of the values of Eelast derived from calculated and measured lattice and elastic constants for (Th0.5Ce0.5)O2 suggests that the discrepancy between calculations and measured values of DHmix originates from an underestimation of the elastic moduli in the DFT–GGA calculations, as discussed in Ref. [4]. More, specifically, the results in Table 2 show that the discrepancy likely arises primarily from the relatively large underestimation of C11 and C44 by the DFT–GGA calculations for CeO2. For (Zr1xCex)O2 the value of Eelast derived from calculated elastic moduli and lattice constants is seen to significantly overestimate the values derived from CE and SQS. This discrepancy may reflect the importance of anharmonic contributions to the elastic energy not considered in the linear-elasticity theory underlying Eqs. (2) and (3). Specifically, the cubic fluorite structure of ZrO2 is known to be dynamically unstable (i.e., it has imaginary phonon frequencies) at zero temperature in DFT calculations [39], such that the force constants in the compounds considered for the parameterization of the CE may be softer than expected based on the estimate from the elastic constants underlying Eqs. (2) and (3). Additionally, (Zr1xCex)O2 features a much larger cation size mismatch, relative to (Th1xCex)O2. Both effects are expected to give rise to relatively larger displacements in (Zr1xCex)O2, as observed in the relaxed SQS structure at x = 0.5: the relaxed bond lengths in (Zr0.5Ce0.5)O2 are distorted by roughly 4% relative to the average value imposed by the fluorite structure, as compared to a much smaller value of approximately 0.4% in (Th0.5Ce0.5)O2. The larger bond strains and potentially soft force constants in (Zr1xCex)O2 are thus expected to lead to larger anharmonic contributions to the elastic energy in this system relative to (Th1xCex)O2. We note that if the size mismatch between the cations is increased further, the bond distortions become so large that the fluorite structure is no longer stable. Specifically, in analogous SQS calculations for (Th1xZrx)O2, where the cation size mismatch is Dr=r ¼ 22% (compared to 8% and 14% for (Th1xCex)O2 and (Ce1xZrx)O2, respectively) the bond distortions in the relaxed structures are as large as 15% and the nearest-neighbor coordination number for the cations deviates from the value of 8 characteristic of the fluorite structure. Such large distortions and changes in the coordination numbers are not found in the calculations for (Th1xCex)O2 and (Ce1xZrx)O2 using the PBE–GGA functional. We do note, however, that the inclusion of Hubbard-U corrections for Ce 4f orbitals (using U = 4 eV) increased the magnitude of the atomic

296

V. Alexandrov et al. / Journal of Nuclear Materials 444 (2014) 292–297

Table 2 Calculated lattice parameter a (Å), bulk modulus B (GPa), and elastic constants C11, C12, C44 (GPa) for dioxide compounds in the cubic fluorite structure. Compound

a

B

C11

C12

C44

Method

Reference

ThO2

5.621 5.61 5.597a 5.466 5.43 5.411b 5.45 5.50 5.466 5.484 5.126 5.111

187.3 198 – 173.2 184 – 222.4 200 197 140.1 250.1 194

343.7 376.0 367 317.2 353.8 403 389.3 363.6 256.5 250.4 505.8 417

109.1 109.8 106 101.2 99.1 105 138.9 118.8 167.9 87.0 122.2 82

73.9 68.1 79.7 50.7 51.0 60 71.3 57.4 59.2 55.3 63.1 47

PBE PBE (FP-LMTO) Ultrasound velocity expt. PBE PBE (FP-LMTO) Brillouin scattering expt. LDA+U (U = 4) GGA+U (U = 4) GGA+U (U = 4) GGA+U (U = 4) PBE Ultrasonic velocity expt.c

This work Ref. [30] Ref. [32] This work Ref. [30] Ref. [34] Ref. [35] Ref. [36] Ref. [35] Ref. [37] This work Ref. [38]

CeO2

UO2 NpO2 PuO2 AmO2 ZrO2 ZrO2 a b c

X-ray diffraction Ref. [31]. X-ray diffraction Ref. [33]. Lattice and elastic constants of pure cubic ZrO2 were determined for single crystals of yttria-stabilized zirconia by extrapolating to 0 mol% Y2O3.

displacements in relaxation calculations for a SQS model of (Ce0.5Zr0.5)O2 and led to changes in coordination numbers similar to what was found for (Th1xZrx)O2. For (Zr1xCex)O2 the measured values of DHmix and the values of Eelast derived from measured elastic moduli (which for ZrO2 are based on an extrapolation from the region of stability of the cubic phase in Y2O3–ZrO2 solid solutions [38]) agree to within 2.1 kJ/mol, which is close to the quoted uncertainties of the regular-solutionmodel fits to the experimental data (2.0 kJ/mol). Thus, the level of agreement between the elasticity theory, DFT calculations, and measurements can be considered as very good. Unlike the situation for the (Th1xCex)O2, it is difficult to assess the origin of the discrepancy between DFT calculations and experimental measurements, as the DFT–GGA results underestimate the elastic constants in CeO2 while overestimating those in ZrO2. One possibility is that the instability of the short-wavelength phonons for ZrO2 in the DFT calculations may lead to an overestimation of the elastic relaxation energy in both the CE and SQS calculations. Despite the discrepancies noted in the previous paragraph, the magnitudes of Eelast are clearly comparable to the calculated and measured values of DHmix for each of the cases presented in Table 1. The results thus support the interpretation that a dominant contribution to the mixing energy arises from strain energy for both of the systems considered in this study. Further, the above analysis suggests that the linear elasticity theory provides more accurate results for the (ThxCe1x)O2 system where the bond strains are smaller in magnitude, and the end member compounds are dynamically stable in the cubic fluorite structure. Given the accuracy of the elasticity theory in predicting the mixing enthalpies in the (ThxCe1x)O2 system, we further employ this formalism to estimate the elastic contributions to the mixing energies in a series of related fluorite-based solid solutions which

feature comparable cation size mismatch, and which are of direct relevance in the context of understanding the thermochemistry of nuclear fuel materials. The results of this analysis are presented in Fig. 2, which plots the elasticity-theory predictions of mixing enthalpies for (Th1xMx)O2 and (U1xMx)O2 fluorite-structured actinide dioxide solid solutions, with M = U, Th, Ce, Np, Pu and Am. These results are based on Eqs. (2) and (3) using the previously-published calculated structural and elastic parameters listed in Table 2. It is seen that the mixing energies are predicted to be largest in magnitude for the (Th1xUx)O2 system, with the values for the technologically important (U1xPux)O2 being more than an order of magnitude smaller. Assuming that the elastic energies given in Fig. 2 represent a dominant contribution to the DHmix, they can be used within the framework of mean-field theory to estimate the critical temperatures for the (incoherent) miscibility gaps in these systems. We note that a comparison of mean-field theory and Monte-Carlo thermodynamic-integration calculations was discussed in Ref. [4] for the (Th1xCex)O2 system. The mean-field predictions based on the results in Fig. 2 suggest maximum critical temperatures that are less than 450 °C, for (Th1xCex)O2, and considerably smaller for all of the other solid solutions investigated. This critical temperature for (Th1xCex)O2 is in reasonable agreement with the previous calculations [6] predicting a miscibility gap below 350 °C. Overall, these results thus suggest that for the temperatures relevant to the operation of nuclear fuels homogeneous solid solutions of the systems considered in Fig. 2 are expected to be thermodynamically stable with respect to phase separation. We end this section by noting that the above predictions for solid solutions containing U, Ce and Pu assume that these cations adopt 4+ oxidation states. As discussed in detail in Refs. [40,41] charge transfer between U and Ce can lead to the formation of solid

Fig. 2. Microscopic elasticity theory predictions of the elastic contribution to the mixing enthalpies for ThO2–MO2 (left) and UO2–MO2 (right) solid solutions.

V. Alexandrov et al. / Journal of Nuclear Materials 444 (2014) 292–297

solutions containing a mixture of U5+, U4+, Ce4+ and Ce3+ cations in (U1xCex)O2 solid solutions at finite temperature. Similar effects may arise also in the (U1xPux)O2 system. The formation of solid solutions with mixed cation charge states will lower the mixing free energy, and thus lead to lower values for the miscibility-gap critical temperatures than would be estimated from mean-field theory using the elastic contributions to the mixing energies plotted in Fig. 2. 4. Summary The mixing energetics of fluorite-structured oxide solid solutions have been computed employing DFT–GGA based clusterexpansion (CE) and special-quasirandom-structure (SQS) approaches, combined with continuum-elasticity analyzes. The DFT–GGA results for DHmix obtained by SQS and CE methods show excellent quantitative agreement for random solid solutions in both the (Th1xCex)O2 and (Zr1xCex)O2 systems. Consistent with calorimetry measurements, these results feature DHmix values that are positive at all compositions, approximately symmetric about x = 0.5, and larger in magnitude for (Zr1xCex)O2 relative to (Th1xCex)O2. Quantitatively, the DFT-GGA results show a tendency for underestimating the magnitude of the calorimetry results by approximately 2–3 kJ/mol. A comparison of the CE and SQS results with linear elasticity theory suggests that the dominant contribution to the mixing energy arises from the strain energy associated with cation size mismatch. We thus employ continuum theory to estimate the elastic energy contribution to DHmix for (Th1xMx)O2 and (U1xMx)O2 actinide-dioxide solid solutions, with M = U, Th, Ce, Np, Pu and Am. The maximum values for DHmix obtained in these calculations is 3 kJ/mol for (Th1xUx)O2, with values for all other systems being considerably smaller in magnitude. 5. Acknowledgements This work was supported as part of the Materials Science of Actinides, an Energy Frontier Research Center funded by the US Department of Energy, Office of Science, Office of Basic Energy Sciences under award number DE-SC0001089. This work made use of resources of the National Energy Research Scientific Computing Center, supported by the Office of Science of the US Department of Energy under Contract No. DE-AC02-05CH11231. References [1] R. Devanathan, L. Van Brutzel, A. Chartier, C. Gueneau, A.E. Mattsson, V. Tikare, T. Bartel, T. Besmann, M. Stan, P. Van Uffelen, Energy & Environmental Science 3 (10) (2010) 1406–1426. http://dx.doi.org/10.1039/c0ee00028k. [2] C. Degueldre, J. Bertsch, G. Kuri, M. Martin, Energy & Environmental Science 4 (5) (2011) 1651–1661, http://dx.doi.org/10.1039/c0ee00476f. [3] V. Alexandrov, N. Grønbech-Jensen, A. Navrotsky, M. Asta, Physical Review B 82 (17) (2010) 174115, http://dx.doi.org/10.1103/PhysRevB.82.174115. [4] T.Y. Shvareva, V. Alexandrov, M. Asta, A. Navrotsky, Journal of Nuclear Materials 419 (2011) 72–75. [5] L.C. Shuller, R.C. Ewing, U. Becker, Journal of Nuclear Materials 412 (1) (2011) 13–21, http://dx.doi.org/10.1016/j.jnucmat.2011.01.017. [6] L.C. Shuller-Nickles, R.C. Ewing, U. Becker, Journal of Solid State Chemistry 197 (2013) 550–559. [7] T.A. Lee, C.R. Stanek, K.J. McClellan, J.N. Mitchell, A. Navrotsky, Journal of Materials Research 23 (4) (2008) 1105–1112, http://dx.doi.org/10.1557/ jmr.2008.0143. [8] C. Lombardi, L. Luzzi, E. Padovani, F. Vettraino, Progress in Nuclear Energy 50 (8) (2008) 944–953, http://dx.doi.org/10.1016/j.pnucene.2008.03.006.

297

[9] C. Degueldre, Journal of Alloys and Compounds 444 (2007) 36–41, http:// dx.doi.org/10.1016/j.jallcom.2006.11.203. [10] A.G. Khachaturyan, Theory of Structural Transformations in Solids, John Wiley & Sons, New York, 1983. [11] J.P. Perdew, K. Burke, M. Ernzerhof, Physical Review Letters 77 (18) (1996) 3865–3868, http://dx.doi.org/10.1103/PhysRevLett.77.3865. [12] S.L. Dudarev, G.A. Botton, S.Y. Savrasov, C.J. Humphreys, A.P. Sutton, Physical Review B 57 (3) (1998) 1505–1509, http://dx.doi.org/10.1103/ PhysRevB.57.1505. [13] P.E. Blöchl, Physical Review B 50 (24) (1994) 17953–17979, http://dx.doi.org/ 10.1103/PhysRevB.50.17953. [14] G. Kresse, J. Hafner, Physical Review B 47 (1) (1993) 558–561, http:// dx.doi.org/10.1103/PhysRevB.47.558. [15] G. Kresse, J. Hafner, Physical Review B 49 (20) (1994) 14251–14269, http:// dx.doi.org/10.1103/PhysRevB.49.14251. [16] G. Kresse, J. Furthmüller, Physical Review B 54 (16) (1996) 11169–11186, http://dx.doi.org/10.1103/PhysRevB.54.11169. [17] H.J. Monkhorst, J.D. Pack, Physical Review B 13 (12) (1976) 5188–5192, http:// dx.doi.org/10.1103/PhysRevB.13.5188. [18] P.E. Blöchl, O. Jepsen, O.K. Andersen, Physical Review B 49 (23) (1994) 16223– 16233, http://dx.doi.org/10.1103/PhysRevB.49.16223. [19] M.J. Mehl, B.M. Klein, D.A. Papaconstantopoulos, First principles calculations of elastic properties of metals, intermetallic compounds, Principles and Practice, Vol. 1, John Wiley & Sons, London, 1994 (Chapter 9). [20] A. van de Walle, M. Asta, G. Ceder, Calphad-Computer Coupling of Phase Diagrams and Thermochemistry 26 (4) (2002) 539–553, http://dx.doi.org/ 10.1016/S0364-5916(02)80006-2. [21] A. van de Walle, Calphad-Computer Coupling of Phase Diagrams and Thermochemistry 33 (2) (2009) 266–278, http://dx.doi.org/10.1016/ j.calphad.2008.12.005. [22] J.M. Sanchez, F. Ducastelle, D. Gratias, Physica A 128 (1–2) (1984) 334–350, http://dx.doi.org/10.1016/0378-4371(84)90096-7. [23] D. deFontaine, Solid State Physics – Advances in Research and Applications 47 (47) (1994) 33–176. [24] S.H. Wei, L.G. Ferreira, J.E. Bernard, A. Zunger, Physical Review B 42 (15) (1990) 9622–9649, http://dx.doi.org/10.1103/PhysRevB.42.9622. [25] J. von Pezold, A. Dick, M. Friak, J. Neugebauer, Physical Review B 81 (9) (2010) 094203, http://dx.doi.org/10.1103/PhysRevB.81.094203. [26] Y. Du, M. Yashima, T. Koura, M. Kakihana, M. Yoshimura, Scripta Metallurgica Et Materialia 31 (3) (1994) 327–332, http://dx.doi.org/10.1016/0956716X(94)90291-7. [27] H. Yamauchi, D. deFontaine, Acta Metallurgica 27 (5) (1979) 763–776, http:// dx.doi.org/10.1016/0001-6160(79)90110-X. [28] D.B. Laks, L.G. Ferreira, S. Froyen, A. Zunger, Physical Review B 46 (19) (1992) 12587–12605, http://dx.doi.org/10.1103/PhysRevB.46.12587. [29] J.W. Doak, C. Wolverton, Physical Review B 86 (14) (2012) 144202, http:// dx.doi.org/10.1103/PhysRevB.86.144202. [30] V. Kanchana, G. Vaitheeswaran, A. Svane, A. Delin, Journal of PhysicsCondensed Matter 18 (42) (2006) 9615–9624, http://dx.doi.org/10.1088/ 0953-8984/18/42/008. [31] J.S. Olsen, L. Gerward, V. Kanchana, G. Vaitheeswaran, Journal of Alloys and Compounds 381 (1–2) (2004) 37–40, http://dx.doi.org/10.1016/ j.jallcom.2004.04.099. [32] P.M. Macedo, W. Capps, J.B. Wachtman, Journal of the American Ceramic Society 47 (12) (1964), http://dx.doi.org/10.1111/j.1151-2916.1964.tb13130.x. 651. [33] L. Gerward, J.S. Olsen, L. Petit, G. Vaitheeswaran, V. Kanchana, A. Svane, Journal of Alloys and Compounds 400 (1–2) (2005) 56–61, http://dx.doi.org/10.1016/ j.jallcom.2005.04.008. [34] A. Nakajima, A. Yoshihara, M. Ishigame, Physical Review B 50 (18) (1994) 13297–13307, http://dx.doi.org/10.1103/PhysRevB.50.13297. [35] P. Zhang, B.T. Wang, X.G. Zhao, Physical Review B 82 (14) (2010) 144110, http://dx.doi.org/10.1103/PhysRevB.82.144110. [36] B.T. Wang, H.L. Shi, W.D. Li, P. Zhang, Physical Review B 81 (4) (2010) 045119, http://dx.doi.org/10.1103/PhysRevB.81.045119. [37] Y. Lu, Y. Yang, F. Zheng, B.T. Wang, P. Zhang, Electronic, Mechanical, and Thermodynamic Properties of Americium Dioxide, arXiv:1211.0378 [condmat.mtrl-sci]. [38] H.M. Kandil, J.D. Greiner, J.F. Smith, Journal of the American Ceramic Society 67 (5) (1984) 341–346, http://dx.doi.org/10.1111/j.1151-2916.1984.tb19534.x. [39] M. Sternik, K. Parlinski, Journal of Chemical Physics 122 (6) (2005) 064707, http://dx.doi.org/10.1063/1.1849157. [40] B.E. Hanken, C.R. Stanek, N. Grønbech-Jensen, M. Asta, Physical Review B 84 (8) (2011) 085131, http://dx.doi.org/10.1103/PhysRevB.84.085131. [41] B.E. Hanken, T.Y. Shvareva, N. Grønbech-Jensen, C.R. Stanek, M. Asta, A. Navrotsky, Physical Chemistry Chemical Physics 14 (16) (2012) 5680–5685, http://dx.doi.org/10.1039/c2cp40295e.