Reduced dependence of thermal conductivity on temperature and pressure of multi-atom component crystalline solid solutions

In search of new functional materials, high entropy materials have attracted much attention due to their unique mechanical,^1–3 electrical, and magnetic properties.⁴ These materials are generally formed by mixtures of at least 5 different elements and demonstrate enhanced configurational entropy, particularly at high temperatures.⁵ Along with thermodynamic stability, these materials have shown improved strength and toughness (including oxidation resistance at elevated temperatures) due to their microstructure, which makes them attractive candidates for high temperature, pressure, and harsh environment applications.⁶ An important factor for consideration in these applications is the management of heat and thermal transport properties of multicomponent materials. In this regard, only a few studies have focused on studying the lattice thermal conductivity of ternary and quaternary alloys,^7–9 and the effect of high pressure and temperature on the lattice thermal conductivity of multi-atom component solid solutions has largely been unexplored.

Recently, Rost et al.¹⁰ have demonstrated the ability to fabricate single phase, crystalline oxide solid solutions with up to five or six different elements accommodating the lattice. The crystallographic configuration of their multicomponent solid solutions, which they refer to as “entropy stabilized oxides,” is comprised of an anion sublattice with multiple components filling up the cation species in a rocksalt configuration. With respect to these multicomponent solid solutions, we have previously shown that the phonon thermal transport in non-metallic crystalline alloys with five (or more) different elemental species offers the ability to push the limits of traditionally assumed phonon scattering theories, such as the virtual crystal approximation (VCA) approaches.¹¹ However, critical questions regarding the thermal properties of multicomponent solid solutions still remain unanswered, such as: (i) how does the introduction of mass impurities alter the phonon mode properties and affect the spectral contributions to thermal conductivity, (ii) how does chemical order-disorder transition in these alloys dictate the various scattering mechanisms and thus the thermal conductivities, and lastly, (iii) how do high pressures and temperatures affect the thermal properties of disordered multicomponent alloys in comparison to homogeneous crystalline systems?

Thus, in this study we gauge the effect of mass disorder on the spectral thermal conductivity and its dependence on temperature and pressure of multicomponent solid solutions. We accomplish this by conducting molecular dynamics simulations to predict the thermal conductivities of Lennard-Jones (LJ)-based solid solutions (with one to five component atoms in the crystalline system) at a range of uniaxial strain levels and simulations conducted at different temperatures. We show that for multicomponent alloys, increasing only the mass impurity beyond a 2-component alloy by adding atoms with different masses in the solid solution does not lead to a significant change in the spectral contributions. However, increasing the impurity concentration or changing the local force-field of the impurity atoms in the solid solution has a relatively significant impact on the spectral contributions to thermal conductivity. Furthermore, the thermal conductivity of ordered alloys (comprising a host sublattice) is shown to demonstrate a pronounced temperature dependence in comparison to their disordered counterparts, which is ascribed to the different scattering mechanisms responsible in dictating thermal conductivities in these alloys. In comparison to a homogeneous solid, solid solutions (especially in the chemically disordered state) show reduced pressure dependence as well, which becomes more prominent as the number of components is increased. This is attributed to the fact that while anharmonic effects in homogeneous solids lead to the large temperature and pressure dependence in thermal conductivity, impurity scattering in solid solutions leads to a largely reduced dependence on pressure and temperature.

We employ non-equilibrium molecular dynamics (NEMD) simulations and the recently formulated spectral analysis technique to study the effect of mass impurities on the spectral thermal conductivity for a range of pressure and temperature of multicomponent solid solutions. As we are interested in studying the general effects as opposed to material specific properties, we employ the widely used Lennard-Jones (LJ) potential to describe the interatomic interactions. The 6–12 LJ potential is given as $U(r)=4ε[(σ/r)12−(σ/r)6]$⁠, where U is the interatomic potential, r is the interatomic separation, and σ and ε are the LJ-length and -energy parameters, respectively. The cutoff distance is set to $2.5σ$ for all the simulations and the time step is set to 1 fs throughout the simulations. All the MD simulations are performed using LAMMPS.¹² 

