Atomistic and experimental study on thermal conductivity of bulk and porous cerium dioxide.

Cerium dioxide (CeO2) is a surrogate material for traditional nuclear fuels and an essential material for a wide variety of industrial applications both in its bulk and nanometer length scale. Despite this fact, the underlying physics of thermal conductivity (kL), a crucial design parameter in industrial applications, has not received enough attention. In this article, a systematic investigation of the phonon transport properties was performed using ab initio calculations unified with the Boltzmann transport equation. An extensive examination of the phonon mode contribution, available three-phonon scattering phase space, mode Grüneisen parameter and mean free path (MFP) distributions were also conducted. To further augment theoretical predictions of the kL, measurements were made on specimens prepared by spark plasma sintering using the laser flash technique. Since the sample porosity plays a vital role in the value of measured kL, the effect of porosity on kL by molecular dynamics (MD) simulations were investigated. Finally, we also determined the nanostructuring effect on the thermal properties of CeO2. Since CeO2 films find application in various industries, the dependence of thickness on the in-plane and cross-plane kL for an infinite CeO2 thin film was also reported.

Introduction

The lattice thermal conductivity (k) depends on the mean free path (MFP), i.e. the distance travelled by the phonons before the occurrence of any scattering events. A detailed quantitative understanding of the mode dependent phonon properties and the bulk MFPs in each material is imperative to exercise control over its thermal management. Furthermore, the thermal transport properties of the bulk and the nanostructured (wires, films, porous, and nanocrystalline material) counterparts of a given material differ significantly due to increased boundary scattering in the latter. Understanding this phenomenon becomes necessary when studying thermal conductivity of materials such as Cerium dioxide (CeO2), which finds applications in both bulk and nanostructured form. Bulk CeO2 is frequently used as a surrogate material in the nuclear industry, as an oxygen storage material in automobile exhaust systems[1], and as an electrolyte in solid oxide fuel cells[2]. Also, nanostructured CeO2 is used for diverse applications in various industries. For instance, in electronics, the thin films of CeO2 are among the most prospective buffer layers for high-temperature semiconductor and ferroelectric films deposited on silicon[3,4]. In the nuclear industry, CeO2 thin film deposits have been proposed for controlling the oxidation of nickel alloys[5], and microspheres of CeO2 have been synthesized for space nuclear applications[6]. Unfortunately, to date, accurate atomic-scale investigation on the details of mode-wise thermal transport properties and MFP of CeO2 have not been reported.

Although several studies have reported the k of CeO2, the value remains ambiguous. Experimentally, k of CeO2 has been determined using laser flash technique. Nelson et al.[7] and Khafizov et al.[8] measured the k of 95% dense CeO2 pellet and reported respective values of 6.6 W/mK and 17.5 W/mK (correcting for porosity using a modified Loeb expression), at room temperature. The ab initio theoretical calculations for the k of CeO2 performed by Xiao et al.[9] and the molecular dynamics (MD) study by Khafizov et al.[8] predict the k as 12 W/mK and 19 W/mK respectively. However, Xiao et al.[9] used the analytical Slack model[10], which heavily relied on the fitting parameters that may obscure the underlying physics of heat transfer. Whereas, the MD simulations performed by Khafizov et al.[8] could have over-predicted the k at lower temperatures as MD simulations do not account for intrinsic scattering at these temperatures. Right through the literature, the value of k of CeO2 at room temperature has varied between (6.6 W/mK and 19 W/mK) and thus exists a disparity. To resolve this aberration, we perform an ab initio calculation using density functional theory (DFT)[11], density functional perturbation theory (DFPT)[12] and lattice dynamics in concord with the Boltzmann transport equation (BTE) using ShengBTE[13]. This theoretical approach does not require any assumption on the phonon lifetime to predict the k. BTE has been successfully utilized in the investigations of phonon transport properties for many materials[14-19] with high accuracy. In addition to the ab initio calculations, MD simulations were carried out to understand the effect of porosity on the k of CeO2. Moreover, using the laser flash experiment technique, the thermal conductivity of porous CeO2 pellets prepared by spark plasma sintering (SPS) method has also been reported be a restrictive factor for the heat transfer.

In nanostructured materials, where the phonon MFP is comparable to the grain size, grain boundaries can modulated thermoreflectance technique and an analytical solution of BTE to understand the thermal transport properties of ceria thin films grown by unbalanced magnetron sputtering. The significantly reduced conductivity of these thin films compared to the bulk CeO2 was attributed to the combined effect of the point defects, grain boundaries, and dislocations. However, their study was not able to shed light on the influence of nanostructuring. For nanoscale applications, the characteristic lengths of the nanostructures can be a limiting factor in the thermal transport. To have a comprehensive understanding of this limiting factor, it is critical to know MFP of CeO2. It is also important to know what specific length of the nanostructure is going to be significant for CeO2. Therefore, in this work, the ab initio prediction of the MFP, the relaxation time of the phonons, their mode-wise contribution to thermal conductivity, and the effect of nanostructuring on the reduction of thermal conductivity were investigated. These findings will aid the selection of the size and thickness of the nanostructures in tuning the thermal properties. Furthermore, the effect of nanostructuring by studying the cross plane and in-plane thermal conductivity of CeO2 thin films and the impact of the thickness of the thin films and temperature on the thermal conductivity was also performed. These predicted results not only enable the accurate explanation of the experimental results but also guide further designs and applications.

