# Optimized local modes for lattice dynamical applications

## I Introduction

Although the translational symmetry of a crystalline solid imposes a delocalized basis of Hamiltonian eigenstates (Bloch’s functions), it is sometimes advantageous to consider a transformation to a new set of basis functions with a local character. Beyond the mathematical equivalence (both sets span the same space of states), a local viewpoint is better suited for the analysis of concepts such as bonding which are eminently local in character. Recent work on electronic Wannier functions has shown the usefulness of a local representation in the chemical characterization of a given band subspace,  in the analysis of bonding topology in a disordered system,  and in more formal developments. 

The lattice dynamical problem is formally very similar to the electronic one: a set of Bloch eigenstates (normal modes) represents the collective vibrations of the atoms in the crystal. A basis change to a set of local displacement patterns (lattice Wannier functions or local modes) can in principle be achieved. So far, the main application for these local modes has been in the field of structural phase transitions. Typically, the behavior of a given dispersion branch or set of branches determines the essential instabilities of the system, and the associated degrees of freedom enter into the construction of an effective Hamiltonian which reproduces the relevant physics. Through the use of a localized basis set, the number of coupling terms in the effective Hamiltonian can be relatively small, easing the statistical mechanical treatment and the interpretation of the results. In particular, the anharmonic terms in the effective Hamiltonian can be kept local (on-site), in contrast with what happens in a reciprocal-space description.

This local mode approach has been used extensively in the past to gain an understanding of the behavior of complex systems, but until recently the local variables were treated as dummy degrees of freedom in a semi-empirical model, with their interactions fitted to reproduce the observed phenomena. In the last few years, a new approach, in which the effective Hamiltonian is parametrized on the basis of first principles calculations, has had great success in studies of the phase transition sequences in perovskite oxides.  Central to the parametrization process is the explicit construction of lattice Wannier functions, and two schemes have been proposed to carry it out. Zhong, Vanderbilt, and Rabe  (ZVR) used the structure of the zone-center soft mode in perovskite BaTiO<sub>3</sub> to construct symmetry-adapted highly localized local modes. Subsequently, Rabe and Waghmare  (RW) generalized this approach to reproduce the normal modes at several (typically, high symmetry) points of the Brillouin Zone.

While both the ZVR and RW approaches have been broadly successful in the specific problems for which they were conceived, in this paper we will argue that they are not completely satisfying in some respects. We will present a new procedure to generate lattice Wannier functions, an approach which makes use of the available symmetry information, produces local modes with a high degree of localization, enables a systematic improvement of their quality, and is straightforward to implement.

## II Method

We are interested in describing a relevant subspace $``$ of the full $`3Np`$-dimensional configuration space of a crystal with $`p`$ atoms per unit cell. Typically, we can choose $``$ as a complete band of dispersion branches (complete in the sense that it is invariant under the action of the space group of the crystal ). Associated to a branch $`j`$ is a set of normal modes ($`3Np`$-dimensional vectors) $`\{u_j^𝐤\}`$  which are eigenvectors of the Fourier transform of the force-constant matrix.  (The displacement in the $`\alpha `$ cartesian direction of the atom $`\kappa `$ in cell $`𝐥`$ is given explicitly by $`u_j^𝐤(𝐥,\kappa ,\alpha )`$.) The normal modes transform according to irreducible representations of the little groups $`G^𝐤`$. These representations, considered over the whole BZ, determine the band symmetry. The relevant subspace $``$ is spanned by all the $`\{u_j^𝐤\}`$ in the band, but it is clear that any transformation

$$\stackrel{~}{u}_j^𝐤=\underset{i=1}{\overset{n}{}}M_{ji}^𝐤u_i^𝐤$$
(1)   
will lead to a new basis of extended states which we will call Bloch modes. Here $`n`$ is the band dimension, the number of dispersion branches in the band.

Having thus specified the relevant subspace by means of Fourier space variables $`\{u_j^𝐤\}`$, the problem we tackle is the construction of a new basis $`\{w_j^𝐧\}`$ which is local, as opposed to extended, in character. Mathematically, the $`𝐤`$ label should be exchanged by a local label $`𝐧`$ associated to the different unit cells in the crystal. Translational symmetry takes the form:

