Helmholtz Fermi surface harmonics: an efficient approach for treating anisotropic problems involving Fermi surface integrals

As already mentioned, the construction of the triangular mesh over the FS allows us to define an approximate discrete version of the Helmholtz equation and represents a very important step when defining all the variables involved in the computation, i.e., the triangle areas, internal angles, local curvature, etc.

For this purpose, we have implemented both the marching cube and the marching tetrahedra algorithms [6] (see appendix A). The marching cube algorithm is more popular and generates a smaller amount of triangles compared to the marching tetrahedra, but an important drawback of this method is that it includes non-manifold features (holes). We have found that although the marching tetrahedra method introduces a larger amount of triangle simplexes, the method shows to be much more robust.

The triangulation of the FS with any method produces a set of ${{n}_{t}}$ triangles whose vertices ${{{\mathbf{v}}}_{i}}$ (i = 1,..., ${{n}_{v}}$) are 3D ${\mathbf{k}}$-space vectors obtained by the triangulation process and ${{n}_{v}}$ are the number of vertices. Figure 1 illustrates the situation. Each vertex (${{{\mathbf{v}}}_{i}}$) has ${{N}_{n}}\left( {{{\mathbf{v}}}_{i}} \right)$ nearest neighbour vertices (${{{\mathbf{v}}}_{j}}$) and nearest neighbour triangles, denoted by ${{T}_{j}}\left( {{{\mathbf{v}}}_{i}} \right)$, $j=1,\ldots ,{{N}_{n}}\left( {{{\mathbf{v}}}_{i}} \right)$. The shaded area around each vertex ${{{\mathbf{v}}}_{i}}$ defines the 'barycenter control area', which is the area associated with each vertex (or 3D ${\mathbf{k}}$-space vector on the FS). The control area (${{S}_{i}}$) of a vertex ${{{\mathbf{v}}}_{i}}$ is calculated by considering the barycenter of the neighbouring triangles and the middle points of the vectors connecting the neighbouring triangles, hence the name. This area is the sum of $\frac{1}{3}$ of each neighbouring triangle area,

$$\begin{eqnarray}&&{{S}_{i}}=\mathop{\sum }\limits_{j=1,{{N}_{n}}\left( {{{\mathbf{v}}}_{i}} \right)}A\left[ {{T}_{j}}\left( {{{\mathbf{v}}}_{i}} \right) \right]/3,\end{eqnarray} \tag{ 5 }$$

where $A\left[ {{T}_{j}}\left( {{{\mathbf{v}}}_{i}} \right) \right]$ denotes the area of the nearest neighbour triangles ${{T}_{j}}\left( {{{\mathbf{v}}}_{i}} \right)$ of vertex ${{{\mathbf{v}}}_{i}}$. This barycenter control area is the simplest possible choice and provides a very simple quadrature formula for the integral of a function defined on the FS,

$$\begin{eqnarray}&&{{\int }_{{\mbox{FS}}}}{{{\mbox{d}}}^{2}}{{s}_{k}}f\left( k \right)\simeq \mathop{\sum }\limits_{i}{{S}_{i}}f\left( {{{\mathbf{v}}}_{i}} \right).\end{eqnarray} \tag{ 6 }$$

It is easily demonstrated that the quadrature formula given by (6) for the integral—and scalar product among functions—is absolutely equivalent to the procedure one would obtain by linearly interpolating the function f(k) within each triangle and integrating the linearly interpolated function. Thus, (6) is the generalization of the trapezoidal integration rule in a boundary free surface.

Zoom In Zoom Out Reset image size

Once the triangulated surface is constructed, and with all the above information at hand, the discrete version of the Laplace operator is numerically available only by identifying the nearest neighbour vertices (${{{\mathbf{v}}}_{j}}$) and triangles of a given vertex ${{{\mathbf{v}}}_{i}}$. In the linear approximation described above, the Laplace operator comes as a function only of the control area ${{S}_{i}}$ and the two opposite angles, ${{\alpha }_{i,j}}$ and ${{\beta }_{i,j}}$, of the triangles sharing the edge joining the vertices ${{{\mathbf{v}}}_{i}}$ and ${{{\mathbf{v}}}_{j}}$ (see figure 1). Thus, the two point centered formula for the second derivative in one dimension is generalized by [4, 5]

$$\begin{eqnarray}&&\nabla _{k}^{2}f\left( {\mathbf{k}} \right)\left| _{{\mathbf{k}}={{{\mathbf{v}}}_{i}}} \right.\simeq -\frac{1}{{{S}_{i}}}\mathop{\sum }\limits_{j=1,{{N}_{n}}\left( {{{\mathbf{v}}}_{i}} \right)}{{\Omega }_{i,j}}f\left( {{{\mathbf{v}}}_{j}} \right),\end{eqnarray} \tag{ 7 }$$

where

