Thermal transport in bulk ternary III-V arsenide (III-As) semiconductor alloys was investigated using equilibrium molecular dynamics with optimized Albe-Tersoff empirical interatomic potentials. Existing potentials for binary AlAs, GaAs, and InAs were optimized to match experimentally obtained acoustic-phonon dispersions and temperature-dependent thermal conductivity. Calculations of thermal transport in ternary III-Vs commonly employ the virtual-crystal approximation (VCA), where the structure is assumed to be a random alloy and all group-III atoms (cations) are treated as if they have an effective weighted-average mass. Here, we showed that it is critical to treat atomic masses explicitly and that the thermal conductivity obtained with explicit atomic masses differs considerably from the value obtained with the average VCA cation mass. The larger the difference between the cation masses, the poorer the VCA prediction for thermal conductivity. The random-alloy assumption in the VCA is also challenged because X-ray diffraction and transmission electron microscopy show order in InGaAs, InAlAs, and GaAlAs epilayers. We calculated thermal conductivity for three common types of order (CuPt-B, CuAu-I, and triple-period-A) and showed that the experimental results for In_{0.53}Ga_{0.47}As and In_{0.52}Al_{0.48}As, which are lattice matched to the InP substrate, can be reproduced in molecular dynamics simulation with 2% and 8% of random disorder, respectively. Based on our results, thermal transport in ternary III-As alloys appears to be governed by the competition between mass-difference scattering, which is much more pronounced than the VCA suggests, and the long-range order that these alloys support.

## I. INTRODUCTION

III-V arsenide (III-As) ternary alloys and their superlattices (SLs) are widely used in optoelectronic and thermoelectric devices.^{1,2} In particular, the use of III-As alloys brings about great flexibility in the design of the active region of quantum cascade lasers.^{3,4} Despite the wide popularity, even the bulk thermal properties of these ternary alloys have not been completely characterized to date,^{5} which hinders the analysis of device thermal performance.

Measuring the thermal conductivity of bulk ternary alloys requires growing them strain-free on a substrate. Since the GaAs substrate is lattice matched to Al_{x}Ga_{1–x}As with any Al composition, the thermal conductivity measurement of AlGaAs was carried out early on, though only at room temperature.^{6} Measurements on In_{x}Ga_{1–x}As and In_{x}Al_{1–x}As are more challenging, as they each only have a specific composition (In_{0.53}Ga_{0.47}As and In_{0.52}Al_{0.48}As) that is lattice matched to InP, a common substrate for the growth of mid-infrared quantum cascade lasers.^{7} As a result, In_{x}Ga_{1–x}As and In_{x}Al_{1–x}As often exist in pairs inside a strain-balanced SL on top of InP. As the influence that the interfaces in SLs have on thermal transport is poorly understood,^{8,9} it is very difficult to extract the thermal conductivity of each material from the measurements on SLs. Abrahams *et al.* measured the room-temperature thermal conductivity of In_{x}Ga_{1–x}As with various In compositions, where the samples were formed through zone leveling.^{10} However, this method is not applicable to In_{x}Al_{1–x}As, as they oxidize easily. Koh *et al.* measured the room-temperature thermal conductivity of In_{0.53}Ga_{0.47}As and In_{0.52}Al_{0.48}As to be 4.02 W/m K and 2.68 W/m K, respectively, using both the 3*ω* technique and time-domain thermoreflectance (TDTR).^{11} Sood *et al.* used TDTR to measure the cross-plane thermal conductivity of In_{0.53}Ga_{0.47}As/In_{0.52}Al_{0.48}As SLs with different layer thicknesses and extracted the thermal conductivity of In_{0.53}Ga_{0.47}As and In_{0.52}Al_{0.48}As to be 5.0 W/m K and 1.2 W/m K, respectively. The two sets of results disagree, with the disparity being particularly pronounced for In_{0.52}Al_{0.48}As. Furthermore, both sets of results were obtained for In_{0.53}Ga_{0.47}As and In_{0.52}Al_{0.48}As lattice matched to InP. We still have very limited knowledge about the thermal conductivity of In_{x}Ga_{1–x}As and In_{x}Al_{1–x}As bulk materials with various compositions.

Our recent theoretical work solved the Boltzmann transport equation for phonons under the relaxation-time approximation in order to describe the thermal conductivity of III-As semiconductors.^{12} We accurately captured the thermal conductivity of binary compounds (AlAs, GaAs, and InAs) from 100 K to 600 K with full acoustic phonon dispersions. Additionally, we adopted the virtual crystal approximation^{13} (VCA) to compute the thermal conductivity of ternary alloys Al_{x}Ga_{1–x}As, In_{x}Ga_{1–x}As, and In_{x}Al_{1–x}As at room temperature. For In_{x}Ga_{1–x}As and Al_{x}Ga_{1–x}As, our calculated thermal conductivity agreed reasonably well with experimental values with various *x* compositions.^{6,10,14} For In_{x}Al_{1–x}As, an experimental thermal conductivity value was only available for In_{0.52}Al_{0.48}As.^{11,15} Our calculated thermal conductivity is 2.82 W/m K, which is close to the value of 2.68 W/m K obtained by Koh *et al.*^{11} but far from the value of 1.2 W/m K obtained by Sood *et al.*^{15} Considering that we have benchmarked the technique for binaries, we speculate that the VCA is not entirely suitable for characterizing the thermal conductivity of ternary III-As alloys. (We do note that the VCA has been employed in the analysis of various alloys in recent years.^{16–19})

