Complementary local-global approach for phonon mode connectivities

Sorting and assigning phonon branches (e.g., longitudinal acoustic) of phonon modes is important for characterizing the phonon bands of a crystal and the determination of phonon properties such as the Gr\"uneisan parameter and group velocity. To do this, the phonon band indices (including the longitudinal and transverse acoustic) have to be assigned correctly to all phonon modes across a path or paths in the Brillouin zone. As our solution to this challenging problem, we propose a computationally efficient and robust two-stage hybrid method that combines two approaches with their own merits. The first is the perturbative approach in which we connect the modes using degenerate perturbation theory. In the second approach, we use numerical fitting based on least-squares fits to circumvent local connectivity errors at or near exact degenerate modes. The method can be easily generalized to other condensed matter problems involving Hermitian matrix operators such as electronic bands in tight-binding Hamiltonians or in a standard density-functional calculation, and photonic bands in photonic crystals.


I. INTRODUCTION
Lattice dynamical studies [1,2] are indispensable for a detailed understanding of the thermodynamics, phase stabilities, phase transitions, and thermal properties of materials [3][4][5][6]. The methods used to calculate the properties related to phonons fall into two general categories with their respective advantages. The methods in the first category are versatile since they apply density-functional perturbation theory (DFPT) [7,8] to a unit cell to compute analytic energy derivatives [9,10]. The methods in the second category use small atomic displacements [11][12][13][14][15][16] to numerically compute the interatomic force constants from the induced forces based on the Hellman-Feynman theorem [17]. Common to these two categories of methods is the problem of dealing with small changes in the dynamical matrices for the evaluation of phonon frequency derivatives with respect to strain or volume for the Grüneisen parameters [6,18,19] or with respect to wavevectors (or q vectors) for phonon group velocities. In this case, the perturbation method from quantum mechanics is extremely useful [20][21][22].
However the problem of phonon mode connectivity is of a somewhat different nature than the calculation of frequency derivatives since the goal is to connect phonon modes (to form bands, analogous to the electronic counterpart [23]) with neighboring q points and the connectivities are to be extended throughout the whole path in the Brillouin zone (BZ). When two neighboring phonon modes at q and q + δq are connected, they are both assigned the same band index which groups them into the same phonon branch. To connect two neighboring phonon modes, we need the eigenvalues and eigenvectors of the phonon modes which can be obtained by solving the eigenvalue problem [24] D(q)u q = ω 2 q u q (1) * ganck@ihpc.a-star.edu.sg where D(q) is the dynamical matrix of dimension N ×N , and u q and ω 2 q are the associated eigenvector and eigenvalue, respectively. For a crystal, N = 3n c where n c is the number of atoms in the unit cell. When diagonalized, the Hermitian matrix D(q) in Eq. (1) yields N phonon modes with eigenvalues ω 2 n,q (where ω n,q denotes the mode frequencies) and the corresponding eigenvectors u n,q , for phonon band index n = 1, 2, · · · , N . In the phonon mode connectivity problem, the phonon band indices should be assigned so that the band index of each phonon mode corresponds to the physical phonon branch to which it belongs. For instance, a phonon mode for graphene with n = 1 should belong to the out-of-plane flexural acoustic branch (see Fig. 1).
In practice, the band index of a phonon mode is often assigned by ranking the frequencies of all phonon modes at the same q point such that ω 1,q ≤ ω 2,q ≤ · · · ≤ ω N,q . In a simple crystal, this assignment usually works when q is close to (but not at) the center of the BZ because the phonon modes sharing the same n and q values usually belong to the same physical phonon branch at small q. However, farther away from the BZ, this naive assignment of the band indices by phonon frequency ranking fails because the crossing of phonon branches changes the order of the eigenvalues, resulting in a breakdown of the association between phonon frequency ranking and the physical phonon branch. Around these crossing points, the use of phonon frequency ranking to connect the phonon modes at different q points results in illusory level repulsion [25], as can be seen in Fig. 1 which shows the graphene phonon dispersion curves along the Γ-K direction with incorrect phonon frequency sorting and correct physical phonon branch sorting, with the phonon frequencies ω n,q obtained using the Tersoff interatomic potential [26,27].
One common solution to this problem is to utilize the similarities of eigenvectors between neighboring q points in the BZ [28]. Suppose the band indices for the eigenmodes at q are assigned correctly with respect to the physical phonon branch. We can then assign the band indices for the eigenmodes of a neighboring point q + δq through the eigenvector similarity condition i.e., an eigenvector at q + δq is assigned the band index n if its overlap with an eigenvector at q with the band index n is the closer to unity than that with any other eigenvectors at q + δq. The accuracy of the band index assignment in Eq. (2) depends on the size of δq and may fail when δq is relatively large. The operational assumption behind Eq. (2) is that under a small enough shift q → q + δq, the frequency ω n,q and eigenvector u n,q change smoothly to ω n,q+δq and u n,q+δq , respectively. However, accurate assignment of the band index in Eq. 2 requires the use of a very dense q grid in the BZ, which is computationally inefficient. Another problem with this eigenvector similarity approach is that it may fail at the crossing points that correspond to near or exact degenerate modes. The reason is that an arbitrary linear combination of eigenvectors of a degenerate eigenvalue is also an eigenvector with the same eigenvalue (to within a numerical tolerance). Therefore it is difficult to correctly connect modes based on the eigenvector similarity condition alone. One may bypass this problem by choosing a different δq during computing runtime to avoid the problematic q points although this entails a very complicated implementation where the outcome may still be questionable. It is important to recognize that Eq. (2) is entirely reliant on the local eigenvector similarity of neighboring q points which, near or at degeneracy points, becomes susceptible to assignment errors that can propagate across the BZ. Hence, a more robust approach that circumvents local assignment errors is needed.
To solve this connectivity problem, we first discuss how the eigenvectors and eigenvalues of the dynamical matri-ces at neighboring q points can be related systematically within the framework of perturbation theory. Using this theoretical framework, we show how the overlap condition in Eq. (2) can be derived and generalized. Exploiting the perturbative approach for estimating eigenvalue, we propose a new two-stage hybrid method where in the first stage of the algorithm, an eigenvalue-perturbative approach is used to connect modes of the neighboring q points in the beginning of a path in the BZ. However, this perturbative approach is error-prone near or at degenerate points as it only uses the local information around the q point. Therefore in the second stage of the algorithm, we switch to a more 'global' approach that uses compatibility with the eigenvalues of the phonon modes connected in earlier stages as a criterion for assigning the band index. We call the second stage the fitting approach since it is based on the least-squares fits to the modes connected in the first stage. As we shall see, the hybrid method completely solves the phonon mode connectivity problem by using global information to circumvent local assignment errors.
Our paper is organized as follows. In Section II, we relate the eigenvectors and eigenvalues of the dynamical matrices at neighboring q points through the framework of perturbation theory, and show how the eigenvectors and eigenvalues can be connected locally in a systematic fashion. The origin of the local band index assignment error is explained. We then present a least-squares fit approach to complement the eigenvalue-perturbative approach used for local phonon mode connectivities. In Section III, we illustrate our two-stage hybrid method by presenting the results for three systems: graphene, Bi 2 Se 3 , and MoS 2 . Finally we state our conclusions in Section IV.