To begin, the length and energy parameters are modeled for argon with $σ=3.405$ Å and $ε=10.3$ meV, respectively. The lattice constant is set to $a0=1.56σ$ and the atoms are arranged in a fcc lattice. For all simulations, the size of the computational domains is $10a0×10a0×80a0$ with periodic boundary conditions applied in all directions. When conducting NEMD simulations, periodic boundary conditions are only applied in the x- and y-directions, whereas fixed boundaries with 2 monolayers of atoms at each end are placed in the z-direction. To accommodate mass impurities in a homogeneous crystal, the impurity masses are equally (and randomly) distributed throughout the host crystalline lattice at the prescribed alloy fraction. For our simulations, the mass of the host lattice is set to $mH=20$ g mol⁻¹ and the impurity atoms are introduced at 20 g mol⁻¹ mass increments (where impurity atom masses range from 40 to 100 g mol⁻¹ and are labeled $mA=40$ g mol⁻¹, $mB=60$ g mol⁻¹, $mC=80$ g mol⁻¹, and $mD=100$ g mol⁻¹) from 2 to 5 different atomic masses comprising the impurity atoms distributed in equal amounts. The impurity alloy concentration, x, and the parameters used for the various domains are listed in Table I, where “$*$” represents all atomic species. For the 5-component solid solution with both mass and bond defect, we increase σ and decrease ε by 15% for the impurity atom interactions in comparison to the values for the atoms in the host lattice; perturbations of σ and ε by 15% for the impurity atoms were chosen arbitrarily as shown in Table I; however, we note that in alloys such as FePt, the energy parameter can vary as much as $∼15%$ and the length parameter can vary by as much as $∼12%$⁠; note, even in these metallic FePt alloys, the phonon contribution to thermal conductivity can be non-negligible.³⁸ The interaction parameters between the host and the impurity atoms (⁠$εhi$ and $σhi$ for the energy and length parameters, respectively) are determined by the traditional mixing rule expressions, $εhi=εhhεii$ and $σhi=(σhh+σii)/2$⁠, for the corresponding energy and length parameters, respectively.¹³ 

The computational domains are equilibrated under the Nose-Hoover thermostat and barostat;¹⁴ the number of atoms, volume, and temperature of the simulation is held constant followed by a isothermal-isobaric ensemble with the number of particles, pressure, and temperature of the system held constant for a total of 2 ns at 0 bar pressure. For the NEMD simulations, a fixed amount of energy is added per time step to a warm bath at one end and the equal amount of energy is removed from a cool bath at the other end. This is performed with the length of the baths set to $10a0$ in the z-direction, and the dynamics are carried out under a microcanonical ensemble, or the NVE integration, with the number of particles (N), volume (V), and energy (E) held constant. After $∼2$ ns, a steady-state temperature gradient in the z-direction is established by averaging the temperature for atoms in each monolayer for a total of another 5 ns. Examples of the steady-state temperature profiles for a homogeneous domain with m = 20 g mol⁻¹ and a 5-component solid solution are shown in Fig. 1(a). These temperature profiles are then used to predict the thermal conductivity of the various structures by invoking the Fourier law, $Q=−κ∂T/∂z$⁠, where the applied flux is in the z-direction. To check for convergence of system size on the thermal conductivity predictions, we simulate computational domains with varying lengths in the direction of the applied heat flux. As shown in Fig. 1(b) for 1-component and 5-component structures, a domain size of 80 unit cells (⁠$∼43$ nm) is large enough to predict size-independent thermal conductivities, which suggests that our choice of $10a0×10a0×80a0$ for the system sizes is computationally reasonable and avoids artifacts due to small system sizes. Along with the computational domain size, the applied heat flux does not alter the calculated thermal conductivity as well.