$$\begin{eqnarray}{{\Omega }_{i,j}}=\left\{ \begin{array}{ccccccccccccccc} -\frac{1}{2}\left[ {\rm cot} \left( {{\alpha }_{i,j}} \right)+{\rm cot} \left( {{\beta }_{i,j}} \right) \right] & i\ne j \\ {{\sum }_{i\ne j}}{{\Omega }_{i,j}} & i=j\;. \\ \end{array} \right.\end{eqnarray} \tag{ 8 }$$

Thus, since only nearest neighbour vertices contribute to (7) and (8), the discretized version of the Helmholtz equation for the HFSH and the BHFSH set are given, respectively, by the following two generalized highly sparse eigenvalue problems:

$$\begin{eqnarray}&&\frac{v\left( {{{\mathbf{v}}}_{i}} \right)}{{{S}_{i}}}\mathop{\sum }\limits_{j}{{\Omega }_{i,j}}\ \Phi _{L}^{j}={{\omega }_{L}}\Phi _{L}^{i},\end{eqnarray} \tag{ 9 }$$

and

$$\begin{eqnarray}&&\frac{1}{{{S}_{i}}}\mathop{\sum }\limits_{j}{{\Omega }_{i,j}}\ \Psi _{L}^{j}=\kappa _{L}^{2}\Psi _{L}^{i}.\end{eqnarray} \tag{ 10 }$$

Note that if the area of the control cells of all vertices (${{S}_{i}}$) was equal and the velocity [$v\left( {{{\mathbf{v}}}_{i}} \right)$] was constant in (9), we would then have a regular eigenvalue problem for the eigenfunctions ${{\Phi }_{L}}$ and the eigenvalues ${{\omega }_{L}}$ which are labelled by L. The same applies to (10). Moreover, since the ${{\Omega }_{i,j}}$ matrix is symmetric and the $\frac{{{S}_{i}}}{v\left( {{{\mathbf{v}}}_{i}} \right)}{{\delta }_{i,j}}$ operator is positive definite, the reality of the eigenvalues is automatically guaranteed. The linear problems in (9) and (10) are solved with the aid of the FEAST sparse eigenvalue solver library [7].

As mentioned, in our computational scheme the first step consists of finding the list of triangles sharing a given vertex of the triangulation. Periodic boundary conditions are implemented by imposing that all the vertices located at boundary planes of the BZ and differing by a reciprocal lattice vector ${\mathbf{G}}$ share the same triangular simplexes. Thus, the setup of the periodic boundary conditions implies not only accounting for all triangles shared by a point ${\mathbf{k}}$ (in a given boundary plane) but also adding to this list those triangle simplexes which are neighbouring an equivalent boundary plane and sharing the vertex ${\mathbf{k}}+{\mathbf{G}}$. When a given vertex is located at two boundary planes at the same time (an edge of the BZ) the same procedure described above is imposed twice.

Within the marching tetrahedra algorithm (appendix A), the entire BZ is sampled by a tetrahedral subdivision and the energy bands are explicitly calculated only at the corners of the tetrahedral simplexes. If a given tetrahedral simplex is found entirely above (or below) the Fermi level, the tetrahedron in question does not intersect the FS. When the corners of a tetrahedral simplex lie at both sides of the Fermi level, the band energies are linearly interpolated along the edges of the tetrahedron, allowing for a linear estimation of the positions of the three (or four) corners defining the intersection of the FS and the tetrahedral simplex. In this way, the initial 3D manifold is reduced to a 2D one.

Clearly, if the initial 3D mesh is relatively coarse or/and if the contribution of the second or higher order derivatives of the electron energy are appreciable, the error in the determination of the Fermi vectors are substantial. Once the 2D surface triangulation is obtained by the marching tetrahedra algorithm, we include a second step where the 2D ${\mathbf{k}}$ vectors are allowed to move along the normal to the FS such that the relaxed vector—and up to numerical precision—lies exactly at the FS. In this way, we iteratively improve the quality of the mesh by application of a generalized version of the Newton–Raphson type algorithm represented by the following iterative formula,

$$\begin{eqnarray}&&{{{\mathbf{k}}}_{\left( n+1 \right)}}={{{\mathbf{k}}}_{\left( n \right)}}-\nabla {{\epsilon }_{{{{\mathbf{k}}}_{\left( n \right)}}}}\frac{\left( {{\epsilon }_{{{{\mathbf{k}}}_{\left( n \right)}}}}-{{{\mbox{e}}}_{F}} \right)}{{{\left| \nabla {{\epsilon }_{{{{\mathbf{k}}}_{\left( n \right)}}}} \right|}^{2}}}.\end{eqnarray} \tag{ 11 }$$

In (11) the electron band indices are not shown for simplicity, but this relaxation scheme is applied to all ${\mathbf{k}}$ triangle vertices defining the FS and all bands composing the FS.

Note that the application of (11) is only affordable due to the low computational cost of the electronic band energies and velocities (energy gradients) through the Wannier method (see appendix B). All the cases that have been investigated have shown that 5–7 iterations (n) are more than sufficient to determine the Fermi wave vectors up to our numerical precision ∼${{10}^{-7}}$ eV.

The application of the Newton–Raphson method, as introduced above, improves the mesh quality and the accuracy of the eigenvalue problems represented by equations (9) and (10), but this additional step is not essential in our numerical scheme and may be ignored if a high quality triangulation is obtained by any other method.

Obviously, the refinement of the mesh improves the estimation of FS integrals in any theory or approach. Moreover, several important physical quantities involve the integral of the inverse of the velocity (energy gradient), the simplest of which being the DOS $\rho \left( E \right)$,

$$\begin{eqnarray}&&\rho \left( {{{\mbox{e}}}_{F}} \right)=\frac{2}{{{\Omega }_{{\mbox{BZ}}}}}\mathop{\sum }\limits_{n}\int \frac{{{{\mbox{d}}}^{2}}{{s}_{{\mathbf{k}},n}}}{{{v}_{{\mathbf{k}},n}}}\end{eqnarray} \tag{ 12 }$$

where ${{{\mbox{d}}}^{2}}{{s}_{{\mathbf{k}},n}}$ is the infinitesimal surface element and n denotes the electron band index (the 2 pre-factor assumes spin degeneracy). We approximate the integral in (12) by considering the barycenter control areas (${{S}_{i,n}}$) obtained by triangulation and explicitly calculating the velocities [${{v}_{n}}\left( {{{\mathbf{v}}}_{i}} \right)$] for all the relaxed vertices i and bands n,

$$\begin{eqnarray}&&\rho \left( {{{\mbox{e}}}_{F}} \right)\simeq \frac{2}{{{\Omega }_{{\mbox{BZ}}}}}\mathop{\sum }\limits_{n}\mathop{\sum }\limits_{i}\frac{{{S}_{i,n}}}{{{v}_{n}}\left( {{{\mathbf{v}}}_{i}} \right)}.\end{eqnarray} \tag{ 13 }$$

In the linear tetrahedron method [8] (12) is treated in essence by considering the electron velocities and vertex positions linearly interpolated in the inner volume of each tetrahedron.

In our numerical scheme two improvements are introduced for computing the integral in (12): (i) the electron velocities are explicitly calculated for all the 2D surface ${\mathbf{k}}$ point vertices and (ii) the triangular mesh is iteratively improved by forcing the triangle vertices to lie at the FS. Thus, the above method for estimating Fermi integrals may be considered as a surface specialized nonlinear version of the linear tetrahedron method.

The Gauss–Bonnet theorem is one the most prominent theorems in differential geometry connecting a local property of the surface such as the Gaussian curvature, and a global topological property such as the genus or the Euler characteristic.

The FS is a compact and periodic surface and does not present any boundary. With these restrictions, the Gauss–Bonnet theorem is stated as follows: for a given Fermi sheet n and $K\left( {\mathbf{k}} \right)$ being the Gaussian curvature at point ${\mathbf{k}}$ of the surface, the Euler characteristic χ is given by

$$\begin{eqnarray}&&\frac{1}{2\pi }\int {{{\mbox{d}}}^{2}}\;s_{{\mathbf{k}}}^{n}\;\;K\left( {\mathbf{k}} \right)={{\chi }^{n}}.\end{eqnarray} \tag{ 14 }$$

We have implemented the above formula as a quality test of our approach and we have checked that the result is completely independent of the choice of the BZ.

Numerically, the Gaussian curvature is given by

$$\begin{eqnarray}&&K\left( {\mathbf{k}} \right)\left| _{{\mathbf{k}}={{{\mathbf{v}}}_{i}}} \right.\simeq \frac{\left( 2\pi -\mathop{\sum }\limits_{j=1,{{N}_{n}}\left( {{{\mathbf{v}}}_{i}} \right)}{{\tau }_{i,j}} \right)}{{{S}_{i}}},\end{eqnarray} \tag{ 15 }$$

where ${{\tau }_{i,j}}$ denotes the angle between the triangle edges joining vertices ${{{\mathbf{v}}}_{i}}$ and ${{{\mathbf{v}}}_{j}}$, and ${{{\mathbf{v}}}_{i}}$ and ${{{\mathbf{v}}}_{j+1}}$ (see figure 1) and ${{S}_{i}}$ is defined in (5). Thus, one obtains that numerically (14) can be rewritten as

$$\begin{eqnarray}&&\mathop{\sum }\limits_{i}\left( 1-\frac{1}{2\pi }\mathop{\sum }\limits_{j=1,{{N}_{n}}\left( {{{\mathbf{v}}}_{i}} \right)}{{\tau }_{i,j}} \right)=\chi ,\end{eqnarray} \tag{ 16 }$$

which is the exact Descartes theorem on total angular defect of a polyhedron.

The genus of a single-sheet surface is given by

$$\begin{eqnarray}&&g=1-\chi /2.\end{eqnarray} \tag{ 17 }$$

While χ is additive for multiple-sheet surfaces, g must be considered separately for each sheet. It is worth mentioning that we obtain the above topological numbers by a direct application of (16), and that we obtain an integer number up to our numerical accuracy.

Even though the genus of several FSs such as lithium are trivial (g = 0), the same is not true even for a free-electron-like material such as Cu, where we find $\chi =-6$ (and g = 4).

Although Allen's FSH [1] and the HFSH presented in this article are defined very differently in the sense that the FSH are constructed as explicit orthogonal polynomials and the HFSH set is generated as a solution of the Helmholtz equation defined on the FS ((3) and (9)), both sets present very important similarities in their properties. First, being a solution of a second order differential equation, the HFSH set is orthogonal; and second, our choice for the normalization in (4) is deliberately chosen equal to that introduced in [1].

Allen rewrote the anisotropic Boltzmann and Eliashberg equations in terms of the polynomial FSH set. One concludes that as the scalar product of the HFSH set is defined exactly as for the FSH, the expressions derived by Allen for the Boltzmann transport and the anisotropic superconductivity in terms of FSH are still valid for the HFSH set presented in this work. This is the reason why in this article we concentrate on the calculation and description of the HFSH set and we refer to [1] and [9] for a detailed treatment of the Boltzmann and Eliashberg equations using FSH. In the next part we address the problem of function representation in terms of HFSH.

The HFSH set allows us to efficiently represent any function defined on the FS, but more importantly, the HFSH modes allow us to express integro-differential equations involving quasi-elastic scattering processes much more economically and probably in a physically more transparent way. Good examples of anisotropic functions defined on the FS are for instance the electron lifetime, $\tau \left( {\mathbf{k}} \right)$, the electron mass enhancement, $\lambda \left( {\mathbf{k}} \right)$, the superconducting gap, $\Lambda \left( {\mathbf{k}} \right)$, electron velocity components, etc. Similar to any integral transformation, the smoother the function to be represented, the more efficient the HFSH representation method. In the next subsections, we formalize the problem of functional representation considering the HFSH modes.

Let us consider any of these physical properties as a generic function, $F\left( {\mathbf{k}} \right)$, and let us rewrite this function in terms of the HFSH set,

$$\begin{eqnarray}&&F\left( {\mathbf{k}} \right)=\mathop{\sum }\limits_{L}{{c}_{L}}\left( F \right){{\Phi }_{L}}\left( {\mathbf{k}} \right).\end{eqnarray} \tag{ 18 }$$

We refer to (18) as the HFSH representation of the ${\mathbf{k}}$ dependent function $F\left( {\mathbf{k}} \right)$ defined on the FS.

The scalar product and normalization of the HFSH $\left\{ {{\Phi }_{L}} \right\}$ are defined by (4) in momentum space, thus the components ${{c}_{L}}\left( F \right)$ of the expansion become directly accessible, the coefficients of the expansion having the same units as the original $F\left( {\mathbf{k}} \right)$ function,

$$\begin{eqnarray}&&{{c}_{L}}\left( F \right)=\left\langle {{\Phi }_{L}}|F \right\rangle \equiv \frac{\int \frac{{{{\mbox{d}}}^{2}}{{s}_{{\mathbf{k}}}}}{v\left( {\mathbf{k}} \right)}{{\Phi }_{L}}\left( {\mathbf{k}} \right)F\left( {\mathbf{k}} \right)}{\int \frac{{{{\mbox{d}}}^{2}}{{s}_{{\mathbf{k}}}}}{v\left( {\mathbf{k}} \right)}}.\end{eqnarray} \tag{ 19 }$$

Of course, as the HFSH set is orthogonal, the scalar product of two functions ${{F}_{1}}$ and ${{F}_{2}}$ may be written conveniently in terms of the HFSH components of these functions,

$$\begin{eqnarray}&&\left\langle {{F}_{1}}|{{F}_{2}} \right\rangle =\mathop{\sum }\limits_{L}{{c}_{L}}\left( {{F}_{1}} \right){{c}_{L}}\left( {{F}_{2}} \right)=\frac{\int \frac{{{{\mbox{d}}}^{2}}{{s}_{{\mathbf{k}}}}}{v\left( {\mathbf{k}} \right)}{{F}_{1}}\left( {\mathbf{k}} \right){{F}_{2}}\left( {\mathbf{k}} \right)}{\int \frac{{{{\mbox{d}}}^{2}}{{s}_{{\mathbf{k}}}}}{v\left( {\mathbf{k}} \right)}}.\end{eqnarray} \tag{ 20 }$$

The procedure of transforming a function $F\left( {\mathbf{k}} \right)$ into a discrete set of coefficients ${{c}_{L}}\left( F \right)$ is very similar to ordinary Fourier transformation or the spherical harmonics expansion in the unit sphere. For sufficiently well behaved functions we should expect that a relatively small amount of HFSH modes is enough but, in any case, one has control of the desired accuracy by tuning the energy cutoff.

The product of two functions ${{F}_{1}}\left( {\mathbf{k}} \right){{F}_{2}}\left( {\mathbf{k}} \right)$ defined in ${\mathbf{k}}$ space may be represented as a function of the HFSH components of each of these functions separately [${{c}_{L}}\left( {{F}_{1}} \right)$ and ${{c}_{L}}\left( {{F}_{2}} \right)$] with the following relation that generalizes the convolution theorem in Fourier analysis

$$\begin{eqnarray}&&{{F}_{1}}\left( {\mathbf{k}} \right){{F}_{2}}\left( {\mathbf{k}} \right)=\sum 1,{{L}_{2}},LL{{\Xi }_{L;{{L}_{1}},{{L}_{2}},L}}\;\;{{c}_{{{L}_{1}}}}\left( {{F}_{1}} \right)\;\;{{c}_{{{L}_{2}}}}\left( {{F}_{2}} \right)\;{{\Phi }_{L}}\left( {\mathbf{k}} \right),\end{eqnarray} \tag{ 21 }$$

where the matrix elements ${{\Xi }_{L;{{L}_{1}},{{L}_{2}}}}$ play the same role as the Clebsch–Gordan coefficients for ordinary spherical harmonics,

$$\begin{eqnarray}&&{{\Xi }_{L;{{L}_{1}},{{L}_{2}}}}=\left\langle {{\Phi }_{L}}|{{\Phi }_{{{L}_{1}}}}{{\Phi }_{{{L}_{2}}}} \right\rangle .\end{eqnarray} \tag{ 22 }$$

The coefficients ${{\Xi }_{L;{{L}_{1}},{{L}_{2}}}}$ are symmetric with respect to the permutations of the $L,{{L}_{1}},{{L}_{2}}$ indices, as was also found by Allen for the FSH set [1]. Some elementary properties of ${{\Xi }_{L;{{L}_{1}},{{L}_{2}}}}$ are (the first element of the HFSH set is L = 1)

$$\begin{eqnarray*}&&{{\Xi }_{L;1,1}}={{\delta }_{L,1}}\quad {{\Xi }_{1;{{L}_{1}},{{L}_{2}}}}={{\delta }_{{{L}_{1}},{{L}_{2}}}}.\end{eqnarray*}$$

The $\Xi$ coefficients allow us to expand the product of two $\Phi$ functions in terms of a simple linear combination of HFSH elements and vice-versa,

$$\begin{eqnarray}\begin{array}{rcl} {{\Phi }_{L}}\left( {\mathbf{k}} \right) & = & \mathop{\sum }\limits_{{{L}_{1}},{{L}_{2}}}{{\Xi }_{L;{{L}_{1}},{{L}_{2}}}}{{\Phi }_{{{L}_{1}}}}\left( {\mathbf{k}} \right){{\Phi }_{{{L}_{2}}}}\left( {\mathbf{k}} \right) \\ {{\Phi }_{{{L}_{1}}}}\left( {\mathbf{k}} \right){{\Phi }_{{{L}_{2}}}}\left( {\mathbf{k}} \right) & = & \mathop{\sum }\limits_{L}{{\Xi }_{L;{{L}_{1}},{{L}_{2}}}}{{\Phi }_{L}}\left( {\mathbf{k}} \right). \\ \end{array}\end{eqnarray} \tag{ 23 }$$

The scalar product of pairs of HFSH may also be expressed in terms of the $\Xi$ matrix elements,

$$\begin{eqnarray}&&\left\langle {{\Phi }_{{{L}_{1}}}}{{\Phi }_{{{L}_{2}}}}|{{\Phi }_{{{L}_{2}}}}{{\Phi }_{{{L}_{3}}}} \right\rangle =\mathop{\sum }\limits_{L}{{\Xi }_{L;{{L}_{1}},{{L}_{2}}}}{{\Xi }_{L;{{L}_{2}},{{L}_{3}}}}.\end{eqnarray} \tag{ 24 }$$

In principle, one would calculate all the ${{\Xi }_{L;{{L}_{1}},{{L}_{2}}}}$ coefficients only once and try to reduce all redundancies as much as possible.

Let us now consider the HFSH representation of a matrix element of any physical magnitude with two momentum indexes. The EP matrix elements ${{g}_{{\mathbf{k}},{{{\mathbf{k}}}^{^{\prime} }}}}$ as well as many other physical quantities, such as the non-local self-energy or the impurity scattering matrix elements, need to be represented generally in terms of a pair on electron momenta ${\mathbf{k}}$ and ${{{\mathbf{k}}}^{^{\prime} }}$. In all these cases we would follow a similar procedure,

$$\begin{eqnarray}&&{{g}_{{\mathbf{k}},{{{\mathbf{k}}}^{^{\prime} }}}}=\mathop{\sum }\limits_{L,{{L}^{^{\prime} }}}{{c}_{L,{{L}^{^{\prime} }}}}\left( g \right){{\Phi }_{L}}\left( {\mathbf{k}} \right){{\Phi }_{{{L}^{^{\prime} }}}}\left( {{{\mathbf{k}}}^{^{\prime} }} \right),\end{eqnarray} \tag{ 25 }$$

where

$$\begin{eqnarray}&&{{c}_{L,{{L}^{^{\prime} }}}}\left( g \right)=\frac{\int \frac{{{{\mbox{d}}}^{2}}{{s}_{{\mathbf{k}}}}}{v\left( {\mathbf{k}} \right)}\frac{{{{\mbox{d}}}^{2}}{{s}_{{{{\mathbf{k}}}^{^{\prime} }}}}}{v\left( {{{\mathbf{k}}}^{^{\prime} }} \right)}\ {{g}_{{\mathbf{k}},{{{\mathbf{k}}}^{^{\prime} }}}}\ {{\Phi }_{L}}\left( {\mathbf{k}} \right){{\Phi }_{{{L}^{^{\prime} }}}}\left( {{{\mathbf{k}}}^{^{\prime} }} \right)}{{{\left( \int \frac{{{{\mbox{d}}}^{2}}{{s}_{{\mathbf{k}}}}}{v\left( {\mathbf{k}} \right)} \right)}^{2}}}.\end{eqnarray} \tag{ 26 }$$

Of course the function ${{g}_{{\mathbf{k}},{{{\mathbf{k}}}^{^{\prime} }}}}$ should be reasonably smooth for a good quality representation in terms of the HFSH as described above. When the quantity in question is a scattering amplitude or a complex matrix element, one must first fix the complex arbitrary phases of ${{g}_{{\mathbf{k}},{{{\mathbf{k}}}^{^{\prime} }}}}$. If time reversal and inversion symmetries are both present, these phases are easily removed [10], but more generally, one is forced to fix these phases by a Wannier procedure [10–12]. Alternatively, one could consider the absolute values of the matrix elements $|{{g}_{{\mathbf{k}},{{{\mathbf{k}}}^{^{\prime} }}}}{{|}^{2}}$.

The above algebraic machinery has great potential in restating integro-differential problems defined on the FS. We refer to [1, 2] for a detailed description of the procedure for transforming the Boltzmann transport and the Eliashberg theory for the FSH, which would be completely valid for our HFSH set.

One of the main problems encountered when trying to characterize an anisotropic physical quantity defined on the FS is that the only accessible method is a graphical representation through a colour code, which requires the computation of that quantity on a large amount (of the order of ${{10}^{5}}$) of ${\mathbf{k}}$ vectors defining the FS. However, this method is mainly visual and not quantitative.

The representation of a function $F\left( {\mathbf{k}} \right)$ using the HFSH expansion coefficients ${{c}_{L}}\left( F \right)$ enables a direct quantitatively description of the anisotropy. Indeed, the HFSH method allows us to tabulate numerically any complex anisotropic quantity by means of about 10$^{2}$coefficients, the most representative being those corresponding to the modes with the lowest energy.

Since ${{\Phi }_{L=1}}\left( F \right)=1$, the first expansion coefficient, ${{C}_{L=1}}\left( F \right)$, yields directly the FS average of the function $F\left( {\mathbf{k}} \right)$. In materials with spherical or nearly spherical symmetry, such as Li or Cu (see section 4.1), the next three modes [${{\Phi }_{L}}\left( F \right),L=2,3,4$] are very similar to the ${{p}_{x}}$, ${{p}_{y}}$ and ${{p}_{z}}$ spherical harmonics. In general, we find that a few number ($\sim 20$) of coefficients is sufficient for capturing the most significant part of the anisotropy of any ${\mathbf{k}}$ dependent function.

In this sense, the numerical tabulation appears to be a potentially important application of the HFSH set. Let us consider as a standard example, the anisotropic EP mass enhancement, ${{\lambda}{({\mathbf{k}})}}$, or the momentum dependent lifetime, ${{\tau }{({\mathbf{k}})}}$. Following a standard procedure, if one wanted to compare two different calculations of ${{\lambda}{({\mathbf{k}})}}$, for instance, obtained using two different computation methods, one would need to compare the ${\mathbf{k}}$ dependent data set point by point. This is, obviously, not practical and results on a high symmetry direction might be practically checked, if at all.

In a HFSH mode expansion a list of a few numerical coefficients [${{c}_{L}}\left( F \right)$] would be sufficient to compare, at least, the grossest details of the directional dependence of any magnitude, and the tabulation of any anisotropic quantity would then be easily accessible.

3.8.1. Denoising, filtering and mismatch error analysis

The HFSH is a complete set and the finest details of any quantity are accessible only by increasing the cutoff energy (up to the triangular grid capability). However, it may happen that a quantity calculated on the FS is accompanied by a noisy background, which is a situation that could be standard experimentally.

Let us suppose that we have calculated a physical property $F\left( {\mathbf{k}} \right)$ explicitly for all the vertex ${\mathbf{k}}$ points of our triangular grid describing the FS. Consider a finite cutoff (${{E}_{c}}$) for the expansion of $F\left( {\mathbf{k}} \right)$ in terms of the HFSH set using (18):

$$\begin{eqnarray}&&{{\tilde{F}}^{{{N}_{L}}}}\left( {\mathbf{k}} \right)=\mathop{\sum }\limits_{L=1}^{{{N}_{L}}}{{c}_{L}}\left( F \right){{\Phi }_{L}}\left( {\mathbf{k}} \right),\end{eqnarray} \tag{ 27 }$$

where ${{N}_{L}}$ denotes the number of modes such that ${{\omega }_{L}}<{{E}_{c}}$. We then obtain a smoothed function ${{\tilde{F}}^{{{N}_{L}}}}\left( {\mathbf{k}} \right)$ where the fine details smaller than a wave length

$$\begin{eqnarray}&&{{\lambda }_{c}}\sim 2\pi \sqrt{\frac{{{C}_{L=1}}}{{{\omega }_{L}}}}\end{eqnarray} \tag{ 28 }$$

are filtered.

A measure of the mismatch error when considering only a finite set of HFSH modes, ${{N}_{L}}$, might be obtained by the FS integral

$$\begin{eqnarray}&&\epsilon \left( {{N}_{L}} \right)=\frac{\int \frac{{{{\mbox{d}}}^{2}}{{s}_{{\mathbf{k}}}}}{{{v}_{{\mathbf{k}}}}}\left| F\left( {\mathbf{k}} \right)-{{{\tilde{F}}}^{{{N}_{L}}}}\left( {\mathbf{k}} \right) \right|}{\int \frac{{{{\mbox{d}}}^{2}}{{s}_{{\mathbf{k}}}}}{{{v}_{{\mathbf{k}}}}}\left| F\left( {\mathbf{k}} \right) \right|},\end{eqnarray} \tag{ 29 }$$

which is approximated by a simple linear quadrature formula as in (6).

Our first example applications are bcc-Li and fcc-Cu at ambient pressure, which are textbook examples of free-electron-like systems. The FS of these materials is approximately spherical and the HFSH functions obtained from (9) and (10) should be expected to approximate the ordinary spherical harmonics.

Figures 2(a) and (b) show the calculated FS of bcc-Li and fcc-Cu, respectively, considering a momentum space division of 40$^{3}$ in both cases. The marching tetrahedron algorithm (as introduced in appendix A) produced 27962 (bcc-Li) and 24582 (fcc-Cu) triangle vertices defining the FS, which were relaxed until they were located within a window of $\left| {{\epsilon }_{{\mathbf{k}}}}-{{E}_{F}} \right|<{{10}^{-6}}$ eV. The colour code in figures 2(a) and (b) corresponds to the absolute value of the electron velocity which enters (3) and (9). Although both bcc-Li and fcc-Cu present an almost spherical FS, we observe that the anisotropy of the electron velocity is not negligible.

Zoom In Zoom Out Reset image size

Figures 3 and 4 present the calculated 16 lowest energy HFSH modes for bcc-Li and fcc-Cu, respectively. The states are shown by following the same degeneracy ordering as if they were the usual spherical harmonics, i.e., the first row corresponds to the constant s-like state, second one corresponds to the ${{p}_{x}}$, ${{p}_{y}}$, and ${{p}_{z}}$-like modes, third row to the d-like set, and so on. The correspondence with the ordinary spherical harmonics is direct because of the simplicity of this surface. The energy dependence of the HFSH for bcc-Li and fcc-Cu is shown with filled circles in the inset of figures 3 and 4, respectively. We observe that the first non-trivial three states of p-like character (second row) are found to be degenerate like the spherical harmonics. However, the original five-fold degeneracy of the ordinary d spherical harmonics is broken in both cases generating three degenerate states plus an additional two dimensional degenerate subspace. Indeed, in these systems the crystal symmetry includes a p-like symmetry but not a d-like one and the effect appears as a crystal field effect acting on a spherically symmetric system. Similarly, the seven-fold original degeneracy is lifted in bcc-Li into a set of degenerate subspaces of three, one and three dimensions respectively ($3+3+1$ in bcc-Cu). For comparison we have also shown in the inset of figure 3 and 4 (dashed lines) the degenerate eigenvalues of (3) when the HFSH are taken to be ordinary spherical harmonics

$$\begin{eqnarray}&&\omega =\frac{l\left( l+1 \right){{v}_{F}}}{k_{F}^{2}},\end{eqnarray} \tag{ 30 }$$

where l denotes the angular momentum, ${{v}_{F}}$ is the mean velocity at the FS and ${{k}_{F}}$ is the mean radius of each of the nearly spherical Fermi surfaces displayed in figure 2. Note that in Li (figure 3), the lowest HFSH energies (filled circles) reproduce very well the spherical harmonics eigenvalues (dashed lines) for l = 1 and l = 2. Even for l = 3, the degeneracy of the HFHS eigenvalues is broken around the value given by (30). The lifting of the degeneracy for l = 3 is stronger in the case of Cu (figure 4), where we still find a reasonable agreement between the HFSH energies and (30) for l = 2. However, for l = 1 the ideal spherical harmonics energy deviates from the HFSH eigenvalue, which we attribute to the fact that the necks of the FS of Cu represent a strong perturbation to the sphere, mostly noticeable at long wavelengths.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 3 and 4 demonstrate several general features of the HFSH functions which are common to those of other more complex systems as well. The first mode (L = 1) is a constant function over the surface and its energy is equal to zero. This is clear from (9) (and (10) for BHFHS), as the diagonal elements of the Laplace operator are equal to the sum of the rest of the elements in the same row and therefore the summation is zero and any constant function multiplied by the discrete Laplace operator ($\Omega$ in (8)) is null. Thus, the constant function is always the lowest energy member in any HFSH set.

For higher energies the wavelength of the mode oscillations are shorter and the energy increases, to a very good approximation, linearly with the number of modes. For both, the HFSH and HFSHB sets, we find that as the mode number $L\to \infty$

$$\begin{eqnarray}&&{{\omega }_{L}}\simeq \frac{{{\pi }^{2}}{{v}_{F}}}{S}L\end{eqnarray} \tag{ 31 }$$

and

$$\begin{eqnarray}&&\kappa _{L}^{2}\simeq \frac{{{\pi }^{2}}}{S}L\end{eqnarray} \tag{ 32 }$$

Equations (31) and (32) work very well already for $L\sim 100$ in all materials—and Fermi sheets—treated in this work. These relations allows us to estimate the number of modes required for a given cutoff energy.

Figure 5 shows the application of the HFSH mode analysis to the cartesian velocity component ${{v}_{x}}\left( {\mathbf{k}} \right)$ (top), the multiplication of the x and y components, ${{v}_{x}}\left( {\mathbf{k}} \right){{v}_{y}}\left( {\mathbf{k}} \right)$, (middle) and the modulus of the velocity $|v\left( {\mathbf{k}} \right)|$ (bottom), defined on the FS of fcc-Cu. We show the absolute value of the coefficients ${{c}_{L}}$, defined in (18) and (19), for each magnitude using two different tetrahedral samplings of 20$^{3}$ (dashed red) and 40$^{3}$ (solid black) in the triangulation of the FS. The HFSH modes are real and completely determined up to a sign factor. This is the reason why we show the absolute value of the coefficients.

Zoom In Zoom Out Reset image size

The first component (L = 1) of the modulus of the velocity (bottom) gives just the average value, as this mode is constant. That fcc-Cu is a free-electron-like material is confirmed because the L = 1 component is the strongest contribution by far. Furthermore, the HFSH spectrum quantifies to which extent the modulus of the FS velocity is isotropic (or anisotropic). The strength of the peaks decreases rapidly and only some of the HFSH modes contribute significantly. We identify these modes graphically (L = 1, 16, 34, 72 and 120) in the inset of figure 5 (bottom).

The HFSH mode analysis for the x component of the velocity (top) shows that the first component is equal to zero ${{c}_{L=1}}\left( {{v}_{x}} \right)=0$, while the next two modes ${{c}_{L=2,3}}\left( {{v}_{x}} \right)$ are the most important. This is consistent with the interpretation of the first three non-trivial degenerate modes to be similar to the ordinary p-like spherical harmonics (see figure 4) and the x direction (in our coordinate system) appears basically as a linear combination of the L = 2 and L = 3 HFSH modes. The average values over the FS of the ${{v}_{x}}$ and ${{v}_{x}}{{v}_{y}}$ functions are zero because these functions are obviously odd, and therefore ${{c}_{L=1}}\left( {{v}_{x}} \right)=0$ and ${{c}_{L=1}}\left( {{v}_{x}}{{v}_{y}} \right)=0$, in both cases.

The ${{c}_{L}}$ components of the HFSH for the dense 40$^{3}$ (solid black) and coarser 20$^{3}$ (dashed red) samplings follow each other very closely. Of course, the HFSH energies are slightly reordered going from a coarse to a denser mesh. This is specially clear in the middle panel where we show the results of the ${{v}_{x}}{{v}_{y}}$ which has a more complicated spatial structure comparing to ${{v}_{x}}$ and $|v|$ and the amplitudes at higher energies are stronger. In any case, we find that at $L\sim 180$ the HFSH amplitudes (${{c}_{L}}$) are already about 20 times weaker than the maximum values of ${{c}_{L}}$ at lower L.

Table 2 shows the analysis of the mode number dependence of the mismatch error for the x cartesian component of the velocity (${{v}_{x}}$), the product of the x and y components (${{v}_{x}}{{v}_{y}}$) and the modulus of the velocity ($|v|$) for the two different tetrahedral samplings considered in the marching tetrahedron method (20$^{3}$ and 40$^{3}$) for fcc-Cu. We observe that the mismatch error is reduced by considering a denser triangular mesh, specially when a relatively large number of modes (see the rows corresponding to ${{N}_{L}}=401$ and 701 in table 2) is used. The mismatch error for ${{v}_{x}}{{v}_{y}}$ is the largest because this function shows a richer spatial structure than ${{v}_{x}}$ or $|v|$ (see figure 5). Expanding the functions with ${{N}_{L}}\sim 701$ modes is sufficient in order to obtain an estimated error below 1% with a tetrahedral sampling of 40$^{3}$.

Table 2.  Mismatch error (percentage) obtained using (29) for ${{v}_{x}}$, ${{v}_{x}}{{v}_{y}}$ and modulus of the velocity ($|v|$) in fcc-Cu as a function of the HFSH number of modes (${{N}_{L}}$).

	20$^{3}$			40$^{3}$
${{N}_{L}}$	${{v}_{x}}$	${{v}_{x}}{{v}_{y}}$	$|v|$	${{v}_{x}}$	${{v}_{x}}{{v}_{y}}$	$|v|$
101	6.9	20.5	4.6	6.8	20.0	4.3
201	3.1	6.9	2.1	2.7	5.5	1.8
401	1.5	3.5	1.1	0.9	1.8	0.6
701	1.2	2.6	0.5	0.4	0.9	0.3

In fcc-Cu the number of vertices in the triangulation of the relatively coarse (20$^{3}$) and dense (40$^{3}$) tetrahedral meshes were $N_{v}^{{{20}^{3}}}=6018$ and $N_{v}^{{{40}^{3}}}=24582$, respectively. One can estimate the data-storage saving when using the HFSH-mode representation instead of the usual vertex (${\mathbf{k}}$-space) representation. For example, for a function such as $F={{v}_{x}}{{v}_{y}}$ and with a target error below 1%, a saving factor of the order of $\alpha \left( {{40}^{3}} \right)$∼$1/40$ would be obtained with a 40$^{3}$ tetrahedral mesh.

Hex-MgB$_{2}$ is an important example of strong EP coupling system [15–17] with a large superconducting transition temperature (for a BCS superconductor).

This material presents three bands crossing the Fermi level, corresponding to the three Fermi sheets shown in figure 6. All these surfaces are visualized for the first BZ (blue) and the conventional cell or the polygon enclosed by the three reciprocal space lattice vectors $\left( {{{\mathbf{b}}}_{1}},{{{\mathbf{b}}}_{2}},{{{\mathbf{b}}}_{3}} \right)$ (red). The first band (${{\sigma }_{1}}$) generates a single cylindroid FS (top). The second band (middle) produces a FS which is in turn composed by two sheets: one of them (${{\sigma }_{2}}$) produces another cylindroid shape FS with a larger radius than that found in ${{\sigma }_{1}}$, and the second sheet (${{\pi }_{1}}$) corresponds to a ring shape FS inside the first BZ. The third band generates a single sheet (${{\pi }_{2}}$) which appears as a double ring FS in the first BZ, but translation to the conventional cell shows clearly that this surface is connected in a single sheet.

Zoom In Zoom Out Reset image size

From the point of view of the HFSH analysis hex-MgB$_{2}$ is a very interesting test example because there are several features in this material that could potentially appear in more complex systems: in hex-MgB$_{2}$ we have several bands (three) crossing the Fermi level, some of them even composed by multiple sheets. Thus, this system allows us to demonstrate that the HFSH method is also applicable in a complex multiple Fermi sheet environment following exactly the same procedure as in simpler single-sheet Fermi surfaces, as those found in bcc-Li and fcc-Cu, for example.

The area, the Euler characteristic (χ) and genus (g) for all the Fermi sheets of MgB$_{2}$ are shown in table 1. We emphasize again that in what concerns the DOS and the areas of the different sheets of the FS, the coarser tetrahedral sampling (20$^{3}$) produces practically converged results. Moreover, we stress again that direct integration of (12) (third column) compares very well with the linear tetrahedron method [8], as implemented in the Quantum Espresso code [14], which makes us confident about the quality of the scalar products needed for the HFSH mode analysis (see sections 3.7 and 3.8.1).

The ${{\sigma }_{1}}$ sheet of the FS (top of figure 6) is composed by a single sheet. Periodic boundary conditions impose that any point in the surface ${\mathbf{k}}$ is equivalent to ${\mathbf{k}}+{{{\mathbf{b}}}_{3}}$, thus this surface sheet is exactly a toroid from the topological point of view. This is confirmed with the integer values of the Euler characteristic and genus obtained ($\chi =0$, g = 1) up to numerical precision (${{10}^{-13}}$) for this band. The Euler characteristic is additive for surfaces which are not connected, as those found in the FS corresponding to the second band of MgB$_{2}$(middle in figure 6). In this surface we find that $\chi =0$ (genus g=1) for the central cylindroid (${{\sigma }_{2}}$) and Euler characteristic $\chi =-2$ (genus g = 2) for the outermost ring shape Fermi sheet described within the first BZ (${{\pi }_{1}}$). The third band of MgB$_{2}$ also shows a genus g = 2 which is better understood when one looks at the FS plotted inside the conventional cell (red polyhedra in figure 6). Periodic boundary conditions apply equally well for the conventional cell, so that each neck (there are four) is connected with another by a simple lattice translation. Thus one concludes that the genus is equal to g = 2. A similar reasoning applies for the ${{\pi }_{1}}$ sheet of the second band. Although topology is not our primary interest in this work, these results demonstrate the accuracy reached in the determination of the input data for the HFSH set in (9) and (10).

In figure 7 we show the calculated HFSH set for the three bands crossing the Fermi level for hex-MgB$_{2}$. We show the 16 lowest energy states for each Fermi sheet: first band (top), second band (middle) and third band (bottom). For each band the mode with the lowest energy is that enclosed by the BZ. Then the energy of the modes increases from right to left within the rows. This way, the highest energy mode sits, for each band, at the upper left-hand corner.

Zoom In Zoom Out Reset image size

The HFSH modes found for the first band (top) correspond basically to the solutions of the Helmholtz or stationary Schrödinger equations in a torus except for the weighting factor introduced by the local electron velocity, the deformation and the curvature. As the radius is smaller than the length of the cylindroid the first two non-trivial HFSH (second and third modes in the first row) are varying only in the axial direction while the third non-trivial HFSH mode (bottom left-hand corner) is the first one with transversal variation.

As we have mentioned before, the two sheets of the second band (middle of figure 6) are obviously not connected, thus the discretized version of the Helmholtz equation (9), appears as a block diagonal system. One strategy to solve the equation would be to separately diagonalize a version of (9) for each of these two sheets. However, at this point it is more interesting to treat the system barely and ignoring the block diagonal structure because in many complex systems we may find a situation where an obvious method for separating different sheets may not exist. We observe that in this band (middle of figure 7), the lowest energy HFSH mode (inside the first BZ) corresponds to zero value for the ${{\pi }_{1}}$ Fermi sheet (out ring shape) and a finite constant value for the ${{\sigma }_{2}}$ sheet (cylindroid). The next HFSH mode in energy is also trivial but it is finite and constant for ${{\pi }_{1}}$ and equal to zero for ${{\sigma }_{2}}$. We would obtain the same result if we separately diagonalized each diagonal block of the generalized eigenvalue equations (9) for the HFSH or (10) for BHFSH) and we energetically ordered the summation of both subspaces. Thus, the present example for the second band enables us to demonstrate that the brute force diagonalization works equally well for systems in which a simple inspection is not enough for the separation of the FS in different sheets. Or for those systems where the block diagonalization is not justified for the sake of computational simplicity.

