In this work, we have determined that there is a major flaw in the way phonons in random alloys are conceptually understood and consequently the way their transport of heat has been modeled. The basic problem is that phonons are almost ubiquitously considered plane waves and this basic concept has motivated modeling their transport with the phonon gas model (PGM), whereby phonons are treated like particles with a certain velocity that scatter with each other, boundaries, interfaces and imperfections. However, what we’ve determined in this study is that phonons are only plane waves in the limit of a perfectly pure defect-free crystal. As a result, it is only in this limiting case that the PGM becomes valid. What we’ve discovered is that even at low impurity concentrations ~ 1%, surprisingly phonons cease to be pure plane waves (Figure 1a) as they quickly lose their propagating character and can no longer be modeled properly with the PGM. As seen in Figure 1a, the degree of periodicity of vibrational modes, i.e., Eigenvector Periodicity Parameter (EPP) changes significantly even at low alloy concentrations and by increasing the alloy concentration, the character of propagating vibrational modes changes and they approach to non-propagating character (EPP<0.2).

We’ve proven the change in vibrational mode character by showing that modeling the thermal conductivity (TC) of a crystalline random alloy, with the most state of the art approach, namely the virtual crystal approximation (VCA), leads to major errors for thin films (Figure 3). The VCA exhibits large errors in predicting the thermal conductivity of an alloy thin film (See experimental section for details), because it assumes all phonons are like plane waves/particles that travel and scatter from impurities. The problem is that this assumption is not true, yet it forms the foundation of what a phonon is believed to be. Given this issue, we’ve developed a new method that is rigorous and distinct from the PGM. Our new method does not make this questionable assumption and shows excellent agreement with all of the experimental data, but our new method requires us to rethink how phonons behave entirely.

The consequences of this study are far reaching, and impact a wide variety of areas of science and engineering because alloys and materials with less than 99.9% purity, as well as materials with defect concentrations > 1% and fully amorphous materials are used throughout a wide variety of industries involving heat transfer, ranging from solar cells and superconductors, to LEDs and polymers in food packaging. Essentially, what this study concludes is that:

  1. Phonons cannot, in general, be thought of as plane waves or particles that travel and scatter.
  2. In general, there are three different types of phonons, namely propagons, diffusons and locons, and each has its own distinct physics for how it interacts with other quantum particles (i.e., phonons, electrons, photons, neutrons, etc.).
  3. In a material with ~5% disorder or more, the momentum of phonons cannot be assumed as hk . Therefore, the physics of how phonons in disordered materials contribute to phenomena such as superconductivity should be re-examined and possibly re-conceptualized, given these new results.

Figure 1. eigenvector periodicity parameter for InxGa1-xAs. (a) Eigenvector periodicity parameter for different alloy compositions, (b) average eigenvector periodicity parameter and fraction of propagating modes Indium composition.

Figure 2| TC of In0.53Ga0.47As 1.6 micron film. Temperature dependent TC of In0.53Ga0.47As film and the corresponding theoretical predictions using VCA and GKMA. The error bars were determined based on the standard deviation of GK results at a given temperature. Each labeled curve highlights the respective contributions associated with propagons, diffusons and locons, according to the GKMA and EP methodologies.

Mode character of phonons changes in alloys: In0.53Ga0.47As Example

To see how the mode character changes when EPP decreases, the eigenvectors of six randomly selected vibrational modes with different values are shown in Figure 4.3. Such visual inspection of modes is insightful, because it shows very clearly what aspects of a mode’s propagating character are first lost and how the transition from propagon to diffuson occurs. The panels on the left show the eigenvectors along the length of the computational domain while the panels on the rights show the eigenvectors viewed from the cross section of the computational domain. It can be seen that for EPP=1, the eigenvectors correspond to plane wave modulated vibrations (i.e., propagating waves). In this case, the wavelength can be clearly recognized as the distance over which the wave’s shape repeats. By decreasing the EPP, the vibration of atoms becomes more random and the modes tend to become non-propagating. For example, when EPP=0.6, one can clearly see more random vibrations in the cross-section of the computational domain compared to the mode with EPP=0.8.  However, in both cases, along the length of the computational domain, periodicity is still evident. For EPP=0.4, the eigenvector seems to be completely random when viewed from the cross-section, and one might consider this mode as diffuson. However, this mode can still be considered a propagons, because the transition between propagon and diffuson occur at EPP=0.2. For EPP=0.2, the vibrational directions are fully random and the mode is a diffuson, hence one cannot define an effective wave vector for this type of vibration.