Within the VCA, there are two major assumptions whose validity needs to be questioned when it comes to the calculation of thermal conductivity in ternary alloys A_{x}B_{1–x}C. (1) All cations can effectively be replaced with an effective averaged cation, whose mass is calculated as the weighted average of cation masses, i.e., m_{III} = *x*m_{A} + (1 – *x*)m_{B}. (2) The alloy is random, i.e., cation sites are randomly taken by atom A or atom B, with the frequency proportional to cation abundance [*x* and (1 – *x*), respectively]. Although an effective mass-difference scattering rate is often employed in Boltzmann-equation-based approaches to compensate for the fact that assumption (1) eliminates scattering caused by the cation mass difference,^{13} these rates are calculated under the assumption that the cation mass difference is small with respect to the average cation mass. Therefore, the validity of this perturbative mass-difference-scattering approach becomes suspect when the cation masses are very different. In our case, m_{Al} = 26.98 au, m_{Ga} = 69.72 au, and m_{In} = 114.82 au, and so, it makes sense that the model would work better for In_{x}Ga_{1–x}As and Al_{x}Ga_{1–x}As than it does for In_{x}Al_{1–x}As since the mass difference between Al atoms and In atoms is significant. The second assumption is directly contradicted by a number of X-ray and transmission-electron-microscopy (TEM) experiments conducted on ternary III-V alloys grown epitaxially on a GaAs or an InP substrate.^{20–31} These experiments show that group III atoms in ternary III-As epilayers are arranged with a certain order rather than completely randomly. This ordered structure is also backed up by an observable change in the width of the electronic bandgap of these materials. The most commonly observed types of ordering for In_{0.53}Ga_{0.47}As and In_{0.52}Al_{0.48}As are the CuPt-B order^{20–22,24} and the triple-period-A (TPA) order.^{23,25–29} The CuAu-I order is more common in Al_{0.5}Ga_{0.5}As,^{30,31} but it could also be observed in In_{0.53}Ga_{0.47}As or In_{0.52}Al_{0.48}As under appropriate growth conditions.^{26,31}

In this work, we used molecular dynamics (MD) to study thermal transport in ternary III-As alloys. In a MD simulation, the mass and location of each atom are tracked in real space; therefore, it is straightforward to include individual atom masses and the exact alloy structure explicitly in the simulation. We adopted the Tersoff^{32,33} empirical interatomic potentials (EIPs) to describe the interaction between atoms. The parameters we used in the EIPs were obtained by an optimization technique, starting from the existing values,^{34–36} with the goal of better capturing the phonon properties of binary III-As semiconductors.^{37} We then used the optimized EIPs to study the influence of the mass difference and ordering on the thermal conductivity of ternary III-As alloys. The results show that the atom mass must be explicitly considered, in particular in In-based compounds. Ordering can significantly increase the thermal conductivity, and partial ordering is a likely reason for the higher experimentally obtained thermal conductivities than what the VCA within the Boltzmann transport framework predicts. Moreover, thermal conductivity can vary a great deal with the introduction of even low levels of disorder.

In Sec. II, we discuss the way we compute the bulk thermal conductivity from an MD simulation and our process of optimizing the Tersoff potentials for AlAs, GaAs, and InAs. In Sec. III, we apply the optimized potentials to ternary alloys. We isolate the effect of the mass difference by comparing the thermal conductivity obtained from structures where each cation has its own mass and those where all cations have the effective VCA mass. In Sec. IV, we study the thermal conductivity of bulk ternary III-As alloys with perfect ordering and lattice-matched In_{0.53}Ga_{0.47}As and In_{0.52}Al_{0.48}As that contain small levels of random disorder. Finally, we summarize our findings in Sec. V.

## II. OPTIMIZATION OF TERSOFF POTENTIALS

All the MD simulations in this work were carried out in the LAMMPS^{38} package. We used equilibrium MD together with the Green-Kubo formula to compute the thermal conductivity of compound semiconductors. (Details are fairly standard but have been given in Appendix A, for completeness.)

Most Tersoff potentials (see Appendix B) are parameterized so as to accurately capture the mechanical properties of materials.^{32–36,39–41} However, for the EIPs to be good at describing thermal transport, they also must produce good phonon properties. Lindsay and Broido^{37} optimized the Tersoff potential of carbon for thermal transport. They used a $\chi 2$-minimization, where