For the spectral contributions to thermal conductivity, the heat flux is spectrally resolved by the relation,¹⁵ $Q=∫0∞dω2πq(ω)$⁠, where ω is the angular frequency and $q(ω)$ is the spectral heat current. The details of the procedure are outlined in our previous work in Ref. 16. In general, this heat current between an atom i and j is proportional to the correlation between the interatomic force $F→ij$ between the atoms and the velocities, $qi→j(ω)∝〈F→ij·(v→i+v→j)〉$⁠, where the brackets denote steady-state non-equilibrium ensemble average.^17–19 We tabulate the forces and velocities for the atoms under consideration for a total of 10 ns with 10 fs time intervals under the NVE integration to calculate the spectrally resolved heat flux, $q(ω)$⁠.¹⁶ Note, as the dot product of the forces and velocities includes separate contributions in the x-, y-, and z-directions, the x- and y-components describe the in-plane or transverse mode contributions, and the component in the z-direction describes the out-of-plane or longitudinal mode contribution to the total heat current. It should also be noted that the assumption that all the modes experience the same temperature drop holds for our spectrally decomposed thermal accumulation calculations. The reader is referred to Ref. 20 for a detailed NEMD-based approach that calculates the spectral phonon temperatures and spatial temperatures of all the individual phonon modes in systems using more material specific potentials.

Figure 2(a) shows the NEMD-predicted thermal conductivities for the 2- to 5-component solid solutions for a range of impurity mass concentrations. It is clear from the figure that the addition of impurity masses ultimately reduces the thermal conductivity. However, in accordance with our previous work, the thermal conductivity of solid solutions beyond 4 impurity atoms, within uncertainties, is similar for the entire range of alloy concentrations. It should be noted that the generality of these results holds for cases with different host masses, different percentages of impurity masses introduced in the host lattice, and different mass increments of the impurity atoms, whereby adding more impurity atoms beyond a critical limit does not significantly change the thermal conductivity.¹¹ This can be attributed to the fact that the term proportional to mass disorder [Γ in $τ−1=δ3Γω4/(4πv3)$ where $δ3$ is the atomic volume, v is the velocity, k is the wavevector, and ω is the frequency] dictating thermal conductivity of these materials approaches a constant value even with the increase in the number of components beyond a certain limit that differ in mass alone.¹¹ To reduce the thermal conductivity of these alloys even further, the local strain-field has to be altered. In our case, this can be achieved via changing the energy or length parameter in the LJ-potential as mentioned above. As shown in Fig. 2(a), the thermal conductivity of the 5-component alloys can be lowered significantly (even beyond the mass scattering limit) by perturbing the local strain-field. It is also interesting to note that the minimum in the thermal conductivity shifts to higher concentrations as the number of components is increased beyond the binary alloy. This is consistent with our recent experimental findings that suggests that a larger difference in mass and crystallographic properties of the parent materials results in a more asymmetric thermal conductivity trend with concentration.²¹ 

To study the effect of mass impurities on the vibrational properties of the LJ-based alloys, we calculate the density of states (DOS) of the solid solutions by performing the Fourier transform of the velocity autocorrelation function.^22,23 Figure 2(b) shows the vibrational DOS for a 2-component alloy (with host lattice set to m = 20 g mol⁻¹ and the impurity atoms set to m = 40 g mol⁻¹) at different alloy concentrations. The shaded area represents the DOS for the host lattice only. The increase in alloy concentration results in a monotonous shift of the DOS to lower frequencies due to the increase in the concentration of the heavier impurity atoms that have a lower cutoff frequency compared to the host lattice. Similarly, in Fig. 2(c), we plot the DOS for the computational domains with 50% alloy fraction and increasing number of components from 2 to 5-components for the impurity masses. The increase in the number of components by incrementing the mass of the impurity atoms has a similar effect on the DOS as the spectrum is shifted to lower frequencies with the depletion of high frequency vibrations. However, it should be noted that for the case of the structure with the 5-component solid solution with change in the local force-field, a greater depletion in the high frequency phonons is observed [see Fig. 2(c)].