$$w_j^𝐧(𝐥,\kappa ,\alpha )=w_j^{𝐧+𝐭}(𝐥+𝐭,\kappa ,\alpha ),$$
(2)   
which is trivially satisfied by the standard Wannier function form :

$$w_j^𝐧=\frac{1}{\mathrm{\Omega }}_{BZ}\mathrm{exp}(i\mathrm{𝐤𝐧})\stackrel{~}{u}_j^𝐤𝑑𝐤$$
(3)   
in which $`\mathrm{\Omega }`$ is the volume of the BZ. A high degree of localization means that the displacement $`w_j^𝐧(𝐥,\kappa ,\alpha )`$ should be very small or zero when $`𝐥`$ is a few lattice constants away from $`𝐧`$. The arbitrariness implicit in their definition (Eq. 1) means that the Wannier functions are non-unique, and a relatively large latitude then exists to tune their properties. In particular, the degree of localization has traditionally been the focus of great interest, and recently, Marzari and Vanderbilt  have succeeded in optimizing the matrices appearing in Eq. 1 to construct very localized electronic Wannier functions starting from the Bloch states. A restriction to unitary matrices resulted in an orthonormal basis of Wannier functions, and the optimization process led to symmetric-looking functions, even though no symmetry conditions were explicitly imposed.  In principle, such an approach should work for the vibrational problem, too. However, we prefer to take an alternate route which takes advantage of the knowledge of the band symmetry.

### A Symmetry requirements

As studied extensively in the literature,  one should supplement the translational constraints of Eq. 2 with another set of conditions which represent the transformational properties of the $`w_j^𝐧`$ under the effect of the point symmetry of the crystal. These are most easily discussed by introducing a symmetry-based definition of the center of a mode. Consider a Wyckoff set with representative site $`𝐫`$ and the set $`\widehat{G}_𝐫`$ of operations in $`G`$ that leave $`𝐫`$ invariant.  Given an irreducible representation $`\tau `$ of $`\widehat{G}_𝐫`$ with dimension $`d_\tau `$, any $`d_\tau `$ displacement patterns transforming with $`\tau `$ under the action of $`\widehat{G}_𝐫`$ are said to be centered in $`𝐫`$. It is then notationally more convenient to use a double index to label these patterns: $`w_{𝐫,s}`$ where $`s`$ ranges from 1 to $`d_\tau `$. The action of the elements of the space group $`G`$ on this set generates images at the rest of the positions in the Wyckoff set, i.e., $`d_rd_\tau `$ patterns $`w_{𝐫_i,s}^𝐧`$ per cell, where $`i`$ ranges from 1 to $`d_r`$ (the multiplicity of the Wyckoff set). This set of lattice functions is represented by the pair $`(𝐫,\tau )`$ and define a representation of $`G`$ which is called band representation

A necessary condition for the description of a relevant band subspace by means of these symmetry-adapted local modes is the equivalence of the band symmetry of $``$ and the band representation $`(𝐫,\tau )`$ (in particular this implies $`n=d_rd_\tau `$). More details about the choice of the correct $`(𝐫,\tau )`$ for a given $``$ are presented in the Appendix, where we also discuss the transformation properties of the corresponding $`\{\stackrel{~}{u}_j^𝐤\}`$. Incidentally, since Eq. 3 establishes a correspondence between lattice Wannier functions and Bloch modes, in what follows the latter can be also labeled by the site and representation indexes: $`\stackrel{~}{u}_{𝐫,s}^𝐤`$

### B Practical criterion for localization

A straightforward scheme to obtain lattice Wannier functions can be based on a direct use of Eq. 3, performing the BZ sum by means of any of the standard “special k-points” methods.  The quality of the subspace description can thus be systematically improved by simply using denser k-point sets. This approach can incorporate information about the normal modes throughout the whole Brillouin zone, as opposed to at just one point (as in the ZVR method), or at a very special set of high-symmetry k-points (as in the RW scheme).

