Complementary local-global approach for phonon mode connectivities

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 ${\omega }_{n,{\boldsymbol{q}}}^{2}$ and eigenvectors u_{n, q} , are governed by the equation

$$\begin{eqnarray}&&D({\boldsymbol{q}}){u}_{n,{\boldsymbol{q}}}={\omega }_{n,{\boldsymbol{q}}}^{2}{u}_{n,{\boldsymbol{q}}},\end{eqnarray} \tag{ 3 }$$

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}_{m{\boldsymbol{q}}}^{\dagger }{u}_{n,{\boldsymbol{q}}={\delta }_{{mn}}}$. At q + δ q , the corresponding dynamical matrix and its associated eigenvalues and eigenvectors are similarly given by

$$\begin{eqnarray}&&D({\boldsymbol{q}}+\delta {\boldsymbol{q}}){u}_{m,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}={\omega }_{m,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}^{2}{u}_{m,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}\end{eqnarray} \tag{ 4 }$$

for m = 1, 2, ⋯, N. For orientational purposes, we assume that none of the eigenvalues in equations (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.,

$$\begin{eqnarray}&&| {u}_{m,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}^{\dagger }{u}_{n,{\boldsymbol{q}}}| ={\delta }_{{mn}}+O(| \delta {\boldsymbol{q}}| ),\end{eqnarray} \tag{ 5 }$$

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 equation (5) is that its accuracy diminishes rapidly as ∣δ q ∣ increases.

To relate equations (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

$$\begin{eqnarray}&&V=D({\boldsymbol{q}}+\delta {\boldsymbol{q}})-D({\boldsymbol{q}}).\end{eqnarray} \tag{ 6 }$$

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.,

$$\begin{eqnarray}&&{u}_{n,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}=\displaystyle \sum _{j=0}^{\infty }{\lambda }^{j}{u}_{n,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}^{(j)}\end{eqnarray} \tag{ 7 }$$

where the zeroth-order term (j = 0) in the series is ${u}_{n,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}^{(0)}={u}_{n,{\boldsymbol{q}}}$. Likewise, the eigenvalue ${\omega }_{n,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}^{2}$ can be expanded as

$$\begin{eqnarray}&&{\omega }_{n,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}^{2}=\displaystyle \sum _{j=0}^{\infty }{\lambda }^{j}{\left({\omega }_{n,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}^{(j)}\right)}^{2}\end{eqnarray} \tag{ 8 }$$

where ${\omega }_{n,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}^{(0)}={\omega }_{n,{\boldsymbol{q}}}$.

The higher-order terms in the power series in equation (7) can be obtained using the familiar nondegenerate perturbation theory in quantum mechanics. For instance, the first-order in equation (7) (j = 1) term is

$$\begin{eqnarray}\begin{array}{rcl}{u}_{n,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}^{(1)} & = & \displaystyle \sum _{m\ne n}\displaystyle \frac{{u}_{m,{\boldsymbol{q}}}^{\dagger }{{Vu}}_{n,{\boldsymbol{q}}}}{{\omega }_{n,{\boldsymbol{q}}}^{2}-{\omega }_{m,{\boldsymbol{q}}}^{2}}{u}_{m,{\boldsymbol{q}}}\\ & = & \displaystyle \sum _{m\ne n}\displaystyle \frac{{u}_{m,{\boldsymbol{q}}}^{\dagger }D({\boldsymbol{q}}+\delta {\boldsymbol{q}}){u}_{n,{\boldsymbol{q}}}}{{\omega }_{n,{\boldsymbol{q}}}^{2}-{\omega }_{m,{\boldsymbol{q}}}^{2}}{u}_{m,{\boldsymbol{q}}}.\end{array}\end{eqnarray} \tag{ 9 }$$

The last equality follows because ${u}_{m{\boldsymbol{q}}}^{\dagger }D({\boldsymbol{q}}){u}_{n,{\boldsymbol{q}}}=0$ for m ≠ n. Therefore, to the first order, we obtain

$$\begin{eqnarray}&&{u}_{n,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}={u}_{n,{\boldsymbol{q}}}+\lambda {u}_{n,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}^{(1)}+O({\lambda }^{2}).\end{eqnarray} \tag{ 10 }$$

Similarly, the first-order term in equation (8) is

$$\begin{eqnarray}&&{({\omega }_{n,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}^{(1)})}^{2}={u}_{n,{\boldsymbol{q}}}^{\dagger }{{Vu}}_{n,{\boldsymbol{q}}}.\end{eqnarray} \tag{ 11 }$$

If we set λ = 1, then we can define the estimated eigenvector with band index n at q + δ q as

$$\begin{eqnarray}&&{u}_{n,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}^{{\prime} }={u}_{n,{\boldsymbol{q}}}+{u}_{n,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}^{(1)}\end{eqnarray} \tag{ 12 }$$

and its associated eigenfrequency as

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{n,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}^{2}={\omega }_{n,{\boldsymbol{q}}}^{2}+{\left({\omega }_{n,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}^{(1)}\right)}^{2}.\end{eqnarray} \tag{ 13 }$$

There are two ways that perturbation theory can be used for phonon mode connectivity given equations (12) and (13). The first way is to use equation (12) to sort the eigenmodes at q + δ q with respect to the eigenmodes at q using a generalized version of equation (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 equation (4), then we can assign it the band index n if it satisfies the overlap condition

$$\begin{eqnarray}&&| {u}_{m,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}^{\dagger }{u}_{n,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}^{{\prime} }| ={\delta }_{{mn}}+O({\lambda }^{2}).\end{eqnarray} \tag{ 14 }$$

where ${u}_{n,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}^{{\prime} }$ is as defined in equation (12). If we use only the first term for ${u}_{n,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}^{{\prime} }$ on the right hand side of equation (12) and substitute it to equation (14), we recover equation (2). The accuracy of the overlap condition of equation (14) may be increased by using a higher-order estimate for equation (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 ${{\rm{\Omega }}}_{n,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}^{2}$ from equation (13) to the exact eigenvalues ${\omega }_{m,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}^{2}$ from equation (4). We assigned the band index of n to ${\omega }_{m,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}^{2}$ (i.e., m = n) if the difference between ${{\rm{\Omega }}}_{n,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}^{2}$ and ${\omega }_{m,{\boldsymbol{q}}+\delta {\boldsymbol{q}}}^{2}$ is minimum.

We now provide more details on the implementation of the second way mentioned above. ¹ 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

$$\begin{eqnarray}&&{{\boldsymbol{q}}}_{i}={{\boldsymbol{q}}}_{b}+i({{\boldsymbol{q}}}_{e}-{{\boldsymbol{q}}}_{b})/p,\ i=0,1,\cdots ,\,p\end{eqnarray} \tag{ 15 }$$

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 ${\omega }_{n,{{\boldsymbol{q}}}_{i}}^{2}$, where n = 1, 2, ⋯, N. Due to possible band crossings, we cannot simply join the values of ${\omega }_{n,{{\boldsymbol{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 ${\omega }_{n,{{\boldsymbol{q}}}_{i}}$ and their associated eigenvectors ${u}_{n,{{\boldsymbol{q}}}_{i}}$, n = 1, 2, ⋯, N. From these eigenvalues and eigenvectors, we find the perturbed (or predicted) eigenvalues ${{\rm{\Omega }}}_{n,{{\boldsymbol{q}}}_{i+1}}^{2}$ at q _i+1, correct up to the first order, through ${{\rm{\Omega }}}_{n,{{\boldsymbol{q}}}_{i+1}}^{2}={\omega }_{n,{{\boldsymbol{q}}}_{i}}^{2}+{{\rm{\Delta }}}_{n,{{\boldsymbol{q}}}_{i}}$ where the correction terms ${{\rm{\Delta }}}_{n,{{\boldsymbol{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 equation (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 ${{\rm{\Delta }}}_{n,{{\boldsymbol{q}}}_{i}}$. Specifically this is achieved by first partitioning all N eigenvalues ${\omega }_{n,{{\boldsymbol{q}}}_{i}}^{2}$ 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}_{\alpha \beta }={u}_{\alpha ,{{\boldsymbol{q}}}_{i}}^{\dagger }{V}_{i}{u}_{\beta ,{{\boldsymbol{q}}}_{i}}$. A diagonalization of A gives d first-order corrections ${{\rm{\Delta }}}_{n,{{\boldsymbol{q}}}_{i}}$ for d eigenvalues ${\omega }_{n,{{\boldsymbol{q}}}_{i}}^{2}$ 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 equation (11).

From the ordering deduced from the sorted perturbed eigenvalues ${{\rm{\Omega }}}_{n,{{\boldsymbol{q}}}_{i+1}}^{2}$ and the ordering of the sorted actual eigenvalues ${\omega }_{m,{{\boldsymbol{q}}}_{i+1}}^{2}$ 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 ${\omega }_{n,{{\boldsymbol{q}}}_{i}}^{2}$ and ${\omega }_{m,{{\boldsymbol{q}}}_{i+1}}^{2}$ guided by ${{\rm{\Omega }}}_{n,{{\boldsymbol{q}}}_{i+1}}^{2}$, as can be seen in figure 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.

Zoom In Zoom Out Reset image size

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 figure 2(b). This switching has to be facilitated by the switching of the order of two predicted eigenvalues (as shown on the left of figure 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 figure 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.

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 ₀, q ₁, ⋯, 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 ${\omega }_{n,{{\boldsymbol{q}}}_{i}}^{2}$, 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 [29]. These eigenvalues ${{\rm{\Omega }}}_{n,{{\boldsymbol{q}}}_{s}}^{2}$ predicted from the least-squares fit are sorted and compared with the sorted exact eigenvalues ${\omega }_{m,{{\boldsymbol{q}}}_{s}}^{2}$ 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 ₀ 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 ₀, 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 ₁ to q ₀. 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.