Table 3 shows a summary of the HFSH mode analysis for the modulus of the electron velocity ($|v|$) the x cartesian component of the velocity (${{v}_{x}}$) and the multiplication of the x and y components of the velocity (${{v}_{x}}{{v}_{y}}$). The results are presented in separate columns for each Fermi sheet. We conclude that about 600 HFSH modes are sufficient for a reasonable $(\epsilon \lesssim 5\%)$ representation of the velocities and even of the tensor ${{v}_{x}}{{v}_{y}}$. Already with ∼2000 modes the method is able to capture the fine details $(\epsilon \lesssim 1\%)$ of these testing example functions $|v|$, ${{v}_{x}}$ and ${{v}_{x}}{{v}_{y}}$ for all sheets. Going to each band in detail, we observe that in the first sheet 600 modes are more than sufficient even for a fine detailed description. A tetrahedral sampling of density 40$^{3}$produced 3864 vertices in this sheet, so the saving factor for a target error of $(\epsilon \lesssim 1\%)$ would be of about ${{\alpha }_{1}}\left( {{40}^{3}} \right)\simeq 1/13$ . The Fermi surfaces corresponding to the second and third bands produced ${{N}_{v}}$=12558 and ${{N}_{v}}$=11478 vertices, respectively, in the triangulation process. Following a similar analysis, one obtains that the saving factor would be of the order of ${{\alpha }_{2}}\left( {{40}^{3}} \right)\simeq 1/6$ and ${{\alpha }_{3}}\left( {{40}^{3}} \right)\simeq 1/9$ for the second and third band, respectively.

