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 In0.53Ga0.47As and In0.52Al0.48As, 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.

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 AlxGa1–xAs with any Al composition, the thermal conductivity measurement of AlGaAs was carried out early on, though only at room temperature.6 Measurements on InxGa1–xAs and InxAl1–xAs are more challenging, as they each only have a specific composition (In0.53Ga0.47As and In0.52Al0.48As) that is lattice matched to InP, a common substrate for the growth of mid-infrared quantum cascade lasers.7 As a result, InxGa1–xAs and InxAl1–xAs 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 InxGa1–xAs with various In compositions, where the samples were formed through zone leveling.10 However, this method is not applicable to InxAl1–xAs, as they oxidize easily. Koh et al. measured the room-temperature thermal conductivity of In0.53Ga0.47As and In0.52Al0.48As 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 In0.53Ga0.47As/In0.52Al0.48As SLs with different layer thicknesses and extracted the thermal conductivity of In0.53Ga0.47As and In0.52Al0.48As 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 In0.52Al0.48As. Furthermore, both sets of results were obtained for In0.53Ga0.47As and In0.52Al0.48As lattice matched to InP. We still have very limited knowledge about the thermal conductivity of InxGa1–xAs and InxAl1–xAs 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 approximation13 (VCA) to compute the thermal conductivity of ternary alloys AlxGa1–xAs, InxGa1–xAs, and InxAl1–xAs at room temperature. For InxGa1–xAs and AlxGa1–xAs, our calculated thermal conductivity agreed reasonably well with experimental values with various x compositions.6,10,14 For InxAl1–xAs, an experimental thermal conductivity value was only available for In0.52Al0.48As.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 AxB1–xC. (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., mIII = xmA + (1 – x)mB. (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, mAl = 26.98 au, mGa = 69.72 au, and mIn = 114.82 au, and so, it makes sense that the model would work better for InxGa1–xAs and AlxGa1–xAs than it does for InxAl1–xAs 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 In0.53Ga0.47As and In0.52Al0.48As are the CuPt-B order20–22,24 and the triple-period-A (TPA) order.23,25–29 The CuAu-I order is more common in Al0.5Ga0.5As,30,31 but it could also be observed in In0.53Ga0.47As or In0.52Al0.48As 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 Tersoff32,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 In0.53Ga0.47As and In0.52Al0.48As that contain small levels of random disorder. Finally, we summarize our findings in Sec. V.

All the MD simulations in this work were carried out in the LAMMPS38 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 Broido37 optimized the Tersoff potential of carbon for thermal transport. They used a χ2-minimization, where

χ2=i(ηiηexp,i)2ηexp,i2ζi,
(1)

where i runs over all the physical properties they optimize for; ηi and ηexp,i are the calculated and experimental values of the physical property i, respectively; and ζi 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.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 a0, the cohesive energy Ecoh, and sound velocities vs,TA and vs,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,β,c,d, and h are very effective in adjusting the four physical properties we want to optimize. Therefore, during the χ2 minimization, we only varied these five parameters.

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 AlAs34 and InAs36 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.

TABLE I.

Optimized parameters for the Albe-Tersoff potentials for AlAs and InAs. Parameters not listed are the same as in the original sets.34,36

AlAsInAs
D e = 2.6372 De = 1.9949 
β = 1.6948 β = 1.7660 
c = 1.4145 c = 4.0249 
d = 0.9116 d = 1.0157 
h = –0.6172 h = –0.6096 
AlAsInAs
D e = 2.6372 De = 1.9949 
β = 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 data42–47 (light green solid lines) for AlAs, GaAs, and InAs, respectively. The insets show the calculated phonon dispersions along the ΓX direction (dashed lines) in comparison to measured values48–51 (dots).

FIG. 1.

Temperature-dependent thermal conductivity, as calculated in this work with optimized potentials (dark green dots with error bars) and from the experiment (light green solid lines) for AlAs,44 GaAs,42,43 and InAs.45–47 The insets show the calculated phonon dispersions (dashed lines) along the Γ – X direction along with experimental values (dots).48–51 

FIG. 1.

Temperature-dependent thermal conductivity, as calculated in this work with optimized potentials (dark green dots with error bars) and from the experiment (light green solid lines) for AlAs,44 GaAs,42,43 and InAs.45–47 The insets show the calculated phonon dispersions (dashed lines) along the Γ – X direction along with experimental values (dots).48–51 

Close modal

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.

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 InxGa1–xAs and InxAl1–xAs 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 (mAl = 26.98 au, mGa = 69.72 au, and mIn = 114.82 au). In the second set, the cations were all assigned the same weighted-average VCA mass mavg = xmIn + (1 – x)mGa/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 InxGa1–xAs and InxAl1–xAs. 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

FIG. 2.

Thermal conductivity as a function of the In composition in InxGa1–xAs and InxAl1–xAs random alloys at room temperature. In blue diamonds, we show the calculated thermal conductivity for the case where each atom's mass is included explicitly. The results from the case where all cations take on the averaged (VCA) mass are shown in red dots. Stars show the experimental data.10,11,15

FIG. 2.

Thermal conductivity as a function of the In composition in InxGa1–xAs and InxAl1–xAs random alloys at room temperature. In blue diamonds, we show the calculated thermal conductivity for the case where each atom's mass is included explicitly. The results from the case where all cations take on the averaged (VCA) mass are shown in red dots. Stars show the experimental data.10,11,15

Close modal

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 InxGa1–xAs and InxAl1–xAs. Moreover, the difference between the explicit-mass and VCA is much more pronounced in InxAl1–xAs (where the cation mass difference is larger) than in InxGa1–xAs (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 InxGa1–xAs, both the explicit-mass thermal conductivity and the VCA thermal conductivity are significantly lower than the experimental value; in the case of InxAl1–xAs, 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.

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 Si0.5Ge0.5 using nonequilibrium molecular dynamics (NEMD) with a simple Lennard-Jones potential. Baker and Norris57 later studied both long-range order and short-range order in Si0.5Ge0.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 A0.5B0.5C 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¯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 AxB1–xC with any composition x can be represented in some TPA ordering where each period has the arrangement of AuB1–u/AvB1–v/AwB1–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 A0.5B0.5C, where u = 1, v = 0, and w = 0.5.

FIG. 3.

The crystal structure of ordered ternary alloy A0.5B0.5C with CuPt-B and CuAu-I ordering.

FIG. 3.

The crystal structure of ordered ternary alloy A0.5B0.5C with CuPt-B and CuAu-I ordering.

Close modal
FIG. 4.

A sample crystal structure of ordered ternary alloy A0.5B0.5C with TPA ordering. The 3 planes in a period (Au B1–u/Av B1–v/Aw B1–w) have configurations of u = 1, v = 0, and w = 0.5.

FIG. 4.

A sample crystal structure of ordered ternary alloy A0.5B0.5C with TPA ordering. The 3 planes in a period (Au B1–u/Av B1–v/Aw B1–w) have configurations of u = 1, v = 0, and w = 0.5.

Close modal

We calculated the room-temperature thermal conductivity of perfectly ordered In0.5Ga0.5As and In0.5Al0.5As with both CuPt-B and CuAu-I types of order; the results are shown in Table II. The calculated thermal conductivity for In0.5Ga0.5As and In0.5Al0.5As with random alloy structures (with explicit mass) is also shown for comparison. Since there are no experiments on In0.5Ga0.5As and In0.5Al0.5As, we compared the results with experiments on In0.53Ga0.47As (Ref. 10) and In0.52Al0.48As.11 From Table II, we concluded that the order of either type significantly increases the thermal conductivity of both In0.5Ga0.5As and In0.5Al0.5As 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 In0.53Ga0.47As and In0.52Al0.48As in experiments are (1) not perfectly ordered and (2) have an In content different from 50%.

TABLE II.

Calculated room-temperature thermal conductivity for In0.5Ga0.5As and In0.5Al0.5As with perfectly ordered CuPt-B and CuAu-I structures and for a random alloy. The experimentally measured thermal conductivity for In0.53Ga0.47As (Ref. 10) and In0.52Al0.48As (Ref. 11) is given for comparison.

MaterialCuPt-BCuAu-IRand. AlloyExpt.
In0.5Ga0.5As 13 14 1.5 4.84 (In0.53Ga0.47As) 
In0.5Al0.5As 6.5 6.5 1.25 2.68 (In0.52Al0.48As) 
MaterialCuPt-BCuAu-IRand. AlloyExpt.
In0.5Ga0.5As 13 14 1.5 4.84 (In0.53Ga0.47As) 
In0.5Al0.5As 6.5 6.5 1.25 2.68 (In0.52Al0.48As) 

To directly compare the simulated thermal conductivity with the experimentally measured thermal conductivity for In0.53Ga0.47As and In0.52Al0.48As 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 In0.53Ga0.47As and In0.52Al0.48As 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 In0.53Ga0.47As is around 2%, while in In0.52Al0.48As, 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.

TABLE III.

Calculated room-temperature thermal conductivity for lattice-matched (LM) In0.53Ga0.47As and In0.52Al0.48As with various percentages of additional random disorder. The experimentally measured thermal conductivity values are also listed for comparison.10,11

CuPt-BCuAu-I
MaterialLMLM + 4%LM + 8%LMLM + 4%LM + 8%Expt.
In0.53Ga0.47As 5.0 3.2 2.6 5.7 3.6 2.8 4.84 
In0.52Al0.48As 4.0 3.2 2.6 4.4 3.2 2.8 2.68 
CuPt-BCuAu-I
MaterialLMLM + 4%LM + 8%LMLM + 4%LM + 8%Expt.
In0.53Ga0.47As 5.0 3.2 2.6 5.7 3.6 2.8 4.84 
In0.52Al0.48As 4.0 3.2 2.6 4.4 3.2 2.8 2.68 

To study the thermal conductivity of ordered InxGa1–xAs and InxAl1–xAs 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 experiment10,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,n=2,4,6,8,and10). The jumps exist because these fractions n12 are reducible, which leads to additional symmetry in the system.

FIG. 5.

Calculated room-temperature thermal conductivity of InGaAs and InAlAs with perfect TPA ordering and various In compositions. Red squares and yellow dots show the results from the most symmetric (MS) and least symmetric (LS) structures. Experimental measurements10,11 are shown in blue diamonds.

FIG. 5.

Calculated room-temperature thermal conductivity of InGaAs and InAlAs with perfect TPA ordering and various In compositions. Red squares and yellow dots show the results from the most symmetric (MS) and least symmetric (LS) structures. Experimental measurements10,11 are shown in blue diamonds.

Close modal

In this paper, we studied the thermal conductivity of III-V ternary alloys, in particular InxGa1–xAs and InxAl1–xAs. 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 InxAl1–xAs), 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 In0.5Ga0.5As and In0.5Al0.5As with both CuPt-B and CuAu-I types of order and the corresponding lattice-matched alloys In0.53Ga0.47As and In0.52Al0.48As. The order in ternary alloys considerably raises TC. By adding random disorder to the lattice-matched alloys In0.53Ga0.47As and In0.52Al0.48As, we found that experimental results could be reproduced with levels of disorder close to 2% in In0.53Ga0.47As and 8% in In0.52Al0.48As.

We also studied perfectly ordered TPA alloys InxGa1–xAs and InxAl1–xAs 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.

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.

1. The Green-Kubo formula for thermal conductivity

All the MD simulations in this work were carried out using the LAMMPS38 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 as58 

κ=13kBVT20S(t)S(0)dt,
(A1)

where kB 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 S(t)S(0) is depicted in Fig. 6. In an MD simulation, the time is discretized, and the integration in Eq. (A1) is converted to a sum.

FIG. 6.

Running integral of the autocorrelation function versus time, Eq. (A1), multiplied by the appropriate constants to yield the thermal conductivity in the steady state. The results are presented for a random InGaAs alloy with 50% In and explicit mass in the MD calculation at room temperature.

FIG. 6.

Running integral of the autocorrelation function versus time, Eq. (A1), multiplied by the appropriate constants to yield the thermal conductivity in the steady state. The results are presented for a random InGaAs alloy with 50% In and explicit mass in the MD calculation at room temperature.

Close modal

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 thermostats59,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 × 106 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×10a0×10a0 in size, where a0 is the lattice constant for the material. For ternary alloys, 8a0×8a0×8a0 is enough for the thermal conductivity to converge. We used 9a0×9a0×9a0 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 TMD onto that of a quantum system with temperature TQ61,62

32NkBTMD=bkω(k,b){12+1exp[ω(k,b)kBTQ]1}.
(A2)

On the left-hand side, N is the number of atoms in the system and kB is the Boltzmann constant. On the right-hand side, the summation is over all the phonon branches b and wavevectors k. ω(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 TMD and TQ coincide at higher temperatures but differ at lower temperatures. At room temperature, there is typically a 4%–8% difference between TMD and TQ for III-As. The temperatures listed in this work are the quantum-corrected temperatures TQ. However, we note that MD is a decidedly classical technique, applicable above the Debye temperature.

FIG. 7.

A typical mapping between the MD-calculated temperature TMD and the quantum-corrected temperature TQ (figure shown for GaAs). TMD and TQ coincide at high temperatures but differ a great deal at low temperatures.

FIG. 7.

A typical mapping between the MD-calculated temperature TMD and the quantum-corrected temperature TQ (figure shown for GaAs). TMD and TQ coincide at high temperatures but differ a great deal at low temperatures.

Close modal

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 

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

E=12ijiVij.
(B1)

The pair interaction potential Vij between atom i and atom j is described by the competition between a repulsive term fR and an attractive term fA and is modulated by a cutoff term fC so that only nearest-neighbor interactions are included

Vij=fC(rij)[fR(rij)+bijfA(rij)].
(B2a)

The cutoff function is

fC(r)={1:r<RD,1212sin(π2rRD):RD<r<R+D,0:r>R+D,
(B2b)

where r is the variable. From Eq. (B2b), the cutoff length is R + D, and there is a transition window with width 2D. 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

fR(r)=Aexp(λ1r),
(B2c)
fA(r)=Bexp(λ2r).
(B2d)

To include the influence of the bond angle and length, the attraction term fA is further modified by the bond angle term bij, where

bij=(1+γnξijn)12n,
(B2e)

and

ξij=ki,jfC(rik)g(θijk)exp[λ3m(rijrik)m],
(B2f)
g(θ)=δijk(1+c2d2c2d2+(cosθcosθ0)2).
(B2g)

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 θijk can be calculated with the three atoms' coordinates, and δijk is a scaling parameter that is unity in most cases. γ, n, λ3, and m are the potential parameters. In particular, m can only take the value of 1 or 3. c, d, and cosθ0 are the bond-angle-related potential parameters. Note that cosθ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 De, S, β, and re in place of A, B, λ1, and λ2.64 The relationship between the two sets of parameters can be derived as

A=DeS1exp(β2Sre),
(B2h)
B=SDeS1exp(β2S/re),
(B2i)
λ1=β2Sre,
(B2j)
λ2=β2S/re.
(B2k)

In the case of Albe-Tersoff potential, we work with the parameters De, S, β, and re and then use a script to convert them to A, B, λ1, and λ2 for the input potential file of LAMMPS.

With a given EIP, the phonon frequencies ω(k,b) for wavevector k and branch b can be obtained by diagonalizing the dynamical matrix65 

Dαβ(ij|k)=1mimjlϕαβ(0i;lj)eikx(l),
(C1)

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

ϕαβij=Vij2αβ=V(rij+dα+dβ)V(rijdα+dβ)V(rij+dαdβ)+V(rijdαdβ)4dαdβ.
(C2)

Phonon dispersion curves were obtained by computing the phonon frequencies for multiple ks 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

|vs,b|=|kω(k,b)|,
(C3)

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 (vs,TA) and one for the longitudinal acoustic (LA) branch (vs,LA).

1.
J.
Faist
,
F.
Capasso
,
D. L.
Sivco
,
C.
Sirtori
,
A. L.
Hutchinson
, and
A. Y.
Cho
,
Science
264
,
553
(
1994
).
2.
R.
Köhler
,
A.
Tredicucci
,
F.
Beltram
,
H. E.
Beere
,
E. H.
Linfield
,
A. G.
Davies
,
D. A.
Ritchie
,
R. C.
Iotti
, and
F.
Rossi
,
Nature
417
,
156
(
2002
).
3.
M.
Razeghi
,
N.
Bandyopadhyay
,
Y.
Bai
,
Q.
Lu
, and
S.
Slivken
,
Opt. Mater. Express
3
,
1872
(
2013
).
4.
D.
Botez
,
C.-C.
Chang
, and
L. J.
Mawst
,
J. Phys. D: Appl. Phys.
49
,
043001
(
2016
).
5.
S.
Adachi
,
J. Appl. Phys.
102
,
063502
(
2007
).
6.
M. A.
Afromowitz
,
J. Appl. Phys.
44
,
1292
(
1973
).
7.
Y. M.
Kim
,
M. J. W.
Rodwell
, and
A. C.
Gossard
,
J. Electron. Mater.
31
,
196
(
2002
).
8.
D. G.
Cahill
,
W. K.
Ford
,
K. E.
Goodson
,
G. D.
Mahan
,
A.
Majumdar
,
H. J.
Maris
,
R.
Merlin
, and
S. R.
Phillpot
,
J. Appl. Phys.
93
,
793
(
2003
).
9.
D. G.
Cahill
,
P. V.
Braun
,
G.
Chen
,
D. R.
Clarke
,
S.
Fan
,
K. E.
Goodson
,
P.
Keblinski
,
W. P.
King
,
G. D.
Mahan
,
A.
Majumdar
,
H. J.
Maris
,
S. R.
Phillpot
,
E.
Pop
, and
L.
Shi
,
Appl. Phys. Rev.
1
,
011305
(
2014
).
10.
M.
Abrahams
,
R.
Braunstein
, and
F.
Rosi
,
J. Phys. Chem. Solids
10
,
204
(
1959
).
11.
Y. K.
Koh
,
S. L.
Singer
,
W.
Kim
,
J. M. O.
Zide
,
H.
Lu
,
D. G.
Cahill
,
A.
Majumdar
, and
A. C.
Gossard
,
J. Appl. Phys.
105
,
054303
(
2009
).
12.
S.
Mei
and
I.
Knezevic
,
J. Appl. Phys.
118
,
175101
(
2015
).
13.
14.
S.
Adachi
,
J. Appl. Phys.
54
,
1844
(
1983
).
15.
A.
Sood
,
J. A.
Rowlette
,
C. G.
Caneau
,
E.
Bozorg-Grayeli
,
M.
Asheghi
, and
K. E.
Goodson
,
Appl. Phys. Lett.
105
,
051909
(
2014
).
16.
Z.
Tian
,
J.
Garg
,
K.
Esfarjani
,
T.
Shiga
,
J.
Shiomi
, and
G.
Chen
,
Phys. Rev. B
85
,
184303
(
2012
).
17.
J.
Garg
,
N.
Bonini
,
B.
Kozinsky
, and
N.
Marzari
,
Phys. Rev. Lett.
106
,
045901
(
2011
).
18.
W.
Li
,
L.
Lindsay
,
D. A.
Broido
,
D. A.
Stewart
, and
N.
Mingo
,
Phys. Rev. B
86
,
174307
(
2012
).
19.
J. M.
Larkin
and
A. J. H.
McGaughey
,
J. Appl. Phys.
114
,
023507
(
2013
).
20.
M. A.
Shahid
,
S.
Mahajan
,
D. E.
Laughlin
, and
H. M.
Cox
,
Phys. Rev. Lett.
58
,
2567
(
1987
).
21.
K.
Shin
,
J.
Yoo
,
S.
Joo
,
T.
Mori
,
D.
Shindo
,
T.
Hanada
,
H.
Makino
,
M.
Cho
,
T.
Yao
, and
Y.-G.
Park
,
Mater. Trans.
47
,
1115
(
2006
).
22.
T. S.
Kuan
,
W. I.
Wang
, and
E. L.
Wilkie
,
Appl. Phys. Lett.
51
,
51
(
1987
).
23.
T.
Mori
,
T.
Hanada
,
T.
Morimura
,
K.
Shin
,
H.
Makino
, and
T.
Yao
,
Appl. Surf. Sci.
237
,
230
(
2004
).
24.
J.
Kulik
,
R.
Forrest
,
J.
Li
,
T.
Golding
,
S. C.
Moss
, and
J.
Bai
,
X-Ray and Tem Studies of Short-Range Order in Al1–xInxAs Thin Films
(
Mater. Res. Soc. Symp. Proc.
,
1999
), Vol.
583
, p. 179.
25.
R. L.
Forrest
,
J.
Kulik
,
T. D.
Golding
, and
S. C.
Moss
,
J. Mater. Res.
15
,
45
55
(
2000
).
26.
T.
Suzuki
, in
Spontaneous Ordering in Semiconductor Alloys
, edited by
A.
Mascarenhas
(
Springer
,
US, New York
,
2002
).
27.
T.
Suzuki
,
T.
Ichihashi
, and
T.
Nakayama
,
Appl. Phys. Lett.
73
,
2588
(
1998
).
28.
A.
Gomyo
,
K.
Makita
,
I.
Hino
, and
T.
Suzuki
,
Phys. Rev. Lett.
72
,
673
(
1994
).
29.
S.
Ohkouchi
,
T.
Furuhashi
,
A.
Gomyo
,
K.
Makita
, and
T.
Suzuki
,
Appl. Surf. Sci.
241
,
9
(
2005
).
30.
T. S.
Kuan
,
T. F.
Kuech
,
W. I.
Wang
, and
E. L.
Wilkie
,
Phys. Rev. Lett.
54
,
201
(
1985
).
31.
J. E.
Bernard
,
R. G.
Dandrea
,
L. G.
Ferreira
,
S.
Froyen
,
S.
Wei
, and
A.
Zunger
,
Appl. Phys. Lett.
56
,
731
(
1990
).
32.
33.
34.
M.
Sayed
,
J.
Jefferson
,
A.
Walker
, and
A.
Cullis
,
Nucl. Instrum. Methods Phys. Res., Sect. B
102
,
218
(
1995
).
35.
D.
Powell
,
M. A.
Migliorato
, and
A. G.
Cullis
,
Phys. Rev. B
75
,
115202
(
2007
).
36.
T.
Hammerschmidt
,
P.
Kratzer
, and
M.
Scheffler
,
Phys. Rev. B
77
,
235303
(
2008
).
37.
L.
Lindsay
and
D. A.
Broido
,
Phys. Rev. B
81
,
205441
(
2010
).
38.
S.
Plimpton
,
J. Comput. Phys.
117
,
1
(
1995
).
39.
K.
Nordlund
,
J.
Nord
,
J.
Frantz
, and
J.
Keinonen
,
Comput. Mater. Sci.
18
,
283
(
2000
).
40.
H.
Detz
and
G.
Strasser
,
Semicond. Sci. Technol.
28
,
085011
(
2013
).
41.
J. T.
Titantah
,
D.
Lamoen
,
M.
Schowalter
, and
A.
Rosenauer
,
J. Appl. Phys.
101
,
123508
(
2007
).
42.
A. V.
Inyushkin
,
A. N.
Taldenkov
,
A. Y.
Yakubovsky
,
A. V.
Markov
,
L.
Moreno-Garsia
, and
B. N.
Sharonov
,
Semicond. Sci. Technol.
18
,
685
(
2003
).
43.
A.
Amith
,
I.
Kudman
, and
E. F.
Steigmeier
,
Phys. Rev.
138
,
A1270
(
1965
).
44.
C. A.
Evans
,
D.
Indjin
,
Z.
Ikonic
,
P.
Harrison
,
M. S.
Vitiello
,
V.
Spagnolo
, and
G.
Scamarcio
,
IEEE J. Quantum Electron.
44
,
680
(
2008
).
45.
P. V.
Tamarin
and
S. S.
Shalyt
,
Sov. Phys. Semicond.
5
,
1097
(
1971
).
46.
G.
Le Guillou
and
H. J.
Albany
,
Phys. Rev. B
5
,
2301
(
1972
).
47.
R.
Bowers
,
R. W.
Ure
,
J. E.
Bauerle
, and
A. J.
Cornish
,
J. Appl. Phys.
30
,
930
(
1959
).
48.
L.
Lindsay
,
D. A.
Broido
, and
T. L.
Reinecke
,
Phys. Rev. B
87
,
165201
(
2013
).
49.
B.
Dorner
and
D.
Strauch
,
J. Phys.: Condens. Matter
2
,
1475
(
1990
).
50.
N. S.
Orlova
,
Phys. Status Solidi B
119
,
541
(
1983
).
51.
R.
Carles
,
N.
Saint-Cricq
,
J. B.
Renucci
,
M. A.
Renucci
, and
A.
Zwick
,
Phys. Rev. B
22
,
4804
(
1980
).
52.
P.
Han
and
G.
Bester
,
Phys. Rev. B
83
,
174304
(
2011
).
53.
T.
Suzuki
,
T.
Ichihasfh
, and
M.
Tsuji
,
Triple-Period (TP)-A and CuPt-A Type Ordering in Al0.5In0.5As Grown by Metalorganic-Vapor-Phase-Epitaxy
(
Mater. Res. Soc. Symp. Proc
.,
1999
), Vol.
583
, p.
267
.
54.
J. M.
Ziman
,
Electrons and Phonons: The Theory of Transport Phenomena in Solids
(
Oxford University Press
,
London
,
1960
).
55.
J.
Duda
,
T.
English
,
D.
Jordan
,
P.
Norris
, and
W.
Soffa
,
J. Heat Transfer
134
,
014501
(
2011
).
56.
J. C.
Duda
,
T. S.
English
,
D. A.
Jordan
,
P. M.
Norris
, and
W. A.
Soffa
,
J. Phys.: Condens. Matter
23
,
205401
(
2011
).
57.
C. H.
Baker
and
P. M.
Norris
,
Phys. Rev. B
91
,
180302
(
2015
).
58.
D.
Frenkel
and
B.
Smit
,
Understanding Molecular Simulation: From Algorithms to Applications
, Computational Science Series (
Elsevier Science
,
2001
).
59.
S.
Nosé
,
J. Chem. Phys.
81
,
511
(
1984
).
60.
W. G.
Hoover
,
Phys. Rev. A
31
,
1695
(
1985
).
61.
C. Z.
Wang
,
C. T.
Chan
, and
K. M.
Ho
,
Phys. Rev. B
42
,
11276
(
1990
).
62.
J. R.
Lukes
and
H.
Zhong
,
Int. J. Heat Transfer
129
,
705
(
2007
).
63.
J. E.
Turney
,
A. J. H.
McGaughey
, and
C. H.
Amon
,
Phys. Rev. B
79
,
224305
(
2009
).
64.
K.
Albe
,
K.
Nordlund
,
J.
Nord
, and
A.
Kuronen
,
Phys. Rev. B
66
,
035205
(
2002
).
65.
G.
Srivastava
,
The Physics of Phonons
(
Taylor & Francis
,
London
,
1990
).