As stated in the Introduction, it is highly desirable that the local mode basis functions be as localized as possible, in order to permit the consideration of only a few coupling terms in the effective Hamiltonian. From the point of view of real applications, a basis of Wannier functions which are not localized is not efficient, even if it spans $``$ perfectly. The form of Eq. 3 suggests a very simple heuristic criterion to achieve a high degree of localization for the lattice Wannier functions: choose the $`M^𝐤`$ matrices in such a way that the $`\stackrel{~}{u}^𝐤`$ Bloch vectors at different $`𝐤`$ add their contributions coherently at the center of the Wannier function. Interference effects can then be counted on to automatically dampen the amplitude of the displacements at sites away from the center.

Both the symmetry requirements and the localization condition can be formulated in the following way. Assume a $`(𝐫,\tau )`$ pair has been determined on the basis of band symmetry, and that we focus on the construction of local modes at cell $`𝐧`$=$`\mathrm{𝟎}`$. Consider a set of $`d_\tau `$ ($`3Np`$-dimensional) orthonormal vectors $`\{x_{𝐫,s}\}`$ which are centered in $`𝐫`$, transform with irrep $`\tau `$ and involve atoms in an orbit as close to $`𝐫`$ as possible.  The localization criterion is implemented by requiring that $`\stackrel{~}{u}_{𝐫,s}^𝐤`$ be orthogonal to $`x_{𝐫,t}`$ if $`st`$, and “parallel” (meaning that their scalar product is positive) if $`s=t`$ It can be seen that this condition fixes the form of the $`M^𝐤`$ matrices, up to overall normalization factors. In general, the $`M^𝐤`$ will not be unitary, with the result that two lattice Wannier functions at different cells $`w_{𝐫,s}^𝐧`$ and $`w_{𝐫^{},s^{}}^𝐧^{}`$ will not be orthogonal if the pairs $`(𝐫,s)`$ and $`(𝐫^{},s^{})`$ are not equal.  In the next section we will provide a simple worked example of the new construction scheme and will compare its results to those of other methods.

## III Examples and Discussion

In order to illustrate the scheme presented in the previous section, we will employ a two-dimensional model crystal with two different atoms which occupy the $`1a`$ (white) and $`1b`$ (black) Wyckoff positions of the plane group $`p4mm`$ (See Fig. 1 a)). A simple harmonic model for the force constants (with the couplings among the white atoms considered up to fourth nearest neighbors and the rest to first nearest neighbors, which corresponds to 6 independent parameters) gives the dispersion branches of panel b) in the figure. We will focus our attention on the two optical branches, which form a single band since they are essentially degenerate at the $`\mathrm{\Gamma }`$ and M points. These optical branches transform according to the decompositions

$$\begin{array}{ccc}\mathrm{\Gamma }\hfill & (4mm):\hfill & E\hfill \\ \mathrm{X}\hfill & (2mm):\hfill & B_1+B_2\hfill \\ \mathrm{M}\hfill & (4mm):\hfill & E,\hfill \end{array}$$
(4)   
in irreducible representations of the little co-groups at the high symmetry points. A simple application of the procedure spelled out in the Appendix shows that the band representation compatible with the above band symmetry is that represented by the pair $`(,E)`$, in which $`E`$ is a two-dimensional irreducible representation which turns out to be the vector representation of the point symmetry group at $``$. The set of $`\{x\}`$ vectors is then trivial to construct: as the “$``$” Wyckoff position is occupied, it is just enough to make $`x_{,1}`$ and $`x_{,2}`$ unit vectors attached to the central atom and pointing in the $`x`$ and $`y`$ cartesian directions, respectively. For this crystal structure, the simplest non-trivial set of special k-points is given by $`\{(1/8,1/8);(1/8,3/8);(3/8,3/8)\}`$. The explicit application of the localization criterion proceeds as follows. At each k-point in the set the normal modes are computed and the $`M^𝐤`$ matrices constructed. For example, at $`(1/8,3/8)`$, the normal modes are

$$\begin{array}{cc}u_1^𝐤=\hfill & (\mathrm{};0.23,0.93;\mathrm{})\hfill \\ u_2^𝐤=\hfill & (\mathrm{};0.84,0.19;\mathrm{}),\hfill \end{array}$$
(5)   
where the “$`\mathrm{}`$” refer to displacements on atoms other than the one at the center. The “coherent addition at the center” condition then becomes:

