How to resolve a phonon-associated property into contributions of basic phonon modes

Many properties of materials are associated with phonons. To better understand the phonon-property relation, it is a common practice to decompose the phonon-associated property into the contributions of basic phonon/vibration modes (e.g. longitudinal/transverse acoustic/optical mode), and identify the mode(s) that dominate(s) the property. The existing methods rely on labelling the phonon into one of the basic modes (BMs), however, the vibration characteristics of many phonons are different from the definitions of BMs, indicating these methods may give wrong decomposition results. Here we present a new method based on treating the phonon as a mixture of the BMs. By aligning the wave vectors and then projecting the phonon eigenvector onto the eigenvectors of the BMs, we can obtain the weights of the BMs on the given phonon, which can be used to quantify the contribution of each BM to the property. As an example, we apply this method to unravel the phonon modes that dominate the scattering of the electrons at the conduction band edge in two-dimensional antimony, an emerging semiconductor that has attracted great interest for electronics, but its mobility-limiting factors remain unclear. We find that the electron scattering is dominated by the out-of-plane acoustic phonon mode followed by the longitudinal acoustic mode, which are different from the results of other methods. Our method is generally applicable to different kinds of phonon-related properties and to all crystal materials.


Introduction
Phonons play an important role in many properties of materials, such as thermal conductivity [1][2][3][4][5][6][7], electron/ hole mobility [8][9][10][11][12], thermoelectric efficiency [13,14], superconductivity [15][16][17][18], light adsorption/emission [19][20][21][22], and electron-hole recombination [23,24]. Resolving the phonon-associated property into the contributions of basic phonon/vibration modes is a very useful approach to understand what basic phonon mode(s) govern(s) the property [4, 10-12, 17, 25-27]. For a given phonon, its mode can be identified by comparing the vibration pattern/characteristics (quantified by the eigenvector, e) with the direction of the wave vector (q). For every crystal material, there exist d×N (where d is the vibration degree of freedom, and N is the number of atoms in the unit cell) basic modes (BMs), such as 'longitudinal acoustic' (LA) mode, and 'transverse acoustic' (TA) mode. For example, if all the atoms in the unit cell move along the same direction parallel to q, then it is called LA mode. If they all move along the same direction but perpendicular to the q, then it is defined as TA mode. However, except those with the q located at/near high-symmetry positions of the Brillouin Zone (BZ), most phonons do not possess such well-defined vibration characteristics with respect to the wave vector direction, and thus cannot be classified into any of the BMs. For example, in the two-dimensional (2D) square lattice with only one atom in the unit cell, there are 2×1=2 BMs: LA and TA; only the phonons located at the high symmetry lines exhibit LA/TA modes, while for most phonons, the atomic vibration is neither parallel nor perpendicular to the q, hence it is problematic to categorize them into LA/TA mode (this will be discussed later in detail).
Although the vibration pattern with respect to the q direction for a given phonon is often different from that of the BMs, the existing methods still classify/label it into one of the BMs, and attribute the associated property into the full contribution of that BM. There are two methods commonly used in literature: (1) energy method: one first sorts the BMs according to the energies of long wavelength phonons near the Γ point. Then for a given phonon, one compares its energy with the other phonons with the same q. If it is ranked nth in energy, then it is classified into the nth BM. (2) continuity method [28,29]. For a given phonon, one compares its eigenvector (e, which represents the vibration pattern) with the eigenvectors of the phonons that are close in q, and classify it into the BM of the phonon with the most similar e. This comparison starts from the Γ point, and spans over the BZ so that every phonon is assigned to one of the BMs. Both methods rely on classifying the phonon into one of the BMs, which is valid only for phonons located at/near high-symmetry positions of the BZ.
Here we present a new method that does not have these problems. In this method, the phonon is not classified into one of the BMs; instead, it is treated as a mixture of the BMs. By rotating and projecting its e onto the e of BMs, we can obtain the weights of the BMs on the given phonon and thus their contributions to the property. This method is thus called the 'projection method'. This method has been implemented into an opensource code that can be found at the corresponding author's website. We then apply this new method to unravel the phonon modes that dominate the scattering of the electrons at the conduction band edge in 2D antimony (Sb), which is a promising material for electronics, but its mobility-governing factors remain unclear. We find that the electron scattering is dominated by the out-of-plane acoustic phonon mode followed by the longitudinal acoustic mode, different from the results of other methods.