In order to study the effect of impurity mass on the spectral contributions to thermal conductivity, we calculate the heat current accumulation for these structures as discussed above. Figure 3(a) shows the accumulation for the host lattice and 2-component solid solutions with varying impurity concentrations. For all structures, the largest contribution to thermal conductivity comes from phonons in the mid-frequency range of the DOS. This is intuitive as these frequencies possess the largest population of phonons in the structures as shown in Fig. 2(b). Similar to the host lattice, the 2-component alloys with various impurity concentrations show $∼60%$ contribution from transverse modes and $∼40%$ contribution from longitudinal modes. Increasing the concentration of impurities is shown to slightly shift the accumulation to lower frequencies, which is in line with the shift in the DOS [Fig. 2(b)].

We also calculate the heat current accumulation for alloys with increasing number of components for impurity concentration of 50% as shown in Fig. 3(b). In contrast to the slight shift in accumulation with increasing impurity concentration, increasing the number of components has negligible effect on the spectral contributions to thermal conductivity even though the DOS shifts to lower frequencies [going from 2 to 5-component alloys as shown in Fig. 2(c)]. This suggests that the mean free paths of phonons in these LJ-structures do not significantly change with the addition of more impurity components (and therefore the individual contributions from the normal modes do not change significantly for these alloys) with the increment of impurity mass scattering alone at a fixed impurity concentration. However, the spectral contributions to thermal conductivity can be altered by changing the local force-field of the impurity atoms by changing the energy and length parameters in the LJ potential. For the 5-component alloy with the 15% reduction in ε and 15% increase in σ, the spectral contributions to thermal conductivity shift significantly to lower frequencies. These results suggest that mass impurity alone does not alter the spectral contributions significantly. However, a change in the local strain-field can cause a drastic alteration to the spectral contributions.

As seen from the above results, the overall disorder strength in solid solutions is determined by alloy concentration, mass ratio, and stiffness ratio. In terms of analytical models, the virtual crystal approximation (VCA), which treats the phonon properties of solid solutions as an average property determined from the individual components, is a widely used model to explain some of the experimentally determined thermal conductivities of alloys.^24–27 In this context, Larkin and Mcgaughey,²⁸ through detailed lattice dynamics and molecular dynamics calculations, have shown that the VCA underpredicts the lifetimes of high-frequency modes for similar LJ-based binary alloys. Similarly, our spectral heat flux calculations for the LJ-based alloys suggest that the high-frequency modes contribute a significant amount to the total heat current. Therefore, the use of the VCA for these structures can lead to erroneous results, as pointed out by Larkin and Mcgaughey.²⁸ However, for more realistic alloys such as SiGe (which are stiffer in nature as compared to the LJ-based alloys studied in this work), where the dominant heat carriers are low-frequency vibrations,²⁹ the VCA has been shown to be more appropriate to describe thermal transport since the underprediction of lifetimes of high frequency phonons ultimately does not have significant impact on the predicted thermal conductivity.²⁸ We note that our LJ-based structures are fundamentally different than commonly studied alloys with stiffer bonds and less-anharmonic potentials in which low frequency vibrations could be the dominant heat carriers. Along these lines, anharmonic interactions are more pronounced in the LJ-potential, which can give rise to very different temperature trends in thermal conductivity compared to more realistic systems such as Si defined by the Stillinger-Weber potential, which have been shown to be relatively more harmonic in nature;³⁰ low frequency phonons in the latter contribute much more significantly to thermal conductivity as compared to the former (even at high temperatures). Therefore, we caution that our results presented here should not be directly and quantitatively compared to specific material systems with stronger bonding environments.

Next, we consider the effect of chemical order and disorder on the temperature dependent thermal conductivities of our LJ-based multicomponent solid solutions. To this extent, we construct computational domains with L1₀-type configurations. These configurations form a simple 1 × 1 superlattice in the [001] direction. For these structures, we construct multicomponent alloys by including the impurity atoms to the monolayers next to the host lattice, thus ensuring the periodicity and order in the host lattice. The computational domain with a binary L1₀ and a 5-component L1₀ configurations is shown in Fig. 4(a). For comparison, we have also included the schematic of the computational domain for our homogeneous and disordered 5-component crystalline solid solutions. Similar to the crystallographic configuration of the entropy stabilized oxides experimentally realized in Ref. 10 where the anion sublattice is ordered while the cation sites can accommodate multiple components, our structures also include periodicity in the host sublattice (with every other atom is periodically repeated in the [001] direction). The side view of the 5-component L1₀ structure is shown in Fig. 4(b) for clarity. It should be noted that one major difference in the structures between the entropy stabilized oxides and our LJ-based solid solutions is that the former is in the rocksalt configuration whereas the latter are fcc crystal structures.