$$\begin{array}{c}0.23M_{11}+0.84M_{12}>0\hfill \\ 0.23M_{21}+0.84M_{22}=0\hfill \\ 0.93M_{11}0.19M_{12}=0\hfill \\ 0.93M_{21}0.19M_{22}>0,\hfill \end{array}$$
(6)   
and is satisfied by

$$M=\left(\begin{array}{cc}0.200& 0.980\\ 0.964& 0.264\end{array}\right),$$
(7)   
uniquely defined but for row-specific arbitrary factors. Since M is not unitary, the two optical Bloch vectors $`\stackrel{~}{u}_{,s}`$ at this k-point will not be orthogonal (although they can of course still be chosen to be normalized).

Once this procedure has been performed at every k-point in the set, the integral (sum) in Eq. 3 can be carried out to give the components of the lattice Wannier functions. Since the $`\{\stackrel{~}{u}^𝐤\}`$ determined by the localization criterion also satisfy the symmetry compatibility relations (see Appendix), the Wannier functions are symmetry-adapted. In Fig. 2 c) we show the displacements associated to the local mode $`w_{,1}`$, which transform as the first component of the vector representation $`(E)`$ of $`4mm`$. The degree of localization of these lattice Wannier functions can be gauged by computing the contribution to the total norm from a given shell around the center atom, as presented on Table I. Less than one per cent of the norm is outside the fourth shell (which corresponds roughly to the second-neighbor unit cells). If the integral in Eq. 3 is computed using a denser special-point set, the degree of localization is maintained, as can be seen by comparing the columns labeled “this work (3 k)” and “this work (10 k)” on Table I. This means that the quality of the local modes can be systematically improved while retaining a high degree of localization.

It is enlightening to compare this scheme to that of Rabe and Waghmare.  In the latter the analysis of the symmetry compatibility relations proceeds in the same way, and once the right $`(𝐫,\tau )`$ set has been identified, a series of orthonormal $`x`$ sets is constructed at successive shells centered on $`𝐫`$. The extent of the outermost shell fixes the localization of the Wannier functions by construction, and the actual atomic displacements are determined by fitting to the normal modes computed at a few high-symmetry points of the Brillouin Zone. In essence, the normal-mode information determines the weight assigned to each symmetry-adapted shell, so there is a tradeoff between the extent of the lattice Wannier functions and the amount of information from the real dispersion relations that can be used in the construction procedure. For example, in the PbTiO<sub>3</sub> work, Rabe and Waghmare found that adding information about the normal modes at the X point resulted in a less localized local mode than if only the $`(\mathrm{\Gamma },\mathrm{M},\mathrm{R})`$ set was used.  In contrast, our scheme can deal with the extra k-point without loss of localization: our local modes for PbTiO<sub>3</sub> using four high-symmetry points  are more localized than the best (three point) RW lattice Wannier functions.

It is clear that localization cannot be the main quality criterion for the construction of local modes. If it were, then the ZVR scheme, which uses only one high-symmetry k-point to construct a (very localized) lattice Wannier function, would be the method of choice. In fact, the real test for local mode sets is the degree to which they reproduce the energetics of the relevant subspace $``$. That is, in our case, the degree to which the dispersion relations of the effective Hamiltonian

$$H_{\mathrm{eff}}=H_{\mathrm{eff}}(Q_1,Q_2,\mathrm{},Q_{Nn})$$
(8)   
match the real dispersion branches associated with $``$.

In Eq. 8, the variables $`Q_i`$ are the amplitudes of the local mode variables, so that $`H_{\mathrm{eff}}`$ can be thought of as the “projection” of the complete Hamiltonian into the relevant subspace $``$ (which is typically considered as energetically decoupled from the rest of the configuration space of the crystal). The explicit form of $`H_{\mathrm{eff}}`$ will depend on the detailed structure of the lattice Wannier functions. In particular, the number of distinct coupling coefficients (representing the interaction of modes at different sites) which one should take into account in $`H_{\mathrm{eff}}`$ is determined by the spatial extent of the local modes.

We have constructed effective Hamiltonians for the model crystal for each of the three local-mode construction schemes discussed above (we obtain the coupling between $`w_{,s}^𝐧`$ and $`w_{,s^{}}^𝐧^{}`$ by calculating the energy associated to the crystal when it is distorted by just these modes). The original crystal Hamiltonian involved interactions up to fourth nearest neighbors for white atoms. Since the local modes involve basically displacements of the central white atom, we have kept the same range of interaction in $`H_{\mathrm{eff}}`$, but now referring of course to fourth nearest local modes. This amounts to using ten independent coupling coefficients; a larger number of parameters would not be reasonable in a practical application.

