Influence of mass and charge disorder on the phonon thermal conductivity of entropy stabilized oxides determined by molecular dynamics simulations

High entropy alloys can be defined as having a composition of five or more approximately equimolar principle elements that are randomly arranged on a crystalline lattice.^1–3 These materials were developed primarily within the metal alloy community, where it was recognized that they could be designed such that mixing entropy is a major contributor to stability. Inspired by that work, Rost et al. demonstrated the solid state and vapor synthesis of a new class of high entropy material, entropy stabilized oxides, in which the cation sublattice in a rocksalt structure is randomly populated by five elements.⁴ Following this seminal paper, several other high-entropy ceramics have been reported that have mixed single component/multicomponent sublattices.^5–13 

The original rocksalt composition studied by Rost et al.,⁴ termed J14, was (Mg_0.1Co_0.1Ni_0.1Cu_0.1Zn_0.1)O_0.5. This choice of composition satisfied several criteria thought to be important for both promoting and proving entropy stabilization. These criteria included similar formal charges and ionic radii, limited immiscibility for some of the binary pairs, and some of the lowest energy binary oxide structures not being rocksalt. For example, the wurtzite rather than the rocksalt structure is most stable for ZnO. Subsequent studies demonstrated that other elements can be introduced into the J14 cation sublattice at equimolar concentrations that do not necessarily have radii and formal charges similar to the ions in J14. These elements include Sc, Cr, Sb, Ge, Sn, Li, Ga, and Ca.^5,6,14

Bader charges from Density Functional Theory (DFT) calculations have been used to probe how charge is compensated with the addition of Sc (+3 oxidation state) or Li (+1 oxidation state) into the J14 cation sublattice.¹⁵ It was found that adding Sc reduces a majority of Cu cations and a few Co and Ni cations, while adding Li oxidizes reduces some of the Co cations as well as some Ni and Cu cations. DFT has also been used to validate conclusions based on EXAFS data that a majority of the lattice distortion in J14 are taken up by the oxygen sublattice¹⁶ and that the Cu may show Jahn-Teller distortions involving axial expansion or contraction depending on the relative positions of the Cu cations in the lattice.¹⁷ 

Being an electrical insulator with relatively large mass disorder, J14 should have a relatively small thermal conductivity that is dominated by phonon heat transport. Recent experimental measurements by Braun et al. have confirmed this expectation,¹⁸ but also revealed an unanticipated trend. Traditional heat transport theory suggests that mass scattering should saturate by the addition of the five distinct cations at equimolar concentrations within a high-entropy oxide;¹⁹ yet, experimental measurements showed that the thermal conductivity of J14 is further (and substantially) reduced by the addition of a sixth cation—Sc, Sn, Cr, or Ge—in an equi-molar proportion. Molecular dynamics (MD) simulations using Lennard-Jones potentials suggested that disorder in the inter-atomic forces can potentially explain the further drop in thermal conductivity in these systems,^19–21 a result that is supported by scattering theory and further analysis of the experimental data.¹⁸ Lennard-Jones potentials, however, are not in general a good approximation to forces in oxides, which are dominated by long-range Coulombic forces, and therefore, the physical origin of any inter-atomic force disorder is unclear from these studies.

In this paper, we report the results from a series of molecular dynamics simulations of J14, and J14 + Sc, Sn, Cr, or Ge using the Green-Kubo method^22–24 for calculating phonon-mediated thermal conductivity. The simulations use a Buckingham potential energy function and Coulombic electrostatic forces from atom-centered fixed partial charges. To model the randomness introduced by the charge transfer predicted by DFT within the disordered system, the partial charges were assigned to DFT Bader charges.¹⁵ The remaining parameters in the Buckingham potential were taken unmodified from a literature potential for MgO.²⁵ The simulations show a decrease in thermal conductivity with the introduction of a sixth cation in J14, in agreement with the experiment. Furthermore, by manipulating the masses and charges in the simulations, it is shown that while both mass and charge disorder can lower the thermal conductivity with the addition of a sixth cation, the thermal conductivity lowering can be largely reproduced with disorder in interatomic forces from the charges without an accompanying mass disorder. We also show that when parameterized to the simulations, thermal conductivities from the Bridgman equation, which has been used for ionic liquids where scattering lengths are of the order of inter-atomic distances,²⁶ are significantly lower than those from the MD simulations. This result suggests that the thermal conductivities of these entropy stabilized oxides are still above an amorphous limit for these systems.

Molecular dynamics simulations were carried out using the LAMMPS modeling code²⁷ with thermal conductivity calculated using the Green-Kubo method.^22–24 The interatomic potential was modeled using a pair sum of an exponential-6 function and long-range Coulomb interactions of the form