where *i* runs over all the physical properties they optimize for; *η _{i}* and $\eta \u2009exp\u2009,i$ are the calculated and experimental values of the physical property

*i*, respectively; and

*ζ*is the weighting factor determining the relative importance of each physical property in the optimization process. They assigned the most weight to the phonon frequencies and sound velocities because of their crucial role in thermal transport.

_{i}^{37}

Here, we followed a similar minimization approach to improve Tersoff EIPs for describing thermal transport in III-Vs. The physical properties we optimized for are the lattice constant *a*_{0}, the cohesive energy *E*_{coh}, and sound velocities *v _{s}*

_{,TA}and

*v*

_{s}_{,LA}. Note that the accuracy of phonon dispersions only guarantees the accuracy of the second derivatives of the EIPs (or the IFCs) ( Appendix C) while the finite thermal conductivity of crystals originates from the phonon–phonon scattering, related to the third and higher derivatives of the EIPs. Since the strength of these higher-order interactions is not easily accessible in experiments, we used the temperature-dependent thermal conductivity measured in experiments as an additional gauge because the temperature variation of thermal conductivity is dictated by phonon-phonon interactions. First, we optimized the EIPs to match the phonon dispersion following Lindsay and Broido.

^{37}Then, we used the optimized potential to calculate the temperature-dependent TC, following the procedure described in Appendix A. We were most interested in the temperature range between 100 K and 500 K because most devices operate in this range. A comparison between the calculated and measured thermal conductivities instructed us to further adjust the weights in the optimization process. This process was repeated until we obtained a satisfactory temperature-dependent thermal conductivity from the potentials.

For simplicity, we used existing parameterizations as starting points. To choose the best starting point, we calculate the temperature-dependent TC with existing potentials. For the three binary III-As materials (AlAs, GaAs, and InAs) we are most interested in, the potentials that yield the temperature-dependent thermal conductivity closest to the experiment are from the study by Sayed *et al.*,^{34} Powell *et al.*,^{35} and Hammerschmidt *et al.*,^{36} respectively. Incidentally, all three sets are parameterized in the Albe-Tersoff form. Like Lindsay and Broido,^{37} we tried to adjust just a few parameters. *R* and *D* were left untouched because they dictate the cutoff of the EIP. *m* also stayed constant, as per LAMMPS's requirement. Among the remaining parameters, we found that $De,\u2009\beta ,\u2009c,\u2009d,$ and *h* are very effective in adjusting the four physical properties we want to optimize. Therefore, during the $\chi 2$ minimization, we only varied these five parameters.

### A. Optimized potentials

In the optimization, we found that the parameters for GaAs from the study by Powell *et al.*^{35} yielded very good sound velocities and temperature-dependent thermal conductivity from 100 K to 500 K. Therefore, we adopted this set of parameters as it was. However, the parameters for both AlAs^{34} and InAs^{36} had to be optimized. Table I shows the optimized parameters for AlAs and InAs. As mentioned above, the parameters other than the five listed were kept the same as in the original sets.

AlAs . | InAs . |
---|---|

D = 2.6372 _{e} | D = 1.9949 _{e} |

β = 1.6948 | β = 1.7660 |

c = 1.4145 | c = 4.0249 |

d = 0.9116 | d = 1.0157 |

h = –0.6172 | h = –0.6096 |

AlAs . | InAs . |
---|---|

D = 2.6372 _{e} | D = 1.9949 _{e} |

β = 1.6948 | β = 1.7660 |

c = 1.4145 | c = 4.0249 |

d = 0.9116 | d = 1.0157 |

h = –0.6172 | h = –0.6096 |

The three panels in Fig. 1 show the calculated temperature-dependent thermal conductivity (dark green dots with error bars) in comparison to experimental data^{42–47} (light green solid lines) for AlAs, GaAs, and InAs, respectively. The insets show the calculated phonon dispersions along the $\Gamma \u2212X$ direction (dashed lines) in comparison to measured values^{48–51} (dots).

We note that the optical-phonon dispersions in Fig. 1 are not very accurate. Unfortunately, it is a known feature of short-range potentials that a simultaneous optimization of optical and acoustic branches proves problematic.^{37,52} Our purpose in this work, however, is thermal transport rather than phonon dispersions *per se*. Zone-center acoustic modes have an outsize direct effect on thermal transport near room temperature, owing to their high populations and long lifetimes (both stemming from their low energies), as well as their high group velocities. Therefore, the zone-center acoustic phonon spectrum is weighted heavily in the optimization, at the expense of the optical-phonon branches, which are of relatively lower direct importance for thermal transport at room temperature. (Optical phonons play an important indirect role in thermal transport through their participation in phonon–phonon scattering. Ideally, their dispersion should be accurately modeled.)

