Consider the lattice dynamics (LD) formalism whereby the vibrations of any stable system can be decomposed into a series of eigen mode solutions. In the purely harmonic limit, the solutions to the equations of motion are exact and the resulting eigen modes do not interact. Anharmonic systems can be analyzed by using the harmonic modes as a basis set, however, the anharmonicity causes the amplitudes to become time dependent and they cannot be determined analytically. From this perspective, we can still view the vibrations as arising from the modes obtained in the harmonic limit, but the interactions between modes manifest through the time dependent mode amplitudes, which can be studied through a fully anharmonic numerical simulation technique such as molecular dynamics (MD). It is in this sense that the positions and velocities contain all degrees of anharmonicity, and the harmonic eigen modes simply serve as the basis set from which we can interpret/understand them. Using the LD formalism, we can then use the eigen mode coordinates to write the velocity of every individual atom as,
where denotes the atom, is its mass, is the number of unit cells in the system, and is the polarization vector which gives the direction in which each atom moves. Equations 1 and 2 represent a transformation to and from the normal mode coordinates. Equation 1 gives the modal contributions to every individual atom’s velocity. From Eqs. 1 and 2 one can describe the modal contributions to any quantity that is a direct function of the atomic displacements or velocities (see supplementary materials for details). To study thermal conductivity, we substituted the modal velocity into the heat flux operator , to obtain the modal heat flux
We can then take this expression for the heat flux and substitute it directly into the GK expression for modal thermal conductivity,
This new expression has several useful features. For one, with this method, we are able to calculate each mode’s contribution to thermal conductivity.
The above formula is useful for analyzing wave-vector very defined modes, like phonons in the crystalline solids. However GKMA could be changes slightly to study the amorphous materials from gamma-point lattice dynamics results. Here, we utilize the same concept of projecting the MD trajectory onto the modes obtained from LD, but instead of using propagating phonons with a well-defined wave vector, we conduct a direct modal decomposition of the heat current using eigen vectors for all phonons at gamma point. Towards this end, we changed the Eqns 1-5 accordingly, where the normal mode amplitudes are calculated from,
where the displacement and velocity of each atom and can be obtained from the normal mode coordinates via,
In Eqs. (6-9) n denotes the mode (e.g., the nth solution to the equations of motion), mj is the mass of the jth atom, and is the polarization vector which gives the magnitude and direction of motion for atom j in mode n. Equations (6-9) are not new and are well established within the context of the LD formalism.
The critical enabling insight offered herein is the physical meaning associated with Eqs. (8) and (9). Here, we postulated that Eqs. (8) and (9) essentially state that at every instant, each atom’s position and velocity is composed of their respective contributions from the different collective vibrations in the system. Thus, every atom’s position and velocity are dictated by the respective magnitudes of each normal mode’s amplitude and its time derivative . By thinking of each atom’s position and velocity as being composed of an exact sum of modal contributions at every instant, it was then postulated that if an individual mode’s contribution to the displacement or velocity of an atom is used in an expression for the calculation of any other property that depends on that atom’s position and/or velocity, one would subsequently obtain that mode’s contribution to that property. For example, one could calculate each mode’s contribution to the temperature of the system as discussed in the supplementary materials. Similarly, towards the calculation of thermal conductivity, the modal contributions to the velocity of each atom can be substituted into the heat flux operator derived by Hardy, to obtain each mode’s contribution to the volume averaged heat flux at each time step in a EMD simulation,
where , is the sum of potential and kinetic energy of atom i, V is the volume of the supercell, is the potential energy of atom j, and is the distance between atom i and atom j.
We can then take Eq. 10 and substitute it directly into the Green-Kubo expression for thermal conductivity, to obtain an equation that expresses the thermal conductivity as a direct summation over individual mode contributions,
Furthermore, one can also substitute the summation of modal contributions to the heat flux in both places of the heat flux autocorrelation to obtain the thermal conductivity as a double summation over individual mode-mode heat flux cross-correlation functions,
Equations (11) and (12) are the primary results termed the GKMA method, which now allows us to obtain each mode’s contribution to the total thermal conductivity, and we are guaranteed by Eq. (9) that the summation in Eqs (11) and (12) will exactly recover the total GK thermal conductivity. Also, Eq. (12) allows one to examine how the correlation between pairs of modes contributes to thermal conductivity.