where N is the total number of atoms, ρ is a parameter, r_ij is the distance between atoms i and j, and A_ij and C_ij are parameters that depend on atoms i and j, and q_i is the charge centered on atom i. The charges were set equal to the Bader charges from DFT calculations on supercells as described below. The other parameters A_ij, ρ, and C_ij for the short ranged repulsive interactions were taken unmodified from a literature potential that was parameterized to model the MgO structure and thermal-mechanical properties, including the experimental room temperature lattice constant, thermal expansion coefficient, and bulk modulus.²⁵ The Bader charge from the DFT calculations on MgO was ±1.7q, which is the same charge reported in Ref. 25 from their fitting procedure. This interatomic potential model allows us to capture not only the randomness of the atomic masses, but also the randomness of the oxygen and cation charges (including charge transfer) as given by DFT.

The DFT calculations were carried out on supercells that were sufficiently large to allow a random population of cations on the face-centered cubic sublattice within the rocksalt structure. For the MD simulations, these supercells were replicated with corresponding periodic boundaries to create larger systems. The DFT supercells contained 240 oxygen anions and 240 cations. There are eight atoms in the rocksalt unit cell, which was replicated 3 × 4 × 5 to make a 480 atom supercell with equimolar metal contents. The arrangement of the latter was generated using the special quasi-random structure algorithm.²⁸ J14 has five different cations so that each DFT unit cell had 48 of each cation type. The other compositions contained 6 different cation types and therefore 40 of each type in the DFT unit cell. The DFT calculations were carried out using plane wave Projector Augmented Wave (PAW) pseudopotential methods as implemented in Vienna Ab initio Simulation Package (VASP).^29–33 The computational parameters are described in detail in Ref. 15. The atom positions were held fixed on a rocksalt lattice, while the total energy was minimized with respect to lattice constant. The final atom positions and calculated Bader charges for each DFT supercell are given in the Appendix.

The MD simulations used 4 × 4 × 4 replicates of the DFT unit cells for a total of 30 720 atoms. For the MD simulations, the potential energy was first minimized with respect to atom positions and lattice constant, followed by equilibration at zero pressure and 300 K for 32 ps using the Nose-Hoover method with a 4 fs time step.^34,35 After equilibration, the equations of motion were integrated at constant energy for 20 ns, over which heat current autocorrelation functions were generated for calculating thermal conductivity. Based on the decay of the autocorrelation functions and the change in calculated thermal conductivity with simulation time, it was determined that this procedure produces converged values of both properties.

Given in Table I and plotted in Fig. 1 are thermal conductivities experimentally measured,¹⁸ calculated from the MD simulations, and given by the Bridgman equation discussed in Sec. IV. Both the experiment and MD simulations show a rough order of magnitude drop in thermal conductivity for J14 compared to MgO, with a further but less dramatic drop with the addition of a sixth cation to J14. For the experiment, the latter is reduced from 3 W/m K to an average of about 1.6 W/m K, a drop of about 50%. The corresponding reduction from the MD simulations is from 4.9 W/m K to an average of 3.2 W/m K or a drop of about 35% from the simulated thermal conductivity of J14.

To better understand the origin of the thermal conductivity drop with the addition of a sixth cation, additional MD simulations were carried out in which atomic mass and charge were manipulated so that their effects on thermal conductivity could be delineated and compared. Three additional cases were studied for each composition; one in which the masses were kept the same but the charges on the oxygen and on the cations were replaced with their respective average charges, one in which the charges from the DFT calculations were kept constant but the atomic masses were replaced with a single mass that gives a reduced mass that matches the reduced mass of the original system, and one in which both the charge and mass are replaced with their average charge and a single mass that produces the reduced mass of the original system. These cases are referred to as homogeneous charge/heterogeneous mass (HMC/HTM), heterogeneous charge/homogeneous mass (HTC/HMM), and homogeneous charge/homogeneous mass (HMC/HMM), respectively. The thermal conductivities from MD simulations for each of these systems, plus the original system, heterogeneous charge/heterogeneous mass (HTC/HTM), are given in Table II and plotted as histograms in Fig. 2. The values of the homogeneous charge and mass used in the simulations are also given in Table II.

For each entropy stabilized oxide composition, the thermal conductivity is highest for the HMC/HMM system, with the highest of these systems being J14 + Sc. Plotted in Fig. 3 are calculated thermal conductivities for the J14 HMC/HMM system, where either the mass or the magnitude of the average cation/anion charge has been increased. From these plots, it is clear that in general thermal conductivity increases with increasing average charge and decreasing average mass. The system HMC/HMM J14 + Sc has both the largest magnitude of average charge and the smallest average mass, and hence the highest thermal conductivity of the HMC/HMM systems. These calculations illustrate that in addition to phonon scattering from mass and charge disorder, the average mass and charge also affect thermal conductivity.

From Table II and Fig. 2, it is clear that heterogeneity in either mass or charge alone lowers the thermal conductivity from the HMC/HMM systems for all compositions, and of the two, charge heterogeneity has the largest influence on reducing thermal conductivity. Furthermore, for each composition, charge heterogeneity without mass heterogeneity reduces thermal conductivity to a value equal (or very close) to that for heterogeneity in both mass and charge. In other words, the effects of phonon scattering from heterogeneity in the mass and in the charge are not additive, and the latter makes a larger contribution to the reduction in thermal conductivity that occurs with the addition of a sixth cation to J14.