Table 3.  Percentage of error obtained using (29) for ${{v}_{x}}$, ${{v}_{x}}{{v}_{y}}$ and modulus of the velocity ($|v|$) as a function of the number of modes used in the summation for each of the Fermi crossing bands in MgB$_{2}$. The dash means the error is lower than 0.1%.

	Band 1			Band 2			Band 3
${{N}_{L}}$	${{v}_{x}}$	${{v}_{x}}{{v}_{y}}$	$|v|$	${{v}_{x}}$	${{v}_{x}}{{v}_{y}}$	$|v|$	${{v}_{x}}$	${{v}_{x}}{{v}_{y}}$	$|v|$
101	3.5	3.8	1.6	15.9	27.0	6.4	8.3	12.8	2.8
601	0.2	0.5	0.1	1.9	5.0	0.6	1.5	2.1	0.4
1201	0.1	0.3	—	1.2	3.1	0.4	0.8	1.4	0.2
1801	—	0.2	—	0.9	2.1	0.3	0.2	0.5	0.1
2001	—	—	—	0.3	1.2	—	0.1	0.3	—

As mentioned before, the Bloch functions in (B.1) are maximally flat in reciprocal space -by construction in the Wannier procedure-, thus the Fourier interpolation scheme for any quantity involving the electron wave functions is very effective. The electron band energies and velocities are very efficiently calculated by considering the electronic DFT hamiltonian in the basis of Wannier Bloch functions ${{\psi }_{{\mathbf{k}},\alpha }}$