ZVR-style local modes are very localized and do not couple beyond the fourth neighbor shell, so the considered $`H_{\mathrm{eff}}`$ includes all the existing interactions. This can be seen on panel a) of Figure 3: the dispersion branches computed from $`H_{\mathrm{eff}}`$ match the exact ones at the $`\mathrm{\Gamma }`$ point. However, $`H_{\mathrm{eff}}`$ gives a poor description of the dispersion branches away from $`\mathrm{\Gamma }`$, as it should be expected in view of the construction procedure. (Incidentally, the inverse of Eq. 3 leads to Bloch modes which are not normalized to unity, except at the $`\mathrm{\Gamma }`$ point. The standard analysis of $`H_{\mathrm{eff}}`$ as given would lead to the low-lying dispersion branches in the figure. The higher branches are obtained by considering the corresponding generalized eigenvalue problem.) Panel b) shows that the $`H_{\mathrm{eff}}`$ constructed on the basis of RW local modes gives a good qualitative overall description of the dispersion, but fails to match the exact branches at the $`\mathrm{\Gamma }`$ point (as it should, given that this point was used in the construction scheme). The reason is that the local modes are more extended, and it is necessary to include couplings to further shells (at least up to seventh nearest neighbors)  for the match to be essentially perfect. This means that the RW scheme does not lead to efficient local modes, in the sense stated above.  The situation gets worse if more accuracy is needed in the overall description of the dispersion branches: the local modes turn out to be more extended, and even more coupling terms are needed in $`H_{\mathrm{eff}}`$.

In contrast, the local modes constructed following our heuristic criterion for localization do exhibit good efficiency (the dispersion branches do not change much when couplings to more than fourth nearest neighbors are included) and provide a very good qualitative match of the true branches throughout the BZ (Fig. 3 c)). (It should be noted that our construction scheme does not involve any high-symmetry points, hence the offset of the branches at $`\mathrm{\Gamma }`$ and M. We trade an overall good match for perfect accuracy at a few points.) Since a few coupling terms are enough to take into account the structure of the local modes, and the resulting $`H_{\mathrm{eff}}`$ provides a good fit to the true branches, our lattice Wannier functions are well suited for the local representation of the relevant subspace $``$. Moreover, they can be improved if needed by including more k-points in the integration set, with only a minor sacrifice in the compactness of the effective Hamiltonian.

We find these general conclusions to remain valid when more complicated interaction models are considered.

The practical application of our method of local mode construction to real materials requires the knowledge of the normal modes at general points of the Brillouin Zone. This information is easily obtained with modern linear-response codes without the need for large supercells. In the field of phase transitions, the use of this new scheme should enable the study of more complicated situations than those considered up to now. Competition of instabilities associated to different regions of the BZ or complications derived from anti-crossing phenomena are examples in which this method is bound to be useful. On the other hand, this work might provide an illustration of some of its theoretical underpinnings: the physical interpretation of the band representation associated to a dispersion band  or the symmetry-induced continuity of phonon spectra  are two instances of this.

## IV Conclusions

We have presented a straightforward scheme for the construction of very localized lattice Wannier functions, with explicit consideration of crystal symmetry. The new localization procedure enables a systematic improvement in the description of the relevant physics (by simply using denser sets of special k-points in the BZ integration) while still being quite efficient in regard to the number of coupling parameters needed in the effective Hamiltonian. Besides, the present method is straightforward to implement.

## Acknowledgements

We thank Karin Rabe, Philippe Ghosez, and David Vanderbilt for useful comments. This work was supported in part by the UPV research grant 060.310-EA149/95 and by the Spanish Ministry of Education grant PB97-0598. J.I. acknowledges fellowship support from the Basque regional government and thanks Agustin Válgoma for comments on the manuscript.

## Appendix