Sample fabrication

The In0.53Ga0.47As epitaxial layers were grown on n-type (100) InP substrates by metalorganic chemical vapor deposition (MOCVD) in a Thomas Swan 7×2’’ reactor system equipped with a close-coupled showerhead growth chamber. Epipure™ trimethylindium (TMIn, In(CH3)3) and triethylgallium (TEGa, Ga(C2H5)3) were used as group III precursors while arsine (AsH3) and phosphine (PH3) were used as group V precursors. The growth of InP buffer layer on the top of InP substrate was performed at temperature of Tg=600 °C under a reactor chamber pressure of 100 Torr, and then the temperature was elevated to Tg=650 °C for a growth of In0.53Ga0.47As layer. For material characterizations, the Nomarski optical microscopy and X-ray diffraction (XRD) were utilized to analyze the morphological and structural properties of the epitaxial layers, respectively. The thickness of the In0.53Ga0.47As epitaxial layers in this study was 120 nm and 280 nm as determined by scanning electron microscope (SEM) in Fig. 3.

Figure 3. SEM images of InGaAs sample cross-sections used to determine the alloy layer thicknesses as well as the aluminum transducer.

Thermal Conductivity Measurements

Measurements of thermal conductivity were performed using time domain thermoreflectance (TDTR).While specific TDTR experimental setups (Figure 4) may vary, the basic premise of the technique remains the same.  Typically, the method will employ an ultra-fast laser with a pulse width that is less than 1 ps.  This pulse is split using polarizing optics, and is directed onto two separate paths. A pump path is used to create a periodic heating event on the surface of a sample, and a probe path is used to monitor the temperature decay at the surface via a change in reflectivity. The experimental data are then fitted to a diffusive heat conduction model over timescales on the order of several nanoseconds. Because of the ultra-short laser pulse and high oscillation frequency, TDTR enables accurate measurements of thermal properties of thin films on the order of 100’s of nanometers or less, depending on the specific sample configurations being measured.

Figure 4: Schematic of TDTR setup

Thermal conductivity of In1-xGaxAs alloy using Virtual Crystal Approximation

The Virtual Crystal Approximation was used to take into account the alloy effects in thermal conductivity. In this approach the disordered crystal is replaced with an ordered one with compositionally weighted IFCs, atomic mass and lattice constant according to composition. The mass disorder and anharmonicity are both treated as a perturbation. The net scattering rate of a phonon mode is calculated as the sum of scattering rate due to mass disorder and anharmonicity, according to Matthiessen’s rule:

The first term is the phonon-phonon scattering rate which is calculated in the same way as the IHPCs except that for the virtual crystal, the compositionally weighted mass and lattice constant are used. The second term is harmonic mass disorder scattering rate, which is calculated by using perturbation theory

Finally the size effect is applied using Matthiessen’s rule

In this study, we also employed a second method for calculating the thermal conductivity of the alloy using VCA, whereby fitting parameters where used to describe phonon-phonon and phonon-alloy impurity relaxation times originally. In this approach the thermal conductivity is given by

Where kB is Boltzmann’s constant,  T is temperature, h is Planck’s constant divided by, y is a dimensionless parameter, omega_c is the cut-off frequency and obtained using Debye model. The scattering time for a given frequency is related to individual processes via Mattheissen’s rule

Where first, second, and the last terms are the umklapp, mass disorder, and boundary scattering times, respectively. these are given by

Where L is film thickness and constants A, B, and C are fitting parameters.

Thermal conductivity of In1-xGaxAs alloy using Green Kubo Modal Analysis (GKMA)

In this method, first, the harmonic frequencies and eigenvectors are obtained from a Lattice Dynamics (LD) calculation of the entire atomic supercell. To obtain the modal contributions to the velocity of each atom (e.g., x_i(n,t), atom i , mode n) the atom velocities from molecular dynamics were projected onto the normal mode basis. Then the modal heat flux is calculated by substituting the modal velocity into heat flux operator,

Where Ei is the sum of potential and kinetic energy of atom i , Phi_i  is the potential energy of atom j, V is the volume of the supercell, and rij is the distance between atom i and atom j. Finally, the thermal conductivity of each vibrational mode can be calculated by substituting the modal heat flux into the Green-Kubo expression,

Due to the classical nature of molecular dynamics (MD) simulation, the low-temperature calculations (below Debye Temperature) of thermal conductivity become inaccurate. This is due to the fact that MD doesn’t reproduce the proper mode amplitudes that correspond to quantum occupations. Therefore, we applied quantum heat capacity correction factor.