The ease with which the short-range binary potentials can be employed in alloys, along with an optimization procedure that focuses on the part of the phonon spectrum that is the most important for thermal transport (zone-center acoustic branches) and on the temperature-dependent thermal conductivity, gives us confidence that our optimized binary Tersoff potentials give physically meaningful results for thermal transport in the binaries and ternary alloys (in Sec. III) despite producing the less-than-ideal optical-phonon dispersions in the binaries.

## III. THERMAL CONDUCTIVITY OF TERNARY ALLOYS: THE ROLE OF THE CATION MASS DIFFERENCE

After obtaining the optimized EIPs for binary compounds, we applied these potentials to ternary alloys. We investigated the validity of the first assumption in the VCA (all cation atoms are assigned the same weighted-average mass) in predicting the TC. For this part, we kept the second assumption that all ternary alloys were random alloys. We considered In_{x}Ga_{1–x}As and In_{x}Al_{1–x}As alloys with *x* varying from 0.1 to 0.9. For each *x*, we generated 10 different random simulation cells with the corresponding composition. For each random configuration, we conducted 20 simulations with different starting velocity distributions. Each of the final thermal conductivity value was averaged over the 200 runs.

In order to isolate the effect of the mass difference, we carried out two sets of simulations (all at room temperature). In the first set, all the cation atoms kept their own masses (m_{Al} = 26.98 au, m_{Ga} = 69.72 au, and m_{In} = 114.82 au). In the second set, the cations were all assigned the same weighted-average VCA mass m_{avg} = *x*m_{In} + (1 – *x*)m_{Ga/Al}. Other than this difference, the two sets used the same EIPs and random (not ordered) configurations in the simulation (i.e., when we perform the explicit-mass calculation, the bonds are attributed potentials of a certain cation-ion pair; these potentials remain in place as we follow up with the VCA calculation, in which the explicit cation masses are replaced by the averaged VCA mass).

Figure 2 shows the thermal conductivity data for both In_{x}Ga_{1–x}As and In_{x}Al_{1–x}As. The results from the explicit-mass case and the averaged-mass VCA case are shown in blue diamonds and red dots, respectively. Existing experimental data are shown in stars, for comparison.^{10,11,15}

From Fig. 2, it is found that the explicit-mass thermal conductivity is consistently lower than the VCA thermal conductivity across the range of different In% for both In_{x}Ga_{1–x}As and In_{x}Al_{1–x}As. Moreover, the difference between the explicit-mass and VCA is much more pronounced in In_{x}Al_{1–x}As (where the cation mass difference is larger) than in In_{x}Ga_{1–x}As (where the cation mass difference is smaller). Therefore, it is evident that mass-difference scattering of lattice waves is critical in III-V ternary random alloys. Using an averaged VCA mass for cations will underestimate the mass-difference scattering; the larger the mass difference, the greater the underestimate. Additionally, in the case of In_{x}Ga_{1–x}As, both the explicit-mass thermal conductivity and the VCA thermal conductivity are significantly lower than the experimental value; in the case of In_{x}Al_{1–x}As, the two experimental results are so far apart that the explicit-mass thermal conductivity and VCA thermal conductivity fall in between the two measurements. As a result, we can conclude that there must be another mechanism that competes with mass-difference scattering to influence the thermal conductivity of ternary alloys.

## IV. THERMAL CONDUCTIVITY OF TERNARY ALLOYS: THE ROLE OF ORDER

Experimentalists who grow and characterize III-V epitaxial layers discovered long ago that both long-range and short-range order exist in ternary III-V alloys.^{20,22–31,53} Order leads to changes in the bond length and electronic bandgaps. Since phonons are quanta of lattice waves,^{54} it is intuitive that thermal transport will also change in the presence of ordering. Duda *et al.*^{55,56} studied one particular type of order in Si_{0.5}Ge_{0.5} using nonequilibrium molecular dynamics (NEMD) with a simple Lennard-Jones potential. Baker and Norris^{57} later studied both long-range order and short-range order in Si_{0.5}Ge_{0.5} with Stillinger-Weber potentials. We were interested in III-V ternary alloys with compositions away from the 50%–50% case. We also wanted to use EIPs that would give us quantitatively accurate results for TC. It is noteworthy that order in III-V ternary alloy samples is often localized, i.e., it is common to have “poly-ordering,” where the sample has one type of order in one region and a different type of order in another.^{21,24,53}