A. Local mode connectivity
To explicate the relationship between the eigenvectors and eigenvalues at neighboring q points more systematically, it is necessary to recall that at q, the N × N Hermitian dynamical matrix D(q) [24], and its associated eigenvalues ω 2 n,q and eigenvectors u n,q , are governed by the equation where n denotes the phonon band index for n = 1, 2, · · · , N and u n,q is the corresponding column eigenvector which satisfies the orthonormality condition u † mq u n,q=δmn . At q + δq, the corresponding dynamical matrix and its associated eigenvalues and eigenvectors are similarly given by for m = 1, 2, · · · , N . For orientational purposes, we assume that none of the eigenvalues in Eqs. (3) and (4) are degenerate, i.e., they do not form crossing points. Traditionally, we assume that u n,q and u m,q+δq belong to the same phonon branch if their overlap is close to unity, i.e., because under a small enough wavevector shift q → q + δq, the frequency ω n,q and eigenvector u n,q deform smoothly to ω n,q+δq and u n,q+δq , respectively, for most q points. The limitation of the eigenvector similarity condition in Eq. (5) is that its accuracy diminishes rapidly as |δq| increases.
To relate Eqs. (3) and (4) formally, we treat D(q) and D(q + δq) as the analogs of the unperturbed and perturbed Hamiltonian in perturbation theory, respectively. The analog of the perturbation term in this case is To facilitate an expansion in a dummy variable λ, we consider a general perturbation term λV so that u n,q+δq is expanded as a power series in λ, i.e., where the zeroth-order term (j = 0) in the series is u (0) n,q+δq = u n,q . Likewise, the eigenvalue ω 2 n,q+δq can be expanded as where ω (0) n,q+δq = ω n,q . The higher-order terms in the power series in Eq. (7) can be obtained using the familiar nondegenerate perturbation theory in quantum mechanics. For instance, the first-order in Eq. (7) The last equality follows because u † mq D(q)u n,q = 0 for m = n. Therefore, to the first order, we obtain u n,q+δq = u n,q + λu Similarly, the first-order term in Eq. (8) is If we set λ = 1, then we can define the estimated eigenvector with band index n at q + δq as u n,q+δq = u n,q + u and its associated eigenfrequency as There are two ways that perturbation theory can be used for phonon mode connectivity given Eqs. (12) and (13). The first way is to use Eq. (12) to sort the eigenmodes at q + δq with respect to the eigenmodes at q using a generalized version of Eq. (2). Suppose we have a set of eigenmodes at q with the correct band indices (n = 1, · · · , N ) and their eigenvectors given by u n,q . If we have an eigenmode with eigenvector u m,q+δq satisfying Eq. 4, then we can assign it the band index n if it satisfies the overlap condition where u n,q+δq is as defined in Eq. (12). If we use only the first term for u n,q+δq on the right hand side of Eq. (12) and substitute it to Eq. 14, we recover Eq. (2). The accuracy of the overlap condition of Eq. (14) may be increased by using a higher-order estimate for Eq. (12) at the cost of increasing complexity in implementation.
In the second way, we connect the eigenmodes at q and q + δq by comparing the estimated eigenvalues Ω 2 n,q+δq from Eq. (13) to the exact eigenvalues ω 2 m,q+δq from Eq. (4). We assigned the band index of n to ω 2 m,q+δq (i.e., m = n) if the difference between Ω 2 n,q+δq and ω 2 m,q+δq is minimum.
We now provide more details on the implementation of the second way mentioned above. [29] From two chosen points q b and q e that specify a path along a high symmetry direction in the BZ, we may define (p + 1) q points to uniformly connect q b to q e . The (p + 1) dynamical matrices D(q i ) at q i are evaluated and diagonalized to obtain (p + 1) lists of increasing eigenvalues ω 2 n,q i , where n = 1, 2, · · · , N . Due to possible band crossings, we cannot simply join the values of ω n,q i for a particular n value, for i = 0, 1, · · · , p and identify them as the nth phonon band.
To connect the frequencies at q i to the frequencies at q i+1 (obviously i = 0 for a fresh start), we use the following method. At q i , we have the phonon frequencies ω n,q i and their associated eigenvectors u n,q i , n = 1, 2, · · · , N . From these eigenvalues and eigenvectors, we find the perturbed (or predicted) eigenvalues Ω 2 n,q i+1 at q i+1 , correct up to the first order, through Ω 2 n,q i+1 = ω 2 n,q i + ∆ n,q i where the correction terms ∆ n,q i are to be deduced from the first-order perturbation theory with the perturbation term V i = D(q i+1 ) − D(q i ) as in Eq. (6).
However, due to the presence of possible degenerate eigenvalues at q i , we must apply degenerate first-order perturbation theory to find the correction terms ∆ n,q i . Specifically this is achieved by first partitioning all N Figure 2. A schematic to demonstrate how the perturbed eigenvalues (crosses) obtained from the perturbation theory are used to establish two mode connectivities of actual eigenvalues (circles) from q i to q i+1 . The actual eigenvalues are deduced from the diagonalizations of dynamical matrices at q i and q i+1 . The scenarios of (a) noncrossing and (b) crossing of bands show the importance of the accuracy of the predicted eigenvalues for establishing the correct mode connectivities.
eigenvalues ω 2 n,q i into a number of clusters of eigenvalues, where each cluster consists of d eigenvalues (where in general two clusters may have different d) such that the frequency of any mode in the cluster is close to the frequency of at least one other mode in the cluster by a tolerance τ . Hence, each cluster represents a numerically degenerate subspace for which the first-order correction has to be computed separately. We obtain for each cluster a d × d matrix A with matrix elements A αβ = u † α,q i V i u β,q i . A diagonalization of A gives d first-order corrections ∆ n,q i for d eigenvalues ω 2 n,q i in the cluster. When a mode is singly degenerate (or nondegenerate, i.e., d = 1), a cluster consists of just one distinct eigenvalue and we recover the result of the standard nondegenerate first-order perturbation theory in Eq. (11).
From the ordering deduced from the sorted perturbed eigenvalues Ω 2 n,q i+1 and the ordering of the sorted actual eigenvalues ω 2 m,q i+1 which we obtain from diagonalizing the dynamical matrix at q i+1 , we establish a very reasonable mode connectivity from q i to q i+1 by simply connecting the ω 2 n,q i and ω 2 m,q i+1 guided by Ω 2 n,q i+1 , as can be seen in Fig. 2. When the connectivities between q i and q i+1 , we can use the same perturbative approach to establish connectivities between q i+1 and q i+2 , etc.
We note that the perturbative approach for local connectivity may fail at q i due to the presence of degenerate or near-degenerate modes for the following reason. Suppose the correct scenario of the connectivities of the two modes is to cross as shown in Fig. 2(b). This switching has to be facilitated by the switching of the order of two predicted eigenvalues (as shown on the left of Fig. 2(b)). However, if these two predicted eigenvalues fail to switch due to reasons such as inaccuracies in the dynamical matrices, too large a V i term, or the predicted eigenvalues lying extremely close to one another as in the case of degenerate modes where it leads to an incorrect ordering of predicted eigenvalues as shown on the left of Fig. 2(a), then it will lead to an incorrect noncrossing scenario. This means that the perturbed eigenvalues must be correctly ordered if they are to give the correct mode connectivities.