$$\begin{eqnarray}&&H_{{\mathbf{k}},\alpha ,\beta }^{W}=\mathop{\sum }\limits_{{\mathbf{R}}}{{{\mbox{e}}}^{{\mbox{i}}{\mathbf{kR}}}}\left\langle {{\chi }_{\alpha }}\left( {\mathbf{r}} \right)|{{H}_{{\mbox{DFT}}}}|{{\chi }_{\beta }}\left( {\mathbf{r}}-{\mathbf{R}} \right) \right\rangle .\end{eqnarray} \tag{ B.4 }$$

Note that in (B.4) ${\mathbf{k}}$ is not restricted to the self consistent DFT coarse mesh and any ${\mathbf{k}}$ vector is accessible, making the Fourier interpolation procedure straightforward. Considering any ${\mathbf{k}}$ wave vector in reciprocal space $H_{{\mathbf{k}},\alpha ,\beta }^{W}$ is accessible and the Fourier interpolated energies (${{\epsilon }_{{\mathbf{k}},i}}$) are obtained by an ordinary diagonalization procedure,

$$\begin{eqnarray}&&\mathop{\sum }\limits_{\alpha ,\beta }U_{{\mathbf{k}},i,\alpha }^{-1}H_{{\mathbf{k}},\alpha ,\beta }^{W}{{U}_{{\mathbf{k}},\beta ,j}}={{\delta }_{i,j}}{{\epsilon }_{{\mathbf{k}},i}}.\end{eqnarray} \tag{ B.5 }$$

The velocities are similarly calculated by considering the expectation value of the gradient of the Wannier hamiltonian and the solution eigenvectors in (B.5),

$$\begin{eqnarray}&&{{\vec{v}}_{{\mathbf{k}},i}}=\mathop{\sum }\limits_{\alpha ,\beta }U_{{\mathbf{k}},i,\alpha }^{-1}\left( \mathop{\sum }\limits_{{\mathbf{R}}}{\mbox{i}}\vec{R}\ {{{\mbox{e}}}^{{\mbox{i}}{\mathbf{kR}}}}\left\langle {{\chi }_{\alpha ,{\mathbf{0}}}}|{{H}_{{\mbox{DFT}}}}|{{\chi }_{\beta ,{\mathbf{R}}}} \right\rangle \right){{U}_{{\mathbf{k}},i,\beta }}.\end{eqnarray} \tag{ B.6 }$$