For the sake of simplicity, in this work, we only simulated samples with a single type of order and we applied periodic boundary conditions. Moreover, we focused on three different types of order that are most commonly observed in III-V ternary alloys: the CuPt-B type, the CuAu-I type, and the triple-period-A (TPA) type.^{26} Both CuPt-B and order CuAu-I order yield alloys with the composition A_{0.5}B_{0.5}C within the zinc-blende lattice. Figure 3 shows the crystal structures for CuPt-B order and CuAu-I order. In CuPt-B ordering, A atoms and B atoms take alternate cation planes along the [$1\xaf$11] direction. In CuAu-I ordering, A atoms and B atoms take alternate cation planes along the [100] direction. Note that perfectly ordered alloys do not exist in experiments, and so, the alternating planes are actually A-rich and B-rich planes. TPA order in its ideal form refers to the case where the cation planes along the [111] direction have a repeated pattern involving 3 planes. Therefore, ternary alloy A_{x}B_{1–x}C with any composition *x* can be represented in some TPA ordering where each period has the arrangement of A_{u}B_{1–u}/A_{v}B_{1–v}/A_{w}B_{1–w} and $(u+v+w)/3=x$ (note that the three planes in a period cannot be all equal or the triple period collapses). Figure 4 depicts a typical crystal structure of ternary alloy A_{0.5}B_{0.5}C, where *u* = 1, *v* = 0, and *w* = 0.5.

### A. Thermal conductivity of alloys near the A_{0.5}B_{0.5}C composition with CuPt-B and CuAu-I order

We calculated the room-temperature thermal conductivity of perfectly ordered In_{0.5}Ga_{0.5}As and In_{0.5}Al_{0.5}As with both CuPt-B and CuAu-I types of order; the results are shown in Table II. The calculated thermal conductivity for In_{0.5}Ga_{0.5}As and In_{0.5}Al_{0.5}As with random alloy structures (with explicit mass) is also shown for comparison. Since there are no experiments on In_{0.5}Ga_{0.5}As and In_{0.5}Al_{0.5}As, we compared the results with experiments on In_{0.53}Ga_{0.47}As (Ref. 10) and In_{0.52}Al_{0.48}As.^{11} From Table II, we concluded that the order of either type significantly increases the thermal conductivity of both In_{0.5}Ga_{0.5}As and In_{0.5}Al_{0.5}As alloys compared to the random-alloy case. The CuAu-I type order leads to slightly higher thermal conductivity than the CuPt-B type order in both alloys. All the simulated thermal conductivity values from the perfectly ordered structures were higher than experimental measurements, which is intuitive because In_{0.53}Ga_{0.47}As and In_{0.52}Al_{0.48}As in experiments are (1) not perfectly ordered and (2) have an In content different from 50%.

Material . | CuPt-B . | CuAu-I . | Rand. Alloy . | Expt. . |
---|---|---|---|---|

In_{0.5}Ga_{0.5}As | 13 | 14 | 1.5 | 4.84 (In_{0.53}Ga_{0.47}As) |

In_{0.5}Al_{0.5}As | 6.5 | 6.5 | 1.25 | 2.68 (In_{0.52}Al_{0.48}As) |

Material . | CuPt-B . | CuAu-I . | Rand. Alloy . | Expt. . |
---|---|---|---|---|

In_{0.5}Ga_{0.5}As | 13 | 14 | 1.5 | 4.84 (In_{0.53}Ga_{0.47}As) |

In_{0.5}Al_{0.5}As | 6.5 | 6.5 | 1.25 | 2.68 (In_{0.52}Al_{0.48}As) |

To directly compare the simulated thermal conductivity with the experimentally measured thermal conductivity for In_{0.53}Ga_{0.47}As and In_{0.52}Al_{0.48}As lattice matched to InP, we randomly replace 3% (2%) of Ga (Al) atoms in the perfectly ordered structure with In atoms to obtain lattice-matched simulation cells. Table III shows the thermal conductivity of lattice-matched In_{0.53}Ga_{0.47}As and In_{0.52}Al_{0.48}As from both our simulations and experiments. We see that the thermal conductivity obtained directly from the lattice-matched simulation cells is still quite large compared to the experiment, which could be attributed to the existence of additional disorder that is inevitable in real experimental samples. Here, we investigate one type of disorder: randomness. We take the lattice-matched simulation cells and randomly swap an In atom with a Ga or an Al atom to create disorder. The amount of disorder is categorized by the percentage of swapped In atoms inside the cell. In Table III, the cells with the simulation results are labeled with LM + *d*%, where LM stands for lattice matched and *d*% is the level of disorder. Consistent with the findings in Sec. III, even little disorder (<10%) severely reduces the alloy TC. Comparing the calculation with experiments, the level of disorder in In_{0.53}Ga_{0.47}As is around 2%, while in In_{0.52}Al_{0.48}As, it is approximately 8%. Also, as expected, the calculated lattice-matched thermal conductivity for InGaAs and InAlAs is lower than the corresponding thermal conductivity calculated from the perfectly ordered structures.

. | CuPt-B . | . | . | CuAu-I . | . | . | . |
---|---|---|---|---|---|---|---|

Material . | LM . | LM + 4% . | LM + 8% . | LM . | LM + 4% . | LM + 8% . | Expt. . |