Let us study in more detail the equivalence between the band symmetry emerging from the transformation properties of the Bloch modes and the band representation associated to a $`(𝐫,\tau )`$ set.  This can be done by considering the action of $`G`$ on the $`(𝐫,\tau )`$ set of local modes and, consequently, on the associated Bloch modes $`\stackrel{~}{u}_{𝐫_i,s}^𝐤`$.

In order to proceed, we need to formulate the transformation properties of the modes $`w_{𝐫,s}^\mathrm{𝟎}`$ under the action of $`\{\overline{R}|\overline{𝐯}\}\widehat{G}_r`$. Since we consider $`\{\overline{R}|\overline{𝐯}\}`$ acting on the modes themselves and not on their components, we denote symmetry operations by the associated operators $`O\{\overline{R}|\overline{𝐯}\}`$. We have

$$O\{\overline{R}|\overline{𝐯}\}w_{𝐫,s}^\mathrm{𝟎}=\underset{h=1}{\overset{d_\tau }{}}D_{hs}^\tau (\{\overline{R}|\overline{𝐯}\})w_{𝐫,h}^\mathrm{𝟎}$$
(9)   
where $`𝐃^\tau (\{\overline{R}|\overline{𝐯}\})`$ is the matrix associated to $`\{\overline{R}|\overline{𝐯}\}`$ by irrep $`\tau `$. Now, we consider the rest of elements in the $`(𝐫,\tau )`$ set. They can be mathematically defined as

$$w_{𝐫_i,s}^𝐧:=O\{E|𝐧\}O\{R_i|𝐯_i\}w_{𝐫,s}^\mathrm{𝟎}$$
(10)   
where $`\{E|𝐧\}`$ is a lattice translation and $`\{R_i|𝐯_i\}`$ is one of the $`d_r`$ elements in $`G/G_r`$, which are chosen so that all the $`𝐫_i:=\{R_i|𝐯_i\}𝐫`$ lie in the same cell. The action of any $`\{R|𝐯\}G`$ on an arbitrary $`w_{𝐫_i,s}^𝐧`$ can be decomposed in: a lattice translation, a change of the center and a local transformation. Mathematically, this is expressed as

$`O\{R|𝐯\}(O\{E|𝐧\}O\{R_i|𝐯_i\}w_{𝐫,s}^\mathrm{𝟎})=`$ (11)   
$`O\{E|\{R|𝐯\}(𝐧+𝐫_i)𝐫_j\}O\{R_j|𝐯_j\}O\{\overline{R}|\overline{𝐯}\}w_{𝐫,s}^\mathrm{𝟎}`$ (12)   
where $`\{R_j|𝐯_j\}G/G_r`$ and $`\{\overline{R}|\overline{𝐯}\}\widehat{G}_𝐫`$ are univocally determined. Together with Eq. 9, this expression defines the band representation and, by using the inverse of Eq. 3, it can be written in the basis of Bloch modes. We obtain

$`O\{R|𝐯\}\stackrel{~}{u}_{𝐫_i,s}^𝐤=`$ (13)   
$`\mathrm{exp}(iR𝐤(\{R|𝐯\}𝐫_i𝐫_j)){\displaystyle \underset{h=1}{\overset{d_\tau }{}}}D_{hs}^\tau (\{\overline{R}|\overline{𝐯}\})\stackrel{~}{u}_{𝐫_j,h}^{R𝐤}`$ (14)   
By examining the representations this equation defines in the high symmetry k-stars, it can be easily checked whether the band representation $`(𝐫,\tau )`$ is equivalent to the band symmetry we want to describe. 

Once a convenient $`(𝐫,\tau )`$ set is chosen, Eq. 14 fixes the requirements on Bloch modes so that they lead to symmetry adapted local modes $`w_{𝐫_i,s}^𝐧`$. For pure translations, Eq. 14 reduces to Bloch theorem. Point symmetry determines the transformation properties of the $`d_\tau d_r`$ Bloch modes in each k-point and establishes the relationship of these with those in the rest of the k-star. However, Eq. 14 does not determine the form of the $`M^𝐤`$ matrices completely. For instance, in a general k-star ($`\overline{G}^𝐤=\{E\}`$) no condition is imposed on the choice of Bloch modes in a representative $`𝐤`$, though, once this is done, the modes in the rest of the star are fixed. This is the freedom we use in our construction procedure to get the localization of the modes.