Results and discussion
Our method is composed of three steps (figure 1(a)): (1) Choose reference phonons and identify the BMs. Since the phonons near the Γ point have well-defined vibration patterns with respect to the wave vector direction following the definition of BMs, we thus choose these phonons to identify the BMs. Specifically, we select an arbitrary q (e.g. q r ) near the Γ point (i.e. |q r |∼0), and obtain the e for the corresponding 3 N phonons (e r,i , where i is the band index, and=1, 2, K 3N) . We call these phonons as reference phonons, q r as reference wave vector, and e r,i as reference eigenvectors. Each reference phonon represents one BM, which can be identified by the difference between e r,i and q r direction. Note that the choice of q r would not affect the basic-mode-resolved results as long as |q r |∼0 (in examples shown later, we use |q r |=0.005|G|, where G is the unit cell vector of the reciprocal lattice; we have also tested |q r |=0.01|G| and found that the decomposition results change by <0.01%). (2) Rotate vectors. For a given phonon with wave vector q g and eigenvector e g,j (where j is the band index at q g ), in order to quantify the contribution of the ith BM, one needs to compare the difference between the e g,j and the q g direction, with that between e r,i and the q r direction. Note that the vibration mode is defined by the relation between the e and the direction of q, rather than the e or the q alone; thus one cannot directly compare e g,j with e r,i . Since the direction of q g is often different from that of q r , to help the comparison, we rotate both q g and e g,j so that the q g is aligned with the q r . The e g,j has to be rotated simultaneously to retain the relation between the e g,j and the q g . More specifically, to align the q g and q r , we can rotate the q g through an axis that passes through the origin of the BZ and is perpendicular to both q r and q g . The axis can be characterized by a unit vector u: Note that the rotation should be counterclockwise when viewed along the -u direction, and the rotation angle θ should be: The u and θ then define a rotation matrix R(u, θ) that is applied to rotate the e g,j to be e′ g,j : Schematic illustration of the projection method. q r , e r,1 and e r,2 are the wave vector and eigenvectors of the reference phonons (|q r |∼0). q g and e g are the wave vector and eigenvector of an arbitrary phonon. θ is the angle between q g and q r . The e g is rotated and projected on to e r,1 and e r,2 , giving p 1 and p 2 as the coefficients that can be used to decompose the phonon-related property. In Cartesian reference frame, the rotation matrix is [30]: cos  1  cos  1  cos  sin  1  cos  sin   1  cos  sin  cos  1  cos  1  cos  sin   1  cos  sin  1  cos  sin  cos  1  cos , 4 where the (u x , u y , u z ) are the coordinates of the u. (3) Project the rotated eigenvectors. To quantify the contribution of each BM, the e′ g,j is then projected onto the reference eigenvectors, giving: where * denotes the complex conjugate. The e′ g,j thus can be written as: where i runs from 1 to d×N. Since the phonon eigenvectors at the same q are orthonormal to each other, we have: Therefore, the p ij * p ij can be thought as the weight of the ith BM on the given phonon (with wave vector q g and band index j), which can be further denoted as w i : For a property related to this phonon, we can then use p ij * p ij to characterize the contribution of the ith BM to that property. We can also express the e r,i using the e′ g,j : p pp e e ' , and 1. 9 The physics behind our method is that any arbitrary atomic vibration can be considered to be a superposition of the elementary vibrations, and the mathematics underlying this method is that any arbitrary vector can be expressed by the linear combination of a set of basis vectors. Note that if the phonon-associated property is zero, such as the electron-phonon coupling strength in zigzag graphene nanoribbons with even dimers (due to the selection rule of the parity conservation) [31], then the contribution of each BM is also zero. Similar projection procedures have been applied for other purposes, such as unfolding the phonon dispersions of a supercell into the 1st BZ of the primitive cell [32][33][34][35][36][37].
As the first example, we study a model system-the 2D square lattice with one atom per cell ( figure 1(b)). The interactions up to second-nearest neighbors are considered. By solving phonon secular equation, we can obtain the energies and eigenvectors for any phonons (see the supplementary material available online at stacks.iop. org/JPMATER/2/045005/mmedia). Figure 1(c) shows the phonon dispersion along the Γ-M-X-Γ line, which has two branches. The phonons near the Γ point have the atomic vibration either parallel or perpendicular to the q, which define two BMs: LA mode and TA mode. The phonons located at the Γ-M and X-Γ lines also have such vibration patterns (left and right insets in figure 1(c)), and thus can be classified into LA/TA mode. However, for phonons located at the M-X (and other regions), the vibration is neither parallel nor perpendicular to q (central insets in figure 1(c)), hence, they cannot be categorized into any of the BMs. Using our method, we calculate the weights of the BMs on a given phonon, and plot them in figures 1(d)-(e). At the Γ-M line, the higher-energy branch is fully composed of LA modes and the lower-energy branch is fully composed of TA modes. The M is a singular point, where the two eigenvectors can have arbitrary orientations as long as they are normal to each other (see the supplemental material). The phonons approaching the M from the X-M direction have an equal weight of LA and TA modes. From M to X, the higher-energy branch becomes more LA-dominated and the lower-energy branch becomes more TA-dominated. At the X-Γ line, the phonons are again composed of only LA/TA mode. These weights can also be visualized through color maps, as shown in figures 1(f)-(g). The LA/TA mode is localized onto the higher/lower-energy branch at the Γ-M and X-Γ lines, while it spans to the lower/ higher-energy branch at the M-X line. For comparison, we also show the phonon dispersion colored by the conventional methods (figure S1(a) for the energy method and figure S1(b) for continuity method). Clearly, neither of them captures the real vibration characteristics at the M-X line.
We then employ our method to study a real material: 2D Sb (β-antimonene) [38][39][40][41]. It has a buckled hexagonal structure with 2 atoms per cell, as shown in figure 2. The phonon dispersion and eigenvectors are calculated from first principles (see Methods section). By examining the phonon vibration characteristics with respect to q near the Γ point, we identify 6 BMs: out-of-plane acoustic (ZA), TA, LA, transverse optical (TO), longitudinal optical (LO), and out-of-plane optical (ZO) modes. Then for a given phonon, we can rotate its q and e, and project the rotated e onto the e of the reference phonons, as described earlier. These steps give the weights of the BMs on the given phonon, which are visualized by the color maps in figure 2. These figures show that the phonons near the Γ point comprise of only one BM (as expected), however, most of the other phonons, especially those around the BZ boundaries and corners, are a mixture of different BMs. In other words, the BMs are localized to a single branch near the Γ point, while spread over multiple branches in other regions. It is interesting to find that although there is a phonon energy gap (∼10 meV) between the low energy branches and high energy branches, it does not prevent the mixture of the BMs. In contrast, the conventional methods classify every phonon into one of the BMs (figure S2). It is worth noting that in the continuity method, the TA phonons approaching the Γ point from the K-Γ line have higher-energy than LA phonons, which contradicts with their order at the Γ-M line. This is incorrect because near the Γ point, the energy order between TA and LA phonons with the same wave vector should be fixed, regardless of the direction of the wave vector. This issue may be because the path crosses the K point where the e can have multiple choices.
The weights obtained from our method can be used to decompose the phonon-relevant properties. Here we use the scattering rate of electrons in 2D Sb as an example. There is great interest in finding high-mobility 2D semiconductors for next-generation electronics, and 2D Sb has been experimentally demonstrated to be a promising candidate, while its mobility-governing factors remain unclear. The carrier mobility is affected by many scattering mechanisms, and the scattering by phonons defines the intrinsic mobility. The EPC in monolayer Sb has been studied previously [42,43]. Pizzi et al used the deformation potential theory to estimate the electron coupling with the long wavelength LA phonon but neglected the coupling with other phonons [43]; Lugovskoi et al used the perturbation theory and the Wannier interpolation to study the EPC in highly doped 2D Sb and the superconductivity [42], where the BM contributions to the EPC strength are obtained by the energy To answer these questions, here we focus on the scattering of the electrons at the conduction band minimum (CBM), as this state is most populated when Fermi level is within the band gap, and the scattering mechanisms that dominate this state are likely to be also important for other states. The relaxation time (τ) or the scattering rate (1/τ) for the CBM can be calculated as [16]: where m and k, l and k′, j and q, are the band index and wave vector for the initial electronic state, the final electronic state, and the absorbed/emitted phonon, respectively. 1/t(lk′, jq) is the transition rate (here the m and k are omitted because the initial state is fixed to be CBM). ω is the phonon frequency, v is the electron velocity, f is the Fermi distribution, n is the Bose distribution, G is the reciprocal lattice vector, g is the electron-phonon coupling (EPC) strength, and Ω BZ is the volume/area of the first BZ. The Fermi level is placed at the middle of the band gap, representing an intrinsic/undoped semiconductor. Note that here we assume the phonons and the electrons are not 'hot', i.e. their population follows equilibrium distribution. This assumption is valid for studying charge carrier mobility, while invalid for non-equilibrium systems with hot phonons/electrons [44].
In this example, the property associated with the phonon jq is the transition rate 1/t(lk′,jq). The contribution of the ith BM thus can be calculated as: where w i is the weight of the ith BM on that phonon, which can be obtained from our projection method. By substituting the 1/t with 1/t i in equation (10), we can get the ith BM-resolved scattering rate, 1/τ i . All the relevant quantities are calculated from first principles, as detailed in the Methods section. To identify what BM gives the strongest scattering, we compare their contribution to the total scattering rate in figure 3. It is found that the ZA mode scatters most strongly (contributing to 37.6% of the total scattering rate), followed by the LA mode (21.8%). These are different from the results using the conventional methods: the continuity method gives the TO (23.5%), TA (23.3%) and ZA (21.3%) as the main scattering modes, and the energy method results in ZA (42.5%) and TO (22.9%). In addition to Sb, we have also studied the electron scattering in 2D As, and found similar results (see figure S7 and table S2 in the SM). Note that the projection method and the energy method give the same leading scattering mode (ZA, although with different weights: 37.6% versus 42.5%). This may be because in this example, the scattering is dominated by phonons near Γ point, where the BMs are localized onto individual branches. It is possible that in other systems the property is mainly governed by the phonons far from the Γ point, then a larger difference between the results obtained using the projection method and the energy method can be anticipated, which deserves further separate studies. The dominance of ZA-phonon induced scattering can be explained by the parabolic dispersion of the ZA mode in ultrathin material and the absence of the mirror symmetry in the atomic structure of 2D Sb, which give rise to a large density of ZA phonons that can scatter the electrons [45]. If the ZA and/or LA phonons can be frozen, which may be achieved using adhesive substrate and/or strain, then the electron mobility in 2D Sb can be improved. Note that different decomposition methods do not change the total value of the property. Our method can also be applied to decompose other phonon-associated properties, such as the phonon-phonon scattering rate [46], which is important for thermal transport. In this case, the procedure should be the same as that for the electron-phonon scattering, except that the equation (10) should be replaced by that for phonon-phonon scattering [47]. In summary, we have presented a new method to resolve the phonon-relevant properties onto the contributions of the BMs. Different from conventional methods that classify the phonon into one of the BMs, we treat the phonon as their linear combination. By rotating and projecting the phonon eigenvector onto the eigenvectors of the BMs, we can obtain their weights on the given phonon and thus their contributions to the property. This method can help reveal/clarify the phonon modes that dominate the relevant properties of materials. For example, we apply this method to unravel the phonon modes that dominate the electron scattering in two-dimensional antimony, and find that it is the ZA phonons that contribute most strongly, followed by the LA phonons. Our method is generally applicable to different kinds of phonon-related properties [46], and to all crystal materials, and thus advances the studies of systems where phonons are important.

Methods
The analytic derivation of the phonon dispersion of 2D square lattice can be found in the supplemental material. The phonon dispersion and eigenvectors of 2D Sb are calculated using the Quantum Espresso Package [48,49] with the norm-conserving pseudopotentials and the Perdew-Burke-Ernzerhof [50] exchange-correlation functional. The kinetic energy cutoffs for wave functions and charge density are set to 60 and 240 Ry, respectively. The structure is fully relaxed with a 36×36 k mesh [51] until the magnitude of the force acting on each atom becomes less than 0.0001 Ry Bohr −1 . To obtain the scattering rate, we first use the density functional perturbation theory [52] to calculate the EPC strengths on a coarse grid of 9×9 k and 9×9 q points. Then the EPC strengths on a fine grid of 1000×1000 q points are interpolated using the Wannier interpolation [53] as implemented in the EPW code [54].