Figure 4(c) shows the NEMD-predicted thermal conductivities as a function of normalized temperature (with respect to the melting temperature of the LJ-based solid) for the L1₀ solids with 2- and 5-components. For comparison, we have included the NEMD results for a homogeneous LJ-argon (red squares) crystalline domain, and 2-component (blue circles) and 5-component (green diamonds) alloys in the disordered crystalline configurations. At low temperatures, the thermal conductivity of the 2-component L1₀ alloy demonstrates higher thermal conductivity as compared to its disordered counterpart (2-component alloy). However, at higher temperatures, the thermal conductivity of the ordered 2-component alloy surpasses that of its disordered counterpart. Similarly, the thermal conductivity of the homogeneous LJ-argon at high temperatures can approach the thermal conductivity of the 2-component alloys. Since the ordered and disordered states have the same atomic compositions and lattice constant, they possess the same heat capacity and molar density. Therefore, the difference in their thermal conductivities mainly arises due to the different scattering mechanisms and group velocities of the phonons.³¹ 

With regard to phonon scattering, the total scattering time is usually approximated by the Matthiessen's rule, $τtot−1=τU−1+τimp−1$⁠, where $τU$ is the scattering rate due to (three-phonon) Umklapp scattering and $τimp$ is the scattering due to impurities in the crystal. For disordered alloys, $τimp$ dictates thermal conductivity across much of the temperature range while for the homogeneous crystal $τU$ dictates the thermal conductivity.^31–33 Therefore, due to the enhanced role of $τU$⁠, and the different dominant scattering mechanisms and their dependence on temperature, the thermal conductivity of the ordered alloys (or even the homogeneous case) can be lower than the disordered state. This can also be partially attributed to the reduction in the Brillouin zone size in the ordered alloys (with relatively larger unit cells), which can increase phonon-phonon scattering rates, as discussed by Duda et al.³¹ This is further demonstrated by the 5-component alloys, where at low temperatures the thermal conductivity of the ordered (L1₀ configuration) and the disordered alloys is similar but as the temperature is increased, the thermal conductivity of the ordered state is lower than that of the 5-component disordered alloy. Furthermore, at relatively high temperatures, the thermal conductivity of the ordered state approaches the thermal conductivity of the alloy of impurities in the amorphous state, which is generally considered as the theoretical minimum in thermal conductivity.^34,35 We also include the thermal conductivity of amorphous LJ-argon for comparison. Note, the dashed lines for the crystalline domains represent power law fits (⁠$κ=aT−1+b$⁠, where a and b are variables), which demonstrate the extent of Umklapp dominated thermal conductivity. These observations suggest that by creating an ordered sublattice (as experimentally realized in Rost et al.'s work¹⁰) in a solid solution with multiple components, the thermal conductivity can be lowered as much as the theoretical limit at relatively higher temperatures.

Next, we study the influence of strain on the thermal conductivities of multicomponent alloys, since strain has an apparent effect on the spectral thermal conductivity (as suggested by our earlier discussions). For this purpose, we apply uniaxial strain in the z-direction and perform NEMD simulations at different temperatures. We first consider the effects of strain and pressure on the thermal conductivity of a homogeneous crystal (LJ-argon). Figure 5(a) shows the thermal conductivity of LJ-argon as a function of temperature and uniaxial strain. Consistent with prior results on LJ-argon,^30,36 the thermal conductivity monotonically increases with compression and decreases with tension for the temperature range shown in Fig. 5(a). Also, the thermal conductivity monotonically decreases with temperature for all strain levels. However, it is interesting to note that the temperature dependence is greater at compressive strains due to increasing anharmonicity of the crystal lattice.