B. Global mode connectivity
To overcome the shortcoming of the perturbative approach, we do not use it to connect the phonon modes for all q points between q b and q e . Instead, we only apply the perturbative approach in the first stage to connect the modes for the first s q points, i.e., from q 0 , q 1 , · · · , to q s−1 , in order to construct the heads of the phonon branches, before switching in the second stage to an approach that is based on least-squares fits and exploits the smooth 'global' structure of the phonon bands. Beginning from q s , for each band index n, we use the previous r (obviously r ≤ s) values of ω 2 n,q i , for i = s − 1, s − 2, · · · , s − r, to form a quadratic fit for each n and predict the eigenvalue at q s . The subroutine for the least-squares fit of a polynomial is implemented according to Ref. [30]. These eigenvalues Ω 2 n,q s predicted from the least-squares fit are sorted and compared with the sorted exact eigenvalues ω 2 m,q s to establish the mode connectivities from q s−1 to q s . To construct the remainder of the phonon path, we repeat the same fitting approach stepwise to connect the modes from q s to q s+1 , from q s+1 to q s+2 , and so on until we finally connect all the modes from q p−1 to q p for every band. This completes the first pass of mode connectivities from q 0 to q p .
To eliminate the inappropriate connectivities that may have occurred during the first stage with the perturbative approach, we initiate a second pass of mode connectivities from q p to q 0 , where the last r eigenvalues for modes are taken from q p , q p−1 , · · · , q p−r+1 for fitting purposes, and establish connectivities between q p−r+1 and q p−r . Thereafter, the fitting approach is used to move stepwise from q p−r to q p−r−1 and so on until we finally connect q 1 to q 0 . We note here that since all q paths are independent of one another as far as the band indices are concerned, we have to apply the hybrid method afresh at the beginning of each path.
We now explain why the fitting approach handles the near or exact degenerate modes well. This is largely due to the fact that when a value is taken from the list of very close eigenvalues to be part of an array of r values for a quadratic fit, the predicted eigenvalue is likely to be very insensitive to any value chosen from the set of nearly degenerate eigenvalues since all values in the set are very close to one another. We also point out that our hybrid method is applicable to almost any phonon code (DFPT based or small-displacement based) since it is easy to generate the input dynamical matrices at any q points with the knowledge of interatomic force constants.