In_{0.53}Ga_{0.47}As | 5.0 | 3.2 | 2.6 | 5.7 | 3.6 | 2.8 | 4.84 |

In_{0.52}Al_{0.48}As | 4.0 | 3.2 | 2.6 | 4.4 | 3.2 | 2.8 | 2.68 |

. | CuPt-B . | . | . | CuAu-I . | . | . | . |
---|---|---|---|---|---|---|---|

Material . | LM . | LM + 4% . | LM + 8% . | LM . | LM + 4% . | LM + 8% . | Expt. . |

In_{0.53}Ga_{0.47}As | 5.0 | 3.2 | 2.6 | 5.7 | 3.6 | 2.8 | 4.84 |

In_{0.52}Al_{0.48}As | 4.0 | 3.2 | 2.6 | 4.4 | 3.2 | 2.8 | 2.68 |

### B. Thermal conductivity of alloys with different compositions and TPA order

To study the thermal conductivity of ordered In_{x}Ga_{1–x}As and In_{x}Al_{1–x}As with various *x* (away from 0.5), we needed to implement TPA order. While it was impossible to consider every feasible arrangement of the TPA order, we simulated the two extreme cases that would likely yield the upper and lower limits of the thermal conductivity with TPA order for each *x*. The idea is that with more symmetry comes less phonon scattering, and the resulting thermal conductivity should be higher. Note that here we only focused on the perfectly ordered structures, and therefore, *x* values were limited to multiples of $112$. For each *x*, we generated two TPA cells: most symmetric (MS) and least symmetric (LS). Figure 5 shows the thermal conductivity calculated from the MS cells (red squares) and LS cells (yellow dots) and measured in the experiment^{10,11} (blue diamonds). The results are consistent with our expectations. The MS thermal conductivity is higher than the LS thermal conductivity, while both are higher than the experimental data (simulated structures are perfectly ordered, while experimental samples contain disorder). The thermal conductivity values of MS and LS fall on top of each other when $x=112$ and $x=1112$ because the MS and LS structures are equivalent in these cases. The calculated thermal conductivity follows the general U-shape (TC is the lowest when *x* is close to 0.5 and increases as *x* approaches 0 or 1), except for a sudden thermal conductivity jump at certain *x* values ($x=n12,\u2009n=2,4,6,8,\u2009and\u200910$). The jumps exist because these fractions $n12$ are reducible, which leads to additional symmetry in the system.

## V. CONCLUSION

In this paper, we studied the thermal conductivity of III-V ternary alloys, in particular In_{x}Ga_{1–x}As and In_{x}Al_{1–x}As. Using equilibrium molecular dynamics (EMD), we investigated how the mass difference between cation atoms and the arrangement of cations affect thermal transport. Optimized Able-Tersoff EIPs were employed, owing to their demonstrated success in describing III-V binary materials and their nearest-neighbor cutoffs that lend themselves to application in alloys.

We optimized the Able-Tersoff potentials for the binary compounds by matching the calculation to experimentally obtained phonon dispersions (with particular emphasis on zone-center acoustic phonons) and temperature-dependent thermal conductivity (to ensure that the potentials capture higher orders in the phonon–phonon interaction which are key in thermal transport). The optimized potentials for the binaries were then used to describe ternary alloys with both random and ordered structures. For random alloys, we compared the cases where atom masses were explicitly considered versus where they were all replaced with the averaged VCA mass, as is commonly done within the Boltzmann transport framework. The results showed that explicit atomic mass drastically reduces the thermal conductivity of ternary alloys. The larger the mass difference between the cations in the alloy, the larger the discrepancy between explicit-mass and averaged-mass VCA thermal conductivity results. We concluded that, when cation masses differ a great deal (as is the case for In_{x}Al_{1–x}As), it is essential to include atomic masses explicitly and any calculation that relies on VCA is likely inaccurate.

Moreover, measured thermal conductivities of ternary alloys are higher than calculations with either explicit or VCA mass and random positioning of cations, which led us to look at longer-range order in alloys. We considered perfectly ordered In_{0.5}Ga_{0.5}As and In_{0.5}Al_{0.5}As with both CuPt-B and CuAu-I types of order and the corresponding lattice-matched alloys In_{0.53}Ga_{0.47}As and In_{0.52}Al_{0.48}As. The order in ternary alloys considerably raises TC. By adding random disorder to the lattice-matched alloys In_{0.53}Ga_{0.47}As and In_{0.52}Al_{0.48}As, we found that experimental results could be reproduced with levels of disorder close to 2% in In_{0.53}Ga_{0.47}As and 8% in In_{0.52}Al_{0.48}As.

We also studied perfectly ordered TPA alloys In_{x}Ga_{1–x}As and In_{x}Al_{1–x}As with various compositions *x*. We found that more symmetry in the alloy led to higher TC, while the alloys with the least symmetry still yielded higher thermal conductivity than experiments, indicating the existence of disorder in the experiment.