Figure 5(b) shows the thermal conductivity dependence on pressure for 2-, 3-, and 5-component alloys (at impurity concentrations of 50%) calculated at 17% (square symbols) and 40% (triangle symbols) of the melting temperature of the solid solution; the melting temperature of LJ-argon has been calculated to be $∼87$ K.³⁷ For these multicomponent alloys, the temperature dependence of thermal conductivity (for the range of strain levels studied in this work) decreases as the number of components increases. This is clear from the decreasing separation (with increasing number of components) of the thermal conductivities between the two temperatures shown in Fig. 5(b). Moreover, unlike for the case of LJ-argon, the weak temperature dependence for the 5-component alloy at compressive and tensile strains suggests that scattering due to three-phonon processes and anharmonicity induced Umklapp scattering is not the dominant mechanism dictating thermal conductivity for these structures. It should also be noted that the application of compressive or tensile strains leads to a drastic shift of the spectral contributions to higher and lower frequencies, respectively, for the homogeneous crystal as well as the solid solutions.

Furthermore, to gauge the relative pressure dependence of thermal conductivity for multicomponent alloys, we normalize the thermal conductivity of our LJ-argon, 2-component and 5-component solid solutions by their respective thermal conductivities predicted at ambient pressure [as shown in Fig. 5(c)]. The increment in the number of components reduces the dependence of thermal conductivity on strain. For example, at compressive and tensile strains of 6% for LJ-argon, the thermal conductivity can be increased by $∼85%$ and reduced by $∼63%$⁠, respectively. However, this dependence is reduced for the alloys since for the 5-component alloy, at compressive and tensile strains of 6%, the thermal conductivity can be increased by $∼50%$ and reduced by $∼37%$⁠, respectively.

For the homogeneous crystals, due to increased anharmonicity in the crystal lattice at increased compressive strains, Umklapp processes dictate thermal conductivity, whereas, for the disordered solid solutions, increase in anharmonicity is overshadowed by the effect of impurity mass scattering that leads to a reduced pressure dependence of thermal conductivity for the disordered alloys. It should be noted that the L1₀-type alloys also showed a pronounced dependence on pressure, similar to the homogeneous crystals due to Umpklapp dominated scattering mechanisms. Therefore, from the observations and discussions outlined in the preceding paragraphs, it can be confirmed that the increment in mass scattering by increasing the number of components in a disordered crystalline solid solution can lead to a reduced dependence of thermal conductivity on temperature and pressure.

We have investigated the role of mass impurities on the spectral thermal conductivity (and its dependence on pressure and temperature) for LJ-based multicomponent solid solutions via molecular dynamics simulations. Our results suggest that for multicomponent alloys, increasing only the mass impurity-based disorder by adding atoms with different masses in the solid solution does not lead to a significant change in the spectral contributions to thermal conductivity. However, increasing the impurity concentration or changing the local force-field of the impurity atoms in the solid solution has a relatively significant impact on the spectral contributions to thermal conductivity. Furthermore, the thermal conductivity of solid solutions comprised of an ordered sublattice for the host atoms (and impurity atoms randomly distributed) is shown to demonstrate a pronounced temperature dependence in comparison to their fully disordered counterparts. We ascribe this to the different scattering mechanisms responsible for dictating thermal conductivities in ordered vs. disordered solid solutions. Moreover, in comparison to a homogeneous solid, solid solutions (especially in the disordered state) show a reduced pressure dependence as well, which becomes more prominent as the number of components is increased. This is attributed to the fact that while anharmonic effects in homogeneous solids lead to the large temperature and pressure dependencies of thermal conductivity, impurity scattering in solid solutions leads to a largely reduced dependence on pressure and temperature.

We acknowledge the financial support from the Office of Naval Research MURI program (Grant No. N00014-15-1-2863) and the National Science Foundation (Grant No. CBET-1706388).