Additional insight can be obtained from further analysis of the data in Fig. 2. Adding Sc to J14 and considering mass scattering only (i.e., comparing the HMC/HTM samples for J14 and J14 + Sc) the lighter mass of Sc causes the thermal conductivity to increase with Sc addition to J14. This is despite any enhanced phonon scattering from the added mass disorder. In contrast, considering only scattering from charge disorder (i.e., comparing the HTC/HMM samples for J14 and J14 + Sc), the extra charge scattering from adding Sc with a +3 oxidation state¹⁵ (together with the associated reduction of the other ions) overwhelms the mass effect, and the thermal conductivity of J14 decreases with Sc addition.

Adding Cr to J14 does not significantly change the average mass or charge of J14. Consequently, in Fig. 2, the thermal conductivities of the HMC/HMM systems for J14 and J14 + Cr are very close to one another. In addition, the thermal conductivities of the HMC/HTM systems for the two compositions are the same, implying that any additional phonon scattering from the added mass disorder with the addition of Cr is negligible. In contrast, adding Cr to J14 results in a HTC/HMM and a HTC/HTM model that has a lower thermal conductivity than the corresponding systems for J14. In other words, despite having close to the same average charge, the enhanced phonon scattering from additional charge disorder results in the reduction of thermal conductivity with the addition of Cr to J14.

Adding Sn to J14 slightly increases the average charge and increases the average mass of J14; the latter should act to decrease the thermal conductivity, while the former should increase the thermal conductivity (cf. Fig. 3). Adding Ge to J14 makes a much smaller change in average charge, relative to the other additions, and a somewhat smaller increase in mass. Figure 2 indicates that the thermal conductivity of the HMC/HMM systems for the addition of Sn or Ge is very similar, and both are lower than the equivalent systems for the addition of Sc or Cr. Like the other systems, disorder in charge is more effective for scattering phonons than mass disorder based on the thermal conductivities of the HTC versus HTM systems.

The results and analysis presented here for adding a sixth cation to J14 provide new insight into a prior experimental study that measured the thermal conductivity of entropy-stabilized oxides and confirm the conclusion from prior simulations that used Lennard-Jones pair potentials that phonon scattering from variations in interatomic forces are needed to explain the drop in thermal conductivity experimentally measured for the addition of a sixth cation to J14.²⁰ At the same time, the present work has provided some new insights; for example, without phonon scattering from the interatomic forces, addition of Sc to J14 would likely have increased thermal conductivity because of its lighter mass.

To gain additional insight into the results of the MD simulations reported here, the thermal conductivity for each composition was also calculated using the Bridgman equation. This equation has been used for ionic liquids where scattering lengths are of the order of inter-atomic distances;²⁶ hence, it represents an analytic limit of an amorphous structure for ionic materials. This equation gives the thermal conductivity $k$ as a product of the density ρ, heat capacity c_v, unidirectional mean molecular speed v_y, and lattice spacing a

For our calculations, v_y is taken as the sound speed, which is approximated as the longitudinal sound speed given by the square root of the bulk modulus divided by the density, and the interatomic spacing is taken as (V/N)^1/3. The bulk moduli were calculated using the second derivative of a polynomial fit to the potential energy as a function of volume for each system after relaxation of the atom positions and periodic boundaries with respect to potential energy. The heat capacities were calculated from the variance in the potential energy during the MD simulations as described above. Like the MD simulations, because the interatomic potentials have not been specifically fit to the J14 and related compositions, the calculations are not intended to be quantitative. Instead, by parameterizing the Bridgman equations to the MD systems, these calculations are used to quantify how close our model approximates an amorphous limit.

Given in Table III are the calculated properties for each model system that are used in the Bridgman equation, the experiment values.^36–38 The calculated thermal conductivities are given in Table I and Fig. 1. The thermal conductivities from the Bridgman equation are similar for all compositions and significantly lower than those given by the MD simulations. Hence, we conclude that while the addition of a sixth cation to J14 lowers the thermal conductivity largely due to enhanced scattering from charge disorder, the systems are still above this amorphous scattering limit.

Classical molecular dynamics simulations have been used to characterize the influence of mass and charge disorder on the phonon-mediated thermal conductivity of five entropy stabilized oxides, J14 with composition (Mg_0.1Co_0.1Ni_0.1Cu_0.1Zn_0.1)O_0.5, and J14 plus Sc, Sn, Cr, or Ge in equi-molar cation proportions. By manipulating the mass and the atom-centered charges used in the Coulombic contribution to the potential energy, it was shown that the phonon scattering from disorder in atomic charges that enter the interatomic potential can explain prior experimental results that show a lowering of thermal conductivity with the addition of a sixth cation to J14. A comparison of thermal conductivities from the MD simulations to predictions from the Bridgman equation using properties from the simulations suggests that despite the added phonon scattering from disorder in mass and atomic forces, the thermal conductivity is still above an analytic limit for an amorphous system.

See supplementary material for the final atom positions and calculated Bader charges for each DFT supercell.

This work was supported by the U.S. Office of Naval Research MURI program (Grant No. N00014-15-1-2863) and by the National Science Foundation (NSF) Ceramics through Award No. 1610844.