Results and Discussion

Crystal Structure and elastic constants

CeO2 has a fluorite structure with three independent atoms per unit cell and belongs to the space group of Fm-3m (225), as shown in Fig. 1. The equilibrium lattice constants were obtained by minimizing the total energy with respect to the lattice parameter and atomic positions. The detailed description of the optimized parameter and the pseudopotential used for the simulation were provided in section 4.1. Table 1. presents the calculated lattice constant (a), bulk modulus (B) and stiffness constants (C) of CeO2 in comparison with the values from the previous DFT calculations[20-23] and the experiment[24]. The PBEsol functional predicted the lattice constant within an error less than 0.01% compared with the experimental data reported by Nakajima et al.[24] As expected, the new PBEsol functional was able to quantify the ground state structural properties more accurately than previous DFT calculations performed using linear density approximation (LDA)[20], GGA[21], Hartwigsen- Goedecker-Hutter (HGH)[22] and PAW[23] pseudopotentials.

Figure 1

Crystal structure of CeO2 in the conventional cubic unit cell (blue and red represent Ce and O atoms respectively.

Table 1

Structural and mechanical properties of CeO2 compared with previous DFT calculations and the experiment.

a (Å)	B (GPa)	C11 (GPa)	C12 (GPa)	C44 (GPa)	G (GPa)	Y (GPa)	η	References
5.409	192	355	110	62	82	215	0.313	This work (PBEsol)
5.37	203	371	117	68	—	—	—	LDA[20]
5.45	193	—	—	—	—	—	—	PAW[23]
5.37	210	386	124	73	—	—	—	HGH[22]
5.44	181	347	97	51	74	195	0.218	GGA[21]
5.41	204	403	105	60	—	—	—	Exp.[24]

The elastic constants are important parameters that provide information on the properties of a material such as stiffness, strength, mechanical stability, hardness, and ductility or brittleness[25]. We evaluated the single crystal stiffness constants of CeO2 by using a stress-strain method[26,27] with the help of our in-house interface qe-nipy-advanced[28]. The stiffness constants for a mechanically stable cubic structure should satisfy the following Born’s mechanical stability criteria[29] . The listed stiffness constants have satisfied these criteria indicating that the system is mechanically stable.

From the calculated stiffness constants, the polycrystalline bulk modulus, shear modulus, Young’s modulus and Poisson’s ratio (listed in Table 1) were determined using the Voigt-Reuss-Hill approach[30-32]. From Table 1. the bulk modulus values (192 GPa) at zero temperature obtained from PBEsol method is only about 6% lower than the experimental value of 204 GPa[24] at room temperature. Since the calculated structural and mechanical properties are in excellent agreement with the experimental data, the same PBEsol functional pseudopotentials were used for the study of phonon properties and the k of CeO2.

Lattice dynamics

Lattice dynamics is critical for understanding the thermal properties of crystalline solids at finite temperatures. For instance, the phonon densities of states are required to evaluate thermodynamic properties such as thermal expansion coefficient (α), C, entropy (S) and k. Moreover, the fundamental reasons for unique thermal characteristics of a material can be ascertained by analyzing the phonon scattering mechanism through phonon group velocities, and phonon mean free path, relaxation time and phonon scattering phase space[33]. At present, two approaches the linear-response approach[34] and the direct approach[35,36] are widely used to evaluate the phonon dispersion. In Fig. 2a the calculated phonon dispersion spectra using these two approaches the linear-response approach (QE) and the direct approach (Phonopy) are presented along the Γ-X-Γ-L high symmetry points in the Brillouin zone and compared with the known experimental data.

Figure 2

(a) Phonon spectra of CeO2 predicted by the linear response approach using q-point of 6 × 6 × 6 (red line) and the finite displacement (FD) methods (dotted green line) compared with the experiment done by Clausen K. et al.[37] (blue dots). (b) Partial (dotted lines) and total phonon density of states (solid black line).

The predicted phonon spectra using both the linear-response approach and the finite displacement (FD) methods are in good agreement with the experimental data generated by Clausen K et al.[37], for the entire Brillouin zone using inelastic neutron scattering at 293 K. CeO2 has three (n) atoms per primitive unit cell and thus there are in three dimensions (d = 3) three acoustic mode phonons and d (n − 1) = 6 optical phonons. The group theory analysis gives the decomposition of the zone centre modes as doubly degenerate transverse optical (TO) mode, the triply degenerate Raman active mode and the non-degenerate longitudinal optical (LO) mode. The three-zone center frequencies are 8.42 (8.15), 13.03 (13.94), and 15.93(17.84) THz respectively (shown in brackets are the experimental value from the ref.[38]). The partial lifting of degeneracy between the LO and TO phonons at the Brillouin zone centre is due to the polarization effect, which indicates that the CeO2 is a polar material as shown in Fig. 2a The Born effective charge for an insulator is a measure of the change in electronic polarization due to ionic displacements. These charges are essential for elucidating the physical understanding of piezoelectric and ferroelectric properties since they describe the coupling between lattice displacements and the electric field. The Born effective charges,  = 5.793 (5.746)  = 2.896 (2.873), agrees well with the values reported (shown in brackets) by Xiao et al.[9] The calculated dielectric constant value 7.87 (6.0) is in reasonable agreement the work done by Santha et al.[38].

The partial and total phonon density of states is shown in Fig. 2b. The partial density of the state shows that the frequency vibrations lower than 6.0 THz are dominated by the Ce ions and the higher frequency vibrations are mainly contributed by the dynamics of O ions. The phonon dispersion of CeO2 along Γ-X-Γ-L, high symmetry points in the Brillouin zone does not display frequency gap between optical and acoustic phonons. The absence of frequency gap in case of CeO2 indicates that three-phonon scattering is dominant resulting in relatively lower phonon relaxation time thereby leading to lower k. The accurate prediction of phonon dispersion curve of CeO2 using both the linear-response approach and direct approach enables reasonable precision in the projection of phonon-assisted properties.

Three-phonon scattering phase space and Grüneisen parameter

Three-phonon scattering phase space (P3) quantitatively describes the number of scattering channels available for a phonon being scattered[39]. The available P3 gives an insight into the k of a material i.e. larger available P3, will have more channels for scattering and subsequently lower k. The P3 is solely determined from a material’s phonon spectra and is defined as[33],where q is the momenta of phonons, j is the phonon branches, is the normalized factor equal to unrestricted phase space for each type of process. and are the absorption and emission processes respectively, defined by,where N is the total number of q points in the first Brillouin zone. Figure 3a shows the P3 of CeO2 for the acoustic and optical modes. To gain a quantitative understanding of available P3 on k, the total volume in phase space for three phonon processes (P3_total) for CeO2 is compared with our previously reported results of a relatively higher k material such as silicon carbide (420 W/mK at 300 K)[17]. The P3_total (in units of 1/rad/ps) predicted along the same Brillioun path is higher for CeO2 (4.18 × 10−3) than SiC (1.66 × 10−3); clearly indicating that CeO2 has a larger three phonons scattering phase space, more channels for scattering and hence lower thermal conductivity than SiC. Lindsay and Broido[33] have already established the phase space available for the three-phonon scattering for semiconductors. In this work, the possible P3 for a lanthanide oxide was illustrated, and these findings would aid in qualitatively determining the thermal conductivities of new materials with known phonon dispersion.

Figure 3

(a) Three-phonon scattering phase space. (b) Mode Grüneisen parameter (γqj).

Grüneisen parameter (y) illustrates the anharmonicity in a crystal and is defined as the shift in phonon frequency with change in volume. It is an established fact that anharmonicity in an ordered crystal structure determines the strength of each scattering channel and the efficiency of each phonon mode in heat conduction is inversely proportional to square of Grüneisen parameter. To quantitatively examine the anharmonicity of CeO2, we plot the mode Grüneisen parameters at zero temperature (as shown in Fig. 3b) derived from the third-order anharmonic force constants (shown in Eq. (3)) as implemented in the ShengBTE[13] and is defined as:where are the third-order force constants, e is the phonon eigenvectors, n denotes the n primitive cell in the supercell, α is the Cartesian components of x, y, or, z, M refers to the mass of the atomic type n in the primitive cell and r is the position vector of the η atom in the l primitive cell. According to Fig. 3b, the sum of all mode Grüneisen parameters are positive indicating volume expansion. At zero temperature, the average Grüneisen coefficient over the whole frequency range was 2.07. The mode Grüneisen parameter near 7 THz (near the TO mode) is higher, indicating that the scattering of acoustic phonons by the optical mode caused hindrance to heat conduction. Changes in mode Grüneisen parameter with different q-points in the second order force constant and the neighboring atoms in the third order force constant is given as Supplementary Information (S.I) in S.I.1 (Fig. S1).

Theoretical prediction of thermal conductivity and mode contribution

The full iterative solution of BTE for CeO2 is shown in Fig. 4a. As expected, the k of CeO2 decreases with increasing temperature because the phonon-phonon scattering dominates the k at high temperatures. The k at 300 K is 16.71 W/mK. Our predictions were made assuming a defect-free crystal and thus could potentially represent the upper limit of k. The previous experimental value of thermal conductivity for CeO2 appears to be scattered. For example, for ~95% dense pellets at 373 K, Jakub et al.[40] and Nelson et al.[7] reported a value of 9.58 W/mK and 6.2 W/mK, respectively. These scattered values of k reported in the literature are generally justified by pointing to microstructures. This may not be true under all circumstances. In these situations, an accurate theoretical prediction without the usage of fitting parameters can significantly help in arriving at a meaningful physical explanation. The highest reported experimental value for CeO2 for a ~95% dense pellet at room temperature is by Jakub et al.[40] (14.1 W/mK) and our theoretical predictions at 300 K is 16.71 W/mK (red dots) for 100% dense single crystal, are in excellent agreement. Furthermore, the difference in k values obtained through our theoretical prediction and the previously reported simulated value of 12 W/mK at 300 K by Xiao et al.[9] (using fitting parameters) clearly indicate the apparent disadvantages of using fitting parameters. We also noted that Xiao et al. obtained the second order force constant required for their calculation using the finite difference (FD) method. Our k predictions using the second order force constant from the FD method (green dot) has also underpredicted the value of k by 19% in comparison with the k predicted using the second order force constant obtained from DFPT method. Moreover, accurate theoretical predictions enable us to arrive at a rational expression for the porosity correction, as there is a multitude of expressions to choose from, depending on the material and temperature range.

Figure 4

(a) The temperature dependence of the phonon k in CeO2 using BTE and experiment compared to previous theoretical[9] and experimental[7] values of k. (b) The mode-wise k of CeO2 in comparison with the simulated results for UO2 by Pang et al.[42].

Furthermore, we probe the relative contribution of the acoustic and optical phonon modes to the total k of CeO2. Typically, for ceramic materials, there is a tendency to neglect the optical contribution owing to their consistently low group velocities and shorter lifetimes[41]. However, the Uranium dioxide (UO2) with a reported critical optical mode is an exception[42]. In the wake of this, it would be interesting to know the mode-wise thermal conductivity of the popular nuclear fuel surrogate CeO2. To best of author’s knowledge, the mode-dependent k of CeO2 has not received enough attention. Figure 4b shows the simulated contribution of mode-wise k of CeO2 compared with UO2 at room temperature previously calculated by Pang et al.[42]. In CeO2 the optical phonons contribute to about 27% of the total k compared to the acoustic modes. The optical contributions at 300 K can be precisely broken down to 13.1%, 13.7%, and ~1% from the doubly degenerate transverse optical mode, the triply degenerate Raman active mode and the non-degenerate longitudinal optical mode, respectively. The theory predicts that like UO2, its surrogate CeO2 will also have strong optical mode contribution to k. However, in the case of UO2, Pang et al.[42] have experimentally proved that the contribution of the TO branch to thermal conductivity is, in fact, higher than the theoretical prediction. Our theoretical predictions suggest an even more significant contribution of the TO branch towards thermal conductivity of CeO2, and could, therefore, make an experimental validation worthwhile.

Effect of porosity on the thermal conductivity using MD simulations and experiments

So far in this manuscript, the primary focus was in understanding the underlying physics of thermal transport and determining the k of a defect-free single crystal CeO2 using first principles approach. However, manufacturing a defect-free specimen of CeO2 with 100% theoretical density (TD) is impracticable, and most importantly k of a sample with porosity and defects are expected to be lower than 100% dense CeO2. Therefore, it is essential to determine the effect of porosity on k of CeO2 quantitatively. Hence, the effect of porosity on the k of CeO2 has been determined using both simulation (MD approach) and experiment. It has to be noted that evaluation of the effect of porosity using first principles is still in genesis, and for that reason, MD simulations were used. Moreover, the theoretical prediction of k using both BTE and MD approach complement each other well. The BTE calculates the quantum mechanical scattering rates directly by considering only the lowest order of phonon-phonon interactions and hence become less accurate at very high temperatures where higher order interactions become significant. In contrast, the MD simulations based on classical potential are less reliable at a lower temperature. However, at high temperatures, the phonon-phonon interactions of all orders are duly considered.

In this work, the k of CeO2 by MD simulations were calculated on systems of 8 × 8 × 8 (6144 atoms), 10 × 10 × 10 (12000 atoms), 12 × 12 × 12 (20736 atoms) and 20 × 20 × 20 (96000 atoms), using the Embedded Atom Many-body (EAM) potentials developed by Cooper et al.[43]. Table 2 shows that there is no noticeable difference in the k of CeO2 at 300 K with different supercell sizes, suggesting that 8 × 8 × 8 system is sufficient enough to represent all the phonon modes available to reproduce the phonon-phonon scattering present in the bulk CeO2[44]. For investigating the porosity influence on the k of CeO2, two cases with the cell sizes of 10 × 10 × 10 and 20 × 20 × 20 were considered. In each of the supercells, only one pore was introduced by manually removing 5% of atoms, in such a way that for every cerium atom two oxygen atoms are considered, to maintain the charge neutrality of the system. The supercell with vacancies are then made to relax keeping the cell volume same, and then MD calculations were carried out to evaluate the thermal conductivity of porous CeO2 (as explained in section 4.1).

Table 2

The size dependence of k of bulk CeO2 at a temperature of 300 K presented using EMD simulation and the Green-Kubo method.

Size of the supercell	Number of atoms	Thermal conductivity (W/mk) at 300 K
8 × 8 × 8	6144	18.87
10 × 10 × 10	12000	18.35
12 × 12 × 12	20736	18.35
20 × 20 × 20	96000	18.84

Figure 5a compares the experimental data to the predicted k values of CeO2 using MD for a pure crystalline CeO2 (black dots) and porous CeO2 (red dots) modelled on a 20 × 20 × 20 supercell. The k values predicted for pure crystalline CeO2 are in reasonable agreement with the experimental data (orange square dots) provided by Khafizov et al.[8]. As expected, the k of the porous CeO2 is lower than the pure crystalline CeO2. The measurements made by Nelson et al.[7] (navy blue dots) for a specimen of 95% TD was 6.22 W/mK at 100 °C, which is a reduction of almost 55% in the k value when compared to the pure crystalline CeO2. Suzuki et al.[45] (dark green square dot) investigated the reason for this unusual degradation in the k value of a 5% porous CeO2 and found that the purity of CeO2 is an essential factor (in work done by Nelson et al. cerium constituted only 98% of the metallic content of the feedstock). Furthermore, our EMD simulations performed for a 5% porous CeO2 (red dots) are in reasonable agreement with the work done by Suzuki et al.[45]. Besides, we also noted that k values predicted for porous CeO2 is significantly dependent on the size of supercell considered. Even though the porosity was modelled the same way in both the supercells, the k for the case of 10 × 10 × 10 is considerably lower than the case of 20 × 20 × 20. The Fig. 5b clearly shows that compared to the pure crystalline CeO2 at 300 K the k predicted for a porous CeO2 on a 10 × 10 × 10 supercell has reduced by 52%, whereas the same amount of porosity created on a 20 × 20 × 20 supercell decreased the k by 12%. The fact that the k is significantly lower for the 10x10x10 case and can be explained by the increased phonon scattering on the surface of the pores due to the decreased average distance between the pores, which directly can be compared to the mean free path of the phonons, as shown previously by Nichenko et al.[46]

Figure 5

(a) Simulated and experimental k of porous CeO2 as a function of temperature. (b) Influence of porosity on the k of CeO2, by considering the same porosity on a cell size of 10 × 10 × 10 (blue line) and 20 × 20 × 20 (red line) unit cells.

To validate the theoretical predictions and further expand the quantitative understanding of the effect of porosity on thermal conductivity, the experiments were carried out on the specimens prepared by SPS technique. The pellets of CeO2 with varying densities were made by controlling the sintering parameters such as the sintering temperature, pressure, hold time, and the sintering atmosphere. The more detailed description of the sintering conditions and its effect on density and the microstructure of CeO2 are given in our previous work[47]. Therefore for brevity, discussions would be limited to the samples on which the thermal conductivity measurements were made. In this work, the experiments were performed on CeO2 pellets prepared at a sintering temperature of 1000 °C and 1100 °C by maintaining the pressure and hold time at 50 MPa and 10 min respectively. However, for the benefit of the readers, it is important to note that as the sintering temperature increased (>1500 °C), the CeO2 reduced to Ce2O3. Moreover, in our previous work[47], it was observed that at a high sintering temperature (>1100 °C) and reductive atmosphere CeO2 exhibited a range of stoichiometry and such non-stoichiometric oxide had been susceptible to chemical expansion[48] leading to the disintegration of the sintered pellets. It was for this reason that the experimental determination of thermal conductivity was limited to two samples sintered at 1000 °C and 1100 °C.

These sintered pellets were characterized using XRD for the determination of the phase, and the XRD patterns revealed that the pellet sintered at 1000 °C (green line) and 1100 °C (pink line) have a face-centred cubic crystal structure (as illustrated in S.I.2 (Fig. S2)). The density of CeO2 was measured using the Archimedes principle, and the density of these pellets sintered at 1000 °C and 1100 °C were 70% and 75% respectively. As stated earlier in section 4.2, the thermal diffusivity can be determined using the Laser flash apparatus and Fig. 6 shows the variations of thermal diffusivity of these pellets as a function of temperature. Figures 5a and 6 respectively demonstrate that higher the density of the pellet, higher the thermal diffusivity and thermal conductivity. Moreover, we found that the density of the specimens played a vital role in the measured k values (as shown in Fig. 6) at lower temperatures (<1000 °C), however, at a higher temperature, the difference is less significant.

Figure 6

The thermal diffusivity of CeO2 as a function of temperature for specimens sintered at 1000 °C and 1100 °C. The uncertainty in the data point from the laser flash apparatus is ±4.5%.

Nanoscale size effect, cross-plane, and in-plane thermal conductivity

CeO2 nanostructures like wires, films, porous and nanocrystalline materials find applications in energy conversion, sensors, and microelectronics where tailoring of thermal transport property is essential. Moreover, a recent work evinces the capabilities of the ultra-thin CeO2 for oxygen storage[49]. It is established that the nanostructure surfaces significantly reduce the k compared to its bulk counterpart, due to the scattering of the energy carriers. The interplay between the characteristic length and the bulk MFP of the energy carriers is the fundamental physics that determines the dominance of the boundary scattering. The size of the crystalline domain, therefore, acts as limiting length for phonons MFP. Pertaining to this context, a detailed quantitative understanding of the MFP of CeO2 is certainly advantageous.

In Fig. 7 the dependence of normalized cumulative k on phonon MFP at 300 K are presented. The normalized k increases with the increase in MFP, with significant contributions to the k of CeO2 coming from the phonons with MFP between 1 nm to 100 nm. The contribution of phonons mode to k is uneven, and the phonons of MFP below 50 nm constitute 80% of k, indicating that curtailing the size below 50 nm can effectively reduce the k of CeO2. Therefore, this grasp of the MFP of CeO2 will aid the selection of sample size (thickness for thin film, diameter for nanowires and nanoparticles) for diverse technological applications that require a notable deviation from bulk thermal properties.

Figure 7

Thermal conductivity accumulation as a function of phonon MFP at room temperature.

To investigate the influence of nanostructuring on k of CeO2, we demonstrate the thickness dependence in the cross-plane and in-plane k of CeO2 along the (001) and (100) planes respectively at 300 K, as shown in Fig. 8a. For these predictions, we utilized almaBTE[50] a solver of the space-time dependent BTE for phonons in the structured material. The effective in-plane (‖) and cross-plane (⊥) thermal conductivity (k) in relaxation time approximation for a film of thickness L are evaluated as:where S is, the suppression function that considers the additional phonon scattering instigated by the film boundaries, C is the mode contribution to the volumetric heat capacity, v the group velocity, ∧ is the mean free path and ϑ is the angle between the group velocity and transport axis. The sum over q must be interpreted as the combination of a sum of branches and an average over the Brillouin zone. As anticipated, due to the boundary scattering the in-plane and cross-plane k of thin films reduced with the reduction in the thickness. If the thickness of CeO2 is slashed to 10 nm, its cross-plane (in-plane) k is only 4.96 (8.50) W/mK at 300 K. To verify the theoretical findings, we compared the results with experiment. It should be noted that theory predicts the k for the infinite two-dimensional CeO2 and therefore, the presence of defects would further scatter the phonons and reduce the k of CeO2 films. Khafizov et al.[8] reported k of 7.2 W/mK, for a film of long columnar grain size structure with an average column radius of ~450 nm and the film thickness of 12 μm.

Figure 8

(a) In-plane and cross-plane thermal conductivity of CeO2. (b) Cross-plane conductivity as a function of temperature.

Here, the thickness of 12 μm exceeds the condition that the characteristic size of CeO2 must be 50 nm to have a significant reduction in k. However, our simulation predicts a k of 13.0 W/mK for a thickness of 12 μm. The difference between theory and experiment can be primarily associated with the presence of defects and not to the effect of nanostructuring. Therefore, we can surmise that CeO2 thin films below ~50 nm with defects can reduce the k drastically. Metal oxides generally have large Seebeck coefficient at high temperatures and hence can be considered as candidate materials for advanced thermoelectric.

The coating of CeO2 on several alloys has improved the oxidation resistance and can be used in high-temperature applications in industries such as automobiles, aerospace and nuclear, where the knowledge of CeO2 cross-plane thermal conductivity as a function of coating thickness and temperature becomes vital. Therefore, in Fig. 8b we have presented the effect of temperature on the k of thin film and have observed that the k for thin films reduces considerably with an increase in temperature. The presented theoretical work can equip material scientists with vital information required for designing CeO2 thin films of optimal thickness that will aid in the design of experiments.

Summary

We have presented an extensive analysis of the k of CeO2 both in its bulk and nanoforms. The theoretical predictions of k using first principles unified with BTE not only provides an insight into the underlying physics of k of CeO2 but also helps explain the large discrepancy in the k value of CeO2 reported previously. The structural and mechanical properties of CeO2 could be predicted with better precision by the recently developed pseudopotential with PBEsol functional than any other previously reported DFT calculations. The phonon dispersions spectra of CeO2 is evaluated by both direct approach and linear response approach and could describe its polar nature while aptly showing good agreement with the experiment. An investigation of the available three-phonon scattering phase space and mode Grüneisen parameter of CeO2, reveals that the lower k of CeO2 is primarily due to increased scattering and strong anharmonicity respectively. Along with the theoretical investigation of k and its dependence on temperature, we also predict the mode wise contribution to the total k of CeO2. The analysis of mode wise k of CeO2 indicated the notable contribution to the overall thermal conductivity from the optical modes (~30%), akin to UO2. Additionally, we prepared CeO2 pellets of varying density by SPS technique and measured the thermal conductivity of these pellets using laser flash technique. The experimental analysis of the samples reveals the dependence of sintering parameters on the density of the CeO2 sample and its effect on its measured thermal conductivity. Moreover, we also conducted a theoretical study using the classical MD approach to corroborate the effect of porosity on thermal conductivity. Besides the detailed investigation of the thermal conductivity of CeO2 in its bulk form, this article also sheds light on the thermal transport property of nanostructured CeO2. We demonstrate the phonon MFP distribution of CeO2 is critical in the study of nanostructured materials and devices. The contribution of phonons modes to k is uneven, and the phonons of MFP below 65 nm constitute 80% of k, indicating that limiting the size below 65 nm can efficiently reduce the k of CeO2. The in-plane and cross-plane thermal conductivity of CeO2 thin film as a function of film thickness are also reported. These findings have an impact on heat dissipation in nanoelectronics and photonics, as well as the design of nanostructured thermoelectric materials with reduced thermal conductivity. To conclude, this work serves to moderate the existing ambiguity in the thermal conductivity value of CeO2 and provides practical information on CeO2 nanostructuring that will potentially meet the demands of numerous industrial applications.

Methods

Computational details

All the first principles calculations were performed using the DFT as implemented in the open source Quantum ESPRESSO[51] code. The pseudopotential used were norm conserved, and the electronic exchange-correlation is based on the generalized gradient approximation (GGA) with the Perdew, Burke, and Ernzerhof functional for solids (PBEsol)[52]. Geometry optimization evaluates the structural properties at zero Kelvin temperature by minimizing the total energy by varying both cell parameter and atom positions. We obtained the total energy convergence of CeO2, using an electron wave vector grid and the plane wave energy cutoff of 950 eV and 8 × 8 × 8 respectively. The criteria for the electronic energy convergence and the force convergence was respectively set a value of 10−12 eV and 10−7 eV/Å. The ground state structural properties were also assessed using the Birch-Murnaghan equation of state[53]. The elastic constants were calculated using the stress-strain method[26,27]. From the elastic constants, the mechanical properties such as the bulk modulus (B), shear modulus (G), Young’s modulus (Y), and Poisson’s ratio (η) were determined using the Voigt-Reuss-Hill averaging scheme[30-32]. Because we are analyzing a cubic, polycrystalline structure, therefore the elastic moduli can be evaluated assuming an isotropic aggregate of grains with non-isotropic elastic properties. For a cubic symmetry, there are three independent elastic constants: C11, C12, and C44. The bulk modulus for a cubic structure is the same for the Voigt, Reuss, and Hill averages, as shown in Eq. (5).Eq. (6) gives the shear modulus in the Voigt average:while Eq. (7) gives the Reuss average:

The arithmetic mean (Hill) G = (G + G)/2 is taken as shear modulus.

The Young’s modulus is calculated as:The Poisson’s ratio is calculated as:The k of CeO2 is predicted by solving the phonon BTE iteratively using the open source ShengBTE[13] code. The inputs required to solve the BTE are the second order force constants (harmonic force constant) and third order force constants (anharmonic force constants), which are defined as the second and third derivative of the potential energy (V) with respect to the atomic displacements respectively. The potential energy of the crystal with a unit cell characterized by a vector ‘l’ and the atomic positions in each unit cell described by the vector ‘b’, can be expanded in a Taylor series in power of the atomic displacement u (l, b). The first three terms of such a Taylor series expansion is given in Eq. (10) as stated in the ref.[54].The second order force constant () and the third order force constant ( are defined as shown in Eqs (11) and (12) respectively.The dynamical matrices for evaluating the phonon density of states is obtained from the second order force constant as follows,where q is the wave vector and m denote the mass of an atom at a site in the crystal. Further, diagonalizing the dynamical matrix yields the phonon frequencies (ω), which in turn provides the phonon group velocities and heat capacity at a fixed volume (C). In this work, the harmonic force constant is calculated using both linear approach (DFPT)[12] and direct approach (Parlinski-Li-Kawazoe method)[35], based on the supercell approach with finite displacement method as implemented in the Phonopy package[55]. To obtain the converged phonon properties, the calculations of the harmonic force constant were done using a q-point mesh of 8 × 8 × 8 (phonon dispersion for different q-points are shown as S.I.3, (Fig. S3)) and a 5 × 5 × 5 supercell of the primitive cell containing 375 atoms. The calculation of the third order force constants was performed on a 4 × 4 × 4 supercell, and the force cutoff distance was set to the ninth nearest neighboring atoms. The convergence of the k with respect to the number of q points used in the second order force constant, the number of neighboring atoms considered in third order force constant and the number of grid planes along each axis in the reciprocal space for solving the BTE is detailed in Fig. 9. From the cubic force constants, the phonon scattering processes are evaluated using Fermi’s golden rule, and finally, the k is calculated using the iterative solutions of the BTE as implemented in ShengBTE. In the approach implemented in ShengBTE, Eq. (14) gives the expression to compute the lattice thermal conductivity tensor,where Ω is the volume of the primitive cell, α and β are the Cartesian components of x, y, or z, kB is the Boltzmann constant, ω and υ are the angular frequency and group velocity respectively, τ is the relaxation time of mode λ and f0(ωλ) is the Bose-Einstein distribution function. The k presented in this work are the fully iterative solution of the Peierls-equation, and the convergence of k with the number of iterations starting from the zeroth-order approximation (which is equivalent to operating under RTA) and the k at low temperature (<300 K) are shown in S.I.4 (Fig. S4) and S.I.5 (Fig. S5) respectively.

Figure 9

Convergence of the k with respect to (a) the q-points considered for the second order force constants (by maintaining n = 9 and N = 18 × 18 × 18 for all the cases); (b) the number of neighboring atoms considered for calculating the third order force constant (by keeping q = 8 × 8 × 8 and N = 18 × 18 × 18 for all the cases); (c) the number of grid planes along each axis in reciprocal space for solving the BTE (by considering q = 8 × 8 × 8 and n = 9 for all the cases).

For a polar material, the interatomic forces are divided into two additive contributions; the analytic and non-analytic contributions (nac). The analytic contribution accounts for all the forces under the restricted periodic boundary conditions under which the averaged electric field is assumed to be zero. The nonanalytic contribution accounts the additional forces owing to non-zero averaged electric field[56]. The classical Newton’s second law of motion for describing the atomic vibrations for a polar solid is as shown in Eq. (15),where m denotes the atomic mass of the j atom in the primitive unit cell, α (α = x, y, z) is the component of the atomic displacement from its equilibrium position of the j atom in the P reference primitive unit cell within a supercell, as detailed in ref.[56]. The first term on the right-hand side of the Eq. (15) corresponds to the analytical force arising due to the short-range interatomic interactions. The U and E (U) in the first term represent the whole set of atomic displacements and the total energy respectively. In the second term, the dot product between the Born effective charge (Z) and the average of the electric field (E) induced by the atomic vibrations accounts for the non-analytic force due to the long-range coulombic interaction. In this work, CeO2 being a polar material, the non-analytical contribution was considered, and the Born charges and the dielectric constant required to evaluate the non-analytical contribution were calculated using the DFPT. Therefore, the nac contribution needs to be included in the k calculation of CeO2.

For MD simulations, we used the equilibrium classical MD techniques together with the Green-Kubo linear response formalism[57] as implemented in LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) MD simulation code[58]. The Green-Kubo formalism uses the heat current autocorrelation function (HCACF) (shown in Fig. 10) which decay along a direction as described in ref.[57]. To predict the k of CeO2, the Verlet leapfrog algorithm was implemented[59]. Further, the system was first simulated in a constant number of atoms, pressure, and temperature (NPT) ensemble for 4 ns to ensure it reached equilibrium at the desired temperatures, then the ensemble was switched into a constant number of atoms, volume, and temperature (NVT) ensemble and ran for 4 ns. The heat current autocorrelation function (HCACF) were estimated along with an NVE ensemble calculation which generates an 8 ns raw heat current data at every calculation. Finally, the k value was computed by averaging the k over the time range were the fluctuations were minimal as shown in the inset of Fig. 10b.

Figure 10

(a) The HFACF for CeO2 at a temperature of 300 K with a system size of 10 × 10 × 10 and 20 × 20 × 20 unit cells. (b) The HFACF for CeO2 at a temperature of 300 K and 1500 K with a system size of 10 × 10 × 10 unit cells; inset shows the k fluctuations with time and the shaded region indicates the time range over which the HCACF integral is averaged to predict the k.

Finally, to study the effect of the nanostructuring and calculate the in-plane and cross-plane thermal conductivity of thin films of CeO2 we used the open source almaBTE code[50]. Recently, an exact solution to evaluate the BTE in the cross-plane geometries has been obtained. The film conductivity at thickness L as an integral of phonon frequency is given as,where S is a suppression function that contains the thin film physics and k(ω) denotes the bulk spectral conductivity. More detailed description of the procedure can be found in the following refs[50,60]. It has to be noted that the almaBTE code uses the second order force constant predicted by the FD method. Hence the k values for the nanostructured cases are underpredicted by ~19% as shown for the bulk CeO2.

Materials and experimental details

The 99.9% pure CeO2 powder was obtained from ACROS Organics. The as-received powder was observed under the scanning electron microscope, the micrograph as shown in Fig. 11a revealed that the particles were needle-shaped with dimensions around 20 μm long and 5 μm wide. The as-received powder was sintered using SPS (Thermal Technology LLC 10-3 system) in a graphite die-punch setup as shown in Fig. 11b. The powder die contacts were separated by graphite foils to protect the die from contamination and to reduce friction between the powders and die. A 6 mm hole was drilled from the inner surface of the die, and a pyrometer was focused on this hole to record the temperature of the die. The temperature was controlled by regulating the current passing through the die-punch system, and a Data Acquisition System was used to record displacement, temperature and pressure data.

Figure 11

(a) Scanning electron micrograph of as-received CeO2 powders from ACROS Organic. (b) Schematic diagram of SPS sintering of CeO2.

The pellets were sintered by varying the sintering parameters such as temperature, pressure and hold time. The sintered pellets were then subjected to X-Ray Diffraction (XRD) using BRUKER D8 with Chromium K-alpha radiation to determine the phase and the presence of any residual carbon during the SPS sintering. The Archimedes’ method was used to determine the density of each pellet. Finally, the thermal conductivity was calculated using the laser flash apparatus; Laser flash technique records the thermal diffusivity (α) of the specimens using the Parkers relations[61] given as shown in Eq. (17):where L is the thickness of the specimen and t1/2 is the half of the maximum time taken for the signal to reach the detector. From the measured α, the thermal conductivity as a function of temperature (k (T)) can be measured using the relation (18):where ρ(T) and Cp(T) is the density and heat capacity at constant pressure as a function of temperature respectively. In this work the Cp changes as a function of temperature was determined by comparing the maximum value of the temperature rise with that of a reference pellet, using the relation Cp = Q/(dT.m), where Q represent the energy of the pulsed laser beam, m mass of the specimen, and dT is the maximum value of the temperature rise. The reference pellet used was a certified stainless steel. However, the density changes as a function of temperature has been kept constant. The thermal conductivity measurements were made on cylindrical pellets of diameter 12.7 mm and thickness 2–3 mm.

Atomistic and experimental study on thermal conductivity of bulk and porous cerium dioxide