In conclusion, in modeling thermal transport in III-V ternary alloys, it is crucial to include both the explicit masses of atoms and the effects of long-range order. The measured thermal conductivity for III-V ternary alloys is likely a result of the competition between the two: reduction in thermal conductivity stemming from mass-difference scattering associated with random disorder and an increase in thermal conductivity associated with order in the alloy structure. These notions should be incorporated into other techniques for calculating thermal transport in alloys and highlight the importance of critically evaluating the range of validity of even very common approximations, such as the VCA.

## ACKNOWLEDGMENTS

The authors thank Luke Mawst for bringing the issue of long-range order in III-V alloys to our attention. We are also grateful to Dan Botez, Mark Eriksson, Jeremy Kirch, Gabriel Jaffe, and Colin Boyle for useful discussions. This work was supported by the U.S. Department of Energy, Office of Science (Basic Energy Sciences, Division of Materials Sciences and Engineering, Physical Behavior of Materials Program) under Award No. DE-SC0008712. The simulation work was performed using the compute resources and assistance of the UW-Madison Center For High Throughput Computing (CHTC) in the Department of Computer Sciences.

### APPENDIX A: CALCULATION OF THERMAL CONDUCTIVITY FROM EQUILIBRIUM MOLECULAR DYNAMICS

##### 1. The Green-Kubo formula for thermal conductivity

All the MD simulations in this work were carried out using the LAMMPS^{38} package. We used equilibrium molecular dynamics (EMD) together with the Green-Kubo formula to compute the thermal conductivity of compound semiconductors. When a system is equilibrated and then kept at a set temperature, the thermal conductivity can be derived from the fluctuation–dissipation theorem using instantaneous heat flux.^{58} In a three-dimensional (3D) bulk material where all three directions are equivalent, thermal conductivity *κ* is described by the Green-Kubo formula as^{58}

where *k*_{B} is the Boltzmann constant, *V* is the system volume, *T* is the system temperature, and **S**(*t*) is the instantaneous heat flux calculated from the atom velocities and local potentials. The running integral over the heat-current autocorrelation function $\u27e8S(t)S(0)\u27e9$ is depicted in Fig. 6. In an MD simulation, the time is discretized, and the integration in Eq. (A1) is converted to a sum.

For cubic bulk semiconductors, the thermal conductivity is expected to be isotropic. As a result, we used cubic simulation cells with periodic boundary conditions applied in all directions. In a typical simulation, the system was first initialized at the desired temperature *T* by assigning each atom a random velocity that follows the thermal distribution at *T*. Then, the system was equilibrated as an *NPT* ensemble (constant number of atoms, constant pressure at 1 atm, and constant temperature *T*) using the Nosé-Hoover barostats and thermostats^{59,60} for 100 ps. After that, the system was further equilibrated as an *NVE* ensemble (constant number of atoms, constant volume, and constant system energy) for another 100 ps before the heat flux is collected. The instantaneous heat flux was calculated and output into a file at every time step, for 5 × 10^{6} steps. A script was written to post-process the output file and obtain the heat-current autocorrelation function and its running integral. We made sure the integral saturated and extracted the bulk thermal conductivity according to Eq. (A1). The time step was chosen to be small enough to ensure the system stayed stable throughout the simulation—here, the time step was 1 fs for binary materials and 0.1 fs for ternary alloys. (The environment in random ternaries is less stable than in the binaries, and the small time step ensures stability. A larger time step could be chosen for more ordered ternary structures, but for the sake of consistency, we simulate all ternaries with the same short time step.)

For each simulation (given the material and temperature), several random starting velocity distributions were used, and the final result was averaged among different runs. We tested the simulation domain sizes to make sure that there was no size effect. The final results for binary III-As were obtained with a simulation cell $10a0\xd710a0\xd710a0$ in size, where *a*_{0} is the lattice constant for the material. For ternary alloys, $8a0\xd78a0\xd78a0$ is enough for the thermal conductivity to converge. We used $9a0\xd79a0\xd79a0$ cells for alloys with TPA ordering, as the algorithm for generating cells requires that the system size be a multiple of 3 (more details are in Sec. IV).

##### 2. Quantum correction of temperature

In MD simulations, system temperatures are calculated following the rules of classical statistical mechanics.^{61} We adopted a simple quantum-correction procedure for the temperatures by mapping the kinetic energy of an MD system at temperature *T*_{MD} onto that of a quantum system with temperature *T*_{Q}^{61,62}