III. RESULTS
To assess the efficacy of the proposed hybrid method, we carry out density-functional theory (DFT) calculations within the local density approximation as implemented in the Vienna Ab initio Simulation Package (VASP) [31], with projected augmented-wave (PAW) pseudopotentials. The phonon calculations are performed using a small-displacement method [16,32]. The first system is a 2D graphene sheet with lattice constant a = 2.462 Å and a vacuum height of 12 Å. A 4 × 4 × 1 supercell is used. The cutoff energy for the plane-wave basis set is 500 eV. A k mesh of 6 × 6 × 1 is used for electronic relaxation. Figure 3 shows the mode connectivities of graphene along the Γ-K path with only the perturbative approach (i.e., we do not use the fitting approach). It is seen that the upper two bands have wrong connectivities near 1340 cm −1 . The inset in Fig. 3 shows why the perturbation approach fails to predict the correct ordering of eigenvalues at the degenerate point. Here, the eigenvalues in circles are obtained from exact diagonalization and the eigenvalues in crosses are obtained from perturbation theory. The perturbative approach connects well the modes from Γ to a q point corresponding to A. Because of the wrong ordering for the predicted eigenvalues  [14] has been taken account.
at a q point corresponding to B, a wrong assignment of modes occurs that results in visually detectable kinks in Fig. 3. However, when the hybrid approach is used, it is found that the fitting approach is able to overcome the problem highlighted in Fig. 3 and gives correct connectivities in the Γ-K path shown in Fig. 4. It is also observed that other band crossings along the K-M (one intersection point) and M -Γ (two intersection points) paths are correctly produced. Figure 5 shows the mode connectivities for the trigonal Bi 2 Se 3 with a r = 9.621 Å, α r = 24.64 • , which corre-  [14] has been taken account.
sponds to a conventional hexagonal cell of a = 4.105 and c = 27.973 Å. The supercell is 4×4×1 of the conventional hexagonal unit cell. The cutoff energy for the plane-wave basis set is 423.2 eV. A k mesh of 4×4×2 is used for electronic relaxation. The intricate band crossings at about 106 cm −1 , near the midpoint of the Γ-L path, appear to be very well captured by the hybrid method as shown in the inset of Fig. 5. We find that the hybrid method is robust enough to handle the three nearly degenerate points of about 160 cm −1 in the B-Z path. Figure 6 shows the mode connectivities of the hexagonal MoS 2 with a = 3.123 and c = 12.087 Å. The supercell size is 3 × 3 × 2. The cutoff energy for the plane-wave basis set is 700 eV. A k mesh of 2 × 2 × 1 is used for the electronic relaxation. We observe very good connectivities along all paths. For example, along the Γ-M path, the somewhat dispersive bands below 180 cm −1 give rise to many band crossings all of which the hybrid method is capable of resolving. For the high-frequency modes of about 420 cm −1 , four crossing points in a small region are shown to be described correctly. Note that all bands are doubly degenerate in the L-C path, which pose no difficulty for the hybrid method.

IV. CONCLUSION
In this paper we have proposed a two-stage hybrid method that uses local and global information for phonon mode connectivity. In the first stage of the algorithm, we exploit the local continuity of the eigenvalues within each band to construct part of the phonon dispersion path. We connect the actual eigenvalues at q i to actual eigenvalues at q i+1 through the use of first-order degenerate perturbation theory in which the perturbed eigenvalues at q i+1 are predicted from the eigenvectors at q i . In the second stage, we use a global polynomial fit of the partially constructed phonon dispersion path, which is based on leastsquares fits, to predict the eigenvalues and then connect the eigenvalues. This approach avoids local errors associated with degenerate or near-degenerate eigenvalues at or near the crossing points.
We find that a moderately fine q mesh is sufficient for the application of our proposed hybrid method, in contrast with a very fine q mesh used in the commonly adopted eigenvector similarity method. We note that our method is applicable to any phonon code (density-functional perturbation theory based or smalldisplacement based) since the phonon code generates the input dynamical matrices at any q point to our hybrid method. The robustness of our method has been demonstrated for graphene, Bi 2 Se 3 , and MoS 2 . We expect the use of our hybrid method could be be extended to handle mode connectivities in electronic band structures, or photonic modes in photonic crystals.