On the left-hand side, *N* is the number of atoms in the system and *k*_{B} is the Boltzmann constant. On the right-hand side, the summation is over all the phonon branches b and wavevectors **k**. $\u210f\omega (k,b)$ is the corresponding phonon energy. For the right-hand side, we used an approximate isotropic phonon dispersion fitted to the full phonon dispersion, based on our previous work.^{12} Figure 7 shows the mapping between the quantum-corrected temperature and the MD temperature for GaAs (the curve is similar for AlAs and InAs) between 0 K and 500 K. We see that *T*_{MD} and *T*_{Q} coincide at higher temperatures but differ at lower temperatures. At room temperature, there is typically a 4%–8% difference between *T*_{MD} and *T*_{Q} for III-As. The temperatures listed in this work are the quantum-corrected temperatures *T*_{Q}. However, we note that MD is a decidedly classical technique, applicable above the Debye temperature.

The quantum-corrected temperature is employed here because we optimize interatomic potentials against experimentally obtained temperature-dependent conductivity; the temperature correction enables us to formally go below the Debye temperature for fitting purposes. However, we note that the quantum-correction technique has limited applicability, which has been recognized in the literature.^{63}

### APPENDIX B: TERSOFF POTENTIALS

The Tersoff EIP was proposed by Tersoff in 1988 for silicon,^{32} extended to SiC in 1989,^{33} and later parameterized for III-V binary compounds.^{34–36,39} Tersoff potentials have a cutoff distance that limits the atom interactions to nearest neighbors, which is an advantage in our case: the nearest-neighbor interaction makes it easy to apply the optimized potentials for binary III-As to ternary III-As alloys.

In the LAMMPS package, the Tersoff potential is described using 14 parameters in the following form. The total energy of the system is the summation of energy between each pair

The pair interaction potential *V _{ij}* between atom

*i*and atom

*j*is described by the competition between a repulsive term

*f*and an attractive term

_{R}*f*and is modulated by a cutoff term

_{A}*f*so that only nearest-neighbor interactions are included

_{C}The cutoff function is

where *r* is the variable. From Eq. (B2b), the cutoff length is *R *+* D*, and there is a transition window with width 2*D*. Both *R* and *D* are the parameters of the Tersoff potential. Both the repulsive and the attractive terms have an exponential form with *A*, *B*, *λ*_{1}, and *λ*_{2} being the potential parameters

To include the influence of the bond angle and length, the attraction term *f _{A}* is further modified by the bond angle term

*b*, where

_{ij}and

The summation in *ξ _{ij}* goes over all the central atom

*i*'s neighbors within the cutoff except for the neighbor

*j*whose interaction with

*i*we are considering. The bond angle

*θ*can be calculated with the three atoms' coordinates, and

_{ijk}*δ*is a scaling parameter that is unity in most cases.

_{ijk}*γ*,

*n*,

*λ*

_{3}, and

*m*are the potential parameters. In particular,

*m*can only take the value of 1 or 3.

*c*,

*d*, and $cos\u2009\theta 0$ are the bond-angle-related potential parameters. Note that $cos\u2009\theta 0$ is often denoted as

*h*; the form of a cosine simply reminds us that this parameter can only take values between −1 and 1.

A variation of the original Tersoff potential is the Albe-Tersoff potential where the strength and the decay rate in Eqs. (B2c) and (B2d) are expressed through parameters *D _{e}*,

*S*,

*β*, and

*r*in place of

_{e}*A*,

*B*,

*λ*

_{1}, and

*λ*

_{2}.

^{64}The relationship between the two sets of parameters can be derived as

In the case of Albe-Tersoff potential, we work with the parameters *D _{e}*,

*S*,

*β*, and

*r*and then use a script to convert them to

_{e}*A*,

*B*,

*λ*

_{1}, and

*λ*

_{2}for the input potential file of LAMMPS.

### APPENDIX C: OBTAINING PHONON DISPERSION FROM EIP

With a given EIP, the phonon frequencies $\omega (k,b)$ for wavevector **k** and branch b can be obtained by diagonalizing the dynamical matrix^{65}

where $l\u2032j$ represents atom *j* in the $l\u2032$ th unit cell, and the summation is over all the relevant neighbors of the central atom *i*. $x(l\u2032)$ is the relative location of unit cell $l\u2032$ with respect to the unit cell that atom *i* is in. The interatomic force constants (IFCs) between atoms *i* and *j*, denoted by $\varphi \alpha \beta (0i;l\u2032j)$, were calculated using centred finite differences

Phonon dispersion curves were obtained by computing the phonon frequencies for multiple **k**s along a certain direction in the first Brillouin zone (1BZ). The sound velocities were obtained from the acoustic branches near the zone center (Γ point) as

where the subscript b is the branch index. Since the phonon dispersion is isotropic near the Γ point, we used the dispersion curves along the [100] (Γ – X) direction to calculate the sound velocities. The two transverse acoustic (TA) branches are degenerate along the Γ – X direction, and so, we only need one scalar velocity for the TA branches (*v _{s}*

_{,TA}) and one for the longitudinal acoustic (LA) branch (

*v*

_{s}_{,LA}).

## References

_{1–x}In

_{x}As Thin Films

_{0.5}In

_{0.5}As Grown by Metalorganic-Vapor-Phase-Epitaxy