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

We present a new efficient numerical approach for representing anisotropic physical quantities and/or matrix elements defined on the Fermi surface of metallic materials. The method introduces a set of numerically calculated generalized orthonormal functions which are the solutions of the Helmholtz equation defined on the Fermi surface. Noteworthy, many properties of our proposed basis set are also shared by the Fermi Surface Harmonics (FSH) introduced by Philip B. Allen [Physical Review B 13, 1416 (1976)], proposed to be constructed as polynomials of the cartesian components of the electronic velocity. The main motivation of both approaches is identical, to handle anisotropic problems efficiently. However, in our approach the basis set is defined as the eigenfunctions of a differential operator and several desirable properties are introduced by construction. The method demonstrates very robust in handling problems with any crystal structure or topology of the Fermi surface, and the periodicity of the reciprocal space is treated as a boundary condition for our Helmholtz equation. We illustrate the method by analysing the free-electron-like Lithium (Li), Sodium (Na), Copper (Cu), Lead(Pb), Tungsten (W) and Magnesium diboride (MgB2).


Introduction
The Fermi surface (FS) of a metal is a characteristic property of the crystal structure and the material itself. It is found in many physical problems that the exact shape, topology and the precise value of a scalar, vector or tensor variable defined on this surface is crucial for a good description of any physical model. The importance of an accurate and efficient method for storing, filtering or interpolating a scalar or vector physical quantity on a Fermi surface is demonstrated very clearly in the prototype examples of transport and in the electron-phonon (EP) problem. In the EP interaction, the vibrational frequencies are typically about one/two orders of magnitude smaller than the electron energies, and to a good approximation, the interaction is described by scattering events connecting different points of the Fermi surface k and k . The EP theory requires the knowledge of a huge amount of electron-phonon matrix elements g ν k,k connecting different points of the Fermi surface mediated by several phonon branches (ν). For example, one needs typically about n k ∼10 4 data points for a good resolution of the matrix elements defined on the Fermi surface. In a rough estimation for a simple EP problem, considering n p ∼10 phonon branches and about n e ∼10 electron bands, the required amount of data for the matrix elements is then N ∼ 10 10 . Thus, filtering or compressing all this data while keeping accuracy appears very appealing.
The essence of the present method is to perform a linear integral transformation on quantities such as the EP matrix elements g ν k,k →g ν L,L so that the dimension of the new basis set is smaller than the original one (n L ∼10 2 instead of n k ∼10 4 ).
The idea of representing scalar quantities as a function of an orthogonal set defined on the Fermi surface was introduced by Philip B. Allen in 1976 with the so-called Fermi Surface harmonics (FSH) [1,2]. Allen considered a functional set constructed as polynomials of the cartesian components of the electronic velocity, and designed -in principle-to operate with any crystal structure and number of Fermi sheets. The author was able to rewrite the electron self energy, transport equations and the Eliashberg theory of the superconductivity in terms of the FSH set. Although the great potential of the FSH method by Allen appeared promising at first, the above method has not found a systematic application yet, being applied only in relatively simple systems [3]. The weak impact of the above method might be attributed to the several semi-analytic steps involved, the relatively complex treatment of different Fermi sheets, the difficulty to generate the functional set and to the fact that the completeness of the latter was not guaranteed.
We propose a new functional set with a similar spirit and motivation as in Allen's FSH [1,2] but defined very differently, and constructed in such a way that the connection with ordinary Fourier transform in flat space and/or with ordinary spherical harmonics functions in a (curved) sphere is direct. Our proposed functional basis is defined to satisfy a modified version of the Helmholtz equation defined on the Fermi surface. More graphically, in our method the Fermi surface is considered as if it were a vibrating membrane. The standing waves calculated on this surface constitute the new proposed functional set. As these functions are the solution of a second order partial differential equation, the -generalized-orthogonal property of the basis set is recovered by construction, and more importantly, the completeness of the set is automatically guaranteed. Another advantage of the present method, in comparison with the Allen FSH set, is that a definition of an energy cutoff E c for the basis set appears naturally. As in ordinary plane wave theory, the plane wave cutoff (E c ) allows an estimation of the minimal size which is describable by the new basis set. Thus, the sharper the target details defined on the Fermi surface, the larger the cutoff and the number of modes we need for its accurate description. Furthermore, the global and topological properties of the surface are included automatically as demonstrated, for example, by considering a direct application of the Gauss-Bonnet theorem which we utilize as a check.

Mathematical definition and implementation of Helmholtz Fermi Surface harmonics (HFSH)
We define de Bare Helmholtz Fermi Surface harmonics functional basis (BHFSH), as the eigenfunction set of the Laplace-Beltrami operator defined on the Fermi surface and with crystal periodic boundary conditions We stress that ∇ 2 k operates on a 2D surface, not in the embedding 3D k-space, and it must be interpreted as in the ordinary Laplace operator, i.e., as the divergence of gradient. From now on, we will refer to the Laplace-Beltrami operator simply as the Laplace operator. In (1) κ 2 L are the eigenvalues associated to the BHFSH, which satisfy the following orthogonality relation: In many mathematical problems such as when the unique objective is to simply compress a large numerical data set defined on the Fermi surface, the set {Ψ L (k)} might result convenient. However, in many physical theories, the integrals over the Fermi surface contain a weighting factor which is proportional to the local density of sates, i.e. the inverse of the electron velocity. It may be thus desirable to introduce a generalized orthogonality condition incorporating the inverse of the electron velocity as a weighting factor. Therefore, we define the Helmholtz Fermi Surface harmonics (HFSH), as the eigenmodes of a modified version of (1), where we introduce the absolute value of the electron velocity in a given point k of the Fermi surface, v(k): where ω L are the eigenvalues associated to the HFSH set {Φ L (k}.
As the HFSH elements are the solutions of a second order -linear-partial differential equation, the generalized orthogonality property is now which includes naturally the local density of states as a weighting factor. Clearly, the simple plane wave basis set (e ik·r ) and the ordinary spherical harmonics basis set [Y m l (k)] are automatically recovered in the HFSH set, defined as above, because these functions satisfy the Helmholtz equation in their respective spaces, the flat space and the curved unit sphere, respectively. Thus, our approach applied in an approximately spherical Fermi surface with nearly constant density of states leads naturally to the ordinary spherical harmonics.
The Fermi surface lives in the electron momentum space which is periodic, thus in our problem (3) must be solved with the appropriate periodic boundary conditions. We will see that in the simple case of a single Fermi sheet with no intersection with the boundaries of the first Brillouin zone the periodic boundary conditions are automatically satisfied. However, for surfaces or multiple Fermi sheets intersecting the Brillouin zone boundaries, a numerically more elaborated procedure is needed.

Numerical scheme
In this section we describe the numerical algorithm to solve (3), or alternatively (1).
A brief description of the procedure could be the following. In a first step, a triangulated mesh of the Fermi surface is constructed, such that the discretized version of our generalized Helmholtz equation given by (3) can be solved numerically. Once the triangular tesselation of the surface is obtained, a discrete version of the Laplace operator can be derived, as described for instance in [4,5], and in this way, the original partial differential equation is transformed into a generalized sparse eigenvalue problem. The solutions of this linear problem are exactly our proposed HFSH set.

Discretization of the problem: Triangular tesselation of the Fermi surface
As already mentioned, the construction of the triangular mesh over the Fermi surface 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.

Discrete version of the Helmholtz equation in a triangulated surface
The triangulation of the Fermi surface -with any method-produces a set of n t triangles whose vertices v i (i=1,. . . , n v ) are 3D k-space vectors obtained by the triangulation Figure 1. Schematic representation of the triangulation of the Fermi surface. Each vertex, v i , is a 3D k-space vector lying on the Fermi surface and is surrounded by its nearest neighbour vertices v j . The shaded area represents the barycenter control area around vertex v i . The angles entering the discretized version of the Laplace operator in (7), α i,j and β i,j , are the opposite angles of the the edge joining vertices v i and v j shared by two adjacent triangles.
process and n v are the number of vertices. Figure 1 illustrates the situation. Each vertex (v i ) has N n (v i ) nearest neighbour vertices (v j ) and nearest neighbour triangles, denoted by T j (v i ), j = 1, . . . , N n (v i ). The shaded area around each vertex v i defines the "barycenter control area", which is the area associated to each vertex (or 3D kspace vector on the Fermi surface). The control area (S i ) of a vertex 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 1 3 of each neighbouring triangle area, where A[T j (v i )] denotes the area of the nearest neighbour triangles T j (v i ) of vertex 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 Fermi surface, 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. 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 (v j ) and triangles of a given vertex 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, α i,j and β i,j , of the triangles sharing the edge joining the vertices v i and v j (see figure 1). Thus, the two point centered formula for the second derivative in one dimension is generalized by [4,5] where 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: Note that if the area of the control cells of all vertices (S i ) were equal and the velocity [v(v i )] was constant in (9), we would then have a regular eigenvalue problem for the eigenfunctions Φ L and the eigenvalues ω L which are labelled by L. The same applies to (10). Moreover, since the Ω i,j matrix is symmetric and the S i v(v i ) δ 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].

Periodic boundary conditions
As mentioned, in our computational scheme a first step consists on 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 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 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 k + 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.

Fermi surface relaxation: Refinement of the triangular mesh by application of the Newton-Raphson method.
Within the marching tetrahedra algorithm (Appendix A), the entire Brillouin zone 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 Fermi surface. 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 Fermi surface 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 k vectors are allowed to move along the normal to the Fermi surface such that the relaxed vector -and up to numerical precision-lies exactly at the Fermi surface. 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, In (11) the electron band indices are not shown for simplicity, but this relaxation scheme is applied to all k triangle vertices defining the Fermi surface and all bands composing the Fermi surface. 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 to 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.

Density of states (DOS)
Obviously, the refinement of the mesh improves the estimation of Fermi surface 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 density of states (DOS) ρ(E), where d 2 s 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 (v i )] for all the relaxed vertices i and bands n, 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 k point vertices and (ii) the triangular mesh is iteratively improved by forcing the triangle vertices to lie at the Fermi surface. Thus, the above method for estimating Fermi integrals may be considered as a surface specialized non-linear version of the linear tetrahedron method.

Topological characterization: Euler characteristic and the Gauss-Bonnet theorem applied to Fermi surfaces
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 Fermi surface 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(k) being the Gaussian curvature at point k of the surface, the Euler characteristic χ is given by 1 2π 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 Brillouin zone.
Numerically, the Gaussian curvature is given by where τ i,j denotes the angle between the triangle edges joining vertices v i and v j , and v i and v j+1 (see figure 1) and S i is defined in (5). Thus, one obtains that numerically (14) can be rewritten as which is the exact Descartes theorem on total angular defect of a polyhedron. The genus of a single-sheet surface is given by While χ is additive for multiple-sheet surfaces, g must be considered separately for each sheet. It is worth to mention 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 Fermi surfaces 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 χ = −6 (and g = 4).

General properties of the HFSH set
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 Fermi surface [(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 lines we face the problem of function representation in terms of HFSH.
The HFSH set allows us to efficiently represent any function defined on the Fermi surface, but more important, the HFSH modes allow 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 Fermi surface are for instance the electron lifetime, τ (k), the electron mass enhancement, λ(k), the superconducting gap, Λ(k), electron velocity components, etc. Similar to any integral transformation, the smoother the function to be represented the more efficient is 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 (k), and let us rewrite this function in terms of the HFSH set, We refer to (18) as the HFSH representation of the k dependent function F (k) defined on the Fermi surface. The scalar product and normalization of the HFSH {Φ L } are defined by (4) in momentum space, thus the components c L (F ) of the expansion become directly accessible, the coefficients of the expansion having the same units as the original F (k) function, 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, . (20) The procedure of transforming a function F (k) into a discrete set of coefficients c L (F ) 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 the control on the desired accuracy by tuning the energy cutoff.
The product of two functions F 1 (k)F 2 (k) defined in k space may be represented as a function of the HFSH components of each of these functions separately [c L (F 1 ) and c L (F 2 )] with the following relation that generalizes the convolution theorem in Fourier analysis where the matrix elements Ξ L;L 1 ,L 2 play the same role as the Clebsch-Gordan coefficients for ordinary spherical harmonics, The coefficients Ξ L;L 1 ,L 2 are symmetric with respect to the permutations of the L, L 1 , L 2 indices, as it was also found by Allen for the FSH set [1]. Some elementary properties of Ξ L;L 1 ,L 2 are (the first element of the HFSH set is L = 1) The Ξ coefficients allow to expand the product of two Φ functions in terms of a simple linear combination of HFSH elements and vice-versa, Helmholtz Fermi Surface Harmonics: an efficient approach for treating anisotropic problems involving Fermi The scalar product of pairs of HFSH may be expressed also in terms of the Ξ matrix elements, In principle, one would calculate all the Ξ 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 electron-phonon matrix elements g k,k 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 k and k . In all these cases we would follow a similar procedure, where Of course the function g k,k 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 k,k . 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,11,12]. Alternatively, one could consider the absolute values of the matrix elements |g k,k | 2 . The above algebraic machinery has a great potential in restating integro-differential problems defined on the Fermi surface. 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.

Numerical tabulation of an anisotropic quantity
One of the main problems encountered when trying to characterize an anisotropic physical quantity defined on the Fermi surface, 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 k vectors defining the FS. However, this method is mainly visual and not quantitative.
The representation of a function F (k) using the HFSH expansion coefficients c L (F ) enables a direct quantitatively description of the anisotropy. Indeed, the HFSH method allows to tabulate numerically any complex anysotropic quantity by means of about 10 2 coefficients, the most representative being those corresponding to the modes with the lowest energy.
Since Φ L=1 (F ) = 1, the first expansion coefficient, C L=1 (F ), yields directly the FS average of the function F (k). In materials with spherical or nearly spherical symmetry, such as Li or Cu (see section 4.1), the next three modes [Φ L (F ), 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 (∼ 20) of coefficients is sufficient for capturing the most significant part of the anisotropy of any 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 electron-phonon mass enhancement, λ k , or the momentum dependent lifetime, τ k . Following a standard procedure, if one wanted to compare two different calculations of λ k , for instance, obtained using two different computation methods, one would need to compare the 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 (F )] 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 Fermi surface 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 (k) explicitly for all the vertex k points of our triangular grid describing the FS. Consider a finite cutoff (E c ) for the expansion of F (k) in terms of the HFSH set using (18): where N L denotes the number of modes such that ω L < E c . We then obtain a smoothed functionF N L (k) where the fine details smaller than a wave length are filtered.
A measure of the mismatch error when considering only a finite set of HFSH modes, N L might be obtained by the Fermi surface integral which is approximated by a simple linear quadrature formula as in (6).

Test examples for real materials
We have considered materials with one Fermi sheet (bcc-Li, bcc-Na and fcc-Cu), with two Fermi sheets (fcc-Pb and bcc-W) and with three Fermi sheets (hex-MgB 2 ) as testing examples.
The DFT ground state for all of them have been obtained using norm-conserving pseudopotentials with the PBE [14] functional and with a cutoff energy of 30 Ry for bcc-Li, 50 Ry for bcc-Na and fcc-Pb, 55 Ry for bcc-W, 60 Ry for hex-MgB 2 , and 110 Ry for fcc-Cu. Each DFT ground state was in turn used as an input for the Wannier calculations which were carefully converged.  The second column shows the momentum sampling for the tetrahedral division in the marching tetrahedron method for constructing the triangulation of the Fermi surface. The third column displays the DOS obtained using the method presented in this work [using (12)], which includes a relaxation of the surface, as introduced in section 3.4, and an explicitly calculated velocity for each vertex position of the triangulation using the Wannier method, as explained in Appendix B. These results are compared to the linear tetrahedron method [8] (last column), where both, the vertex (without relaxation) and the magnitude of the velocity are included effectively by linear direct interpolation. The agreement between both values is remarkable in all cases and allows us to check that the quadrature formula for the integrals is accurate when a triangular division together with a barycenter control cell (see figure 1) is considered. The fourth column shows the calculated area of the different Fermi sheets for each momentum space sampling. The next two columns show, for each Fermi sheet, the calculated Euler characteristic (χ) and genus (g = 1 − χ/2 for a connected surface). As a convergence test, all these magnitudes were evaluated for the two different sampling densities, 20 3 and 40 3 . Note that all the quantities are found to be practically converged already for a sampling density of 20 3 .
Bcc-Li and bcc-Na present a Fermi surface which is topologically trivial and this is confirmed by a value of the Euler characteristic of χ=2 (g=0) in both cases. The topological classification of the Fermi surfaces is not the goal of this work, but since internal angles of the triangles simplexes enter (7) and (8) in a way which is crucial, a direct computation of the Gaussian curvature, considering the total angular defect in (15), is an essential test of consistency. Regardless of the choice of the unit cell, the tetrahedral sampling, or the complexity of the surface, we obtain in all cases that the Euler characteristic (χ) is recovered as an integer number up to double real numerical precision (∼ 10 −13 ). As an application of a more complex example, we mention the Euler characteristic of the first band of bcc-W [χ = 2 × (6 + 1)], which can be easily understood by inspecting its Fermi surface (not shown). It is composed of 6 disjoined ellipsoidal sheets around the high symmetry point N , and one diamond-shaped sheet centered at H, each of them yielding a genus g = 0.
As for the visualization of the HFSH, we have considered only bcc-Li, fcc-Cu and hex-MgB 2 as representative examples. Bcc-Lithium and fcc-Cu are found to be reasonably close to the ordinary spherical harmonics, although, the Fermi surface of fcc-Cu is topologically not so trivial as that of bcc-Li. The hex-MgB 2 Fermi surface presents three different electron bands crossing the Fermi level, with one of these bands composed by two disconnected pieces. MgB 2 enables us to demonstrate the utility of the method in a more complex situation, but on the same footing as in simpler Fermi surfaces.

Free electron like bcc-Li and fcc-Cu
Our first example applications are bcc-Li and fcc-Cu at ambient pressure, which are textbook examples of free-electron-like systems. The Fermi surface 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 2a and 2b show the calculated Fermi surface 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 | k − E F | < 10 −6 eV. The colour code in figures 2a and 2b 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 Fermi surface, we observe that the anisotropy of the electron velocity is not negligible.

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 figure 3 and  figure 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  (3) is solved analytically in a spherical surface of radius k F and average velocity v F , ω = l(l + 1)v F /k 2 F . The lowest energy state (first row) is a constant function Φ L=1 =1, according to the -defining-normalization condition in (4). In the second row we find three degenerate HFSH modes, approximately resembling the ordinary p x , p y , p z spherical harmonics. In the third and fourth rows we show the calculated HFSH modes which are similar to the d and f states, but in these case the degeneracy is lifted as perfect spherical symmetry is absent.
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 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 Fermi surface of Cu represent a strong perturbation to the sphere, mostly noticeable at long wavelengths. 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 [Ω 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 → ∞ and Equations (31) and (32) work very well already for L ∼ 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.   (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 Fermi surface. 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.
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 Fermi surface 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 (v x ) = 0, while the next two modes c L=2, 3 (v x ) are the most important. This is consistent with the interpretation of the first three nontrivial 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 Fermi surface of the v x and v x v y functions are zero because these functions are obviously odd, and therefore c L=1 (v x ) = 0 and c L=1 (v x v y ) = 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 ∼ 180 the HFSH amplitudes (c L ) are already about 20 times weaker than the maximum values of c L at lower L. 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 ).  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 ∼701 modes is sufficient in order to obtain an estimated error below 1% with a tetrahedral sampling of 40 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 20 3 v = 6018 and N 40 3 v = 24582, respectively. One can estimate the data-storage saving when using the HFSH-mode representation instead of the usual vertex (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 α(40 3 )∼1/40 would be obtained with a 40 3 tetrahedral mesh.
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 Brillouin zone (blue) and the conventional cell or the polygon enclosed by the three reciprocal space lattice vectors (b 1 , b 2 , b 3 ) (red). The first band (σ 1 ) generates a single cylindroid Fermi surface (top). The second band (middle) produces a Fermi surface which is in turn composed by two sheets: one of them (σ 2 ) produces another cylindroid shape Fermi surface with a larger radius than that found in σ 1 , and the second sheet (π 1 ) corresponds to a ring shape Fermi surface inside the first Brillouin zone. The third band generates a single sheet (π 2 ) which appears as a double ring Fermi surface in the first Brillouin zone, but translation to the conventional cell shows clearly that this surface is connected in a single sheet.
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 Fermi surface, 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 [13], 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 σ 1 sheet of the Fermi surface (top of figure 6) is composed by a single sheet. Periodic boundary conditions impose that any point in the surface k is equivalent to k + 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 (χ = 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 Fermi surface corresponding to the second band of MgB 2 (middle in figure 6). In this surface we find that χ = 0 (genus g = 1) for the central cylindroid (σ 2 ) and Euler characteristic χ = −2 (genus g = 2) for the outermost ring shape Fermi sheet described within the  first Brillouin zone (π 1 ). The third band of MgB 2 shows also a genus g = 2 which is better understood when one looks at the Fermi surface 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 π 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.
The HFSH modes found for the first band (top) correspond basically to the solutions of the Helmholtz or stationary Shrö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 π 1 Fermi sheet (out ring shape) and a finite constant value for the σ 2 sheet (cylindroid). The next HFSH mode in energy is also trivial but it is finite and constant for π 1 and equal to zero for σ 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 Fermi surface 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 ( 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 ( 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 ( 1%) would be of about α 1 (40 3 ) 1/13 . The Fermi surfaces corresponding to the second and third bands 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%. 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 α 2 (40 3 ) 1/6 and α 3 (40 3 ) 1/9 for the second and third band, respectively.

Conclusions
We propose a new functional set, the Helmholtz Fermi Surface Harmonics (HFSH), which shows very interesting properties for efficiently representing physical quantities and/or integro-differential equations defined on the Fermi surface. This functional set is defined as the solutions of a Helmholtz type equation defined on top the Fermi surface, and we describe in detail the numerical scheme needed to solve this equation in a curved space including periodic boundary conditions. Although the HFSH presented in this work are defined very differently from the FSH proposed by Allen, all the analytical results for the Boltzmann equation and for the anisotropic Eliashberg equation reported in [1,2] are still valid for the HFSH presented in this work. Despite the fact that topology is not the main goal of this work, direct application of the Gauss-Bonnet theorem (integral of the Gaussian curvature) has showed that the genus and the Euler characteristic of the Fermi surface is easily accessible only with the input data needed to set up the Helmholtz equation on the Fermi surface. Thus, a systematic/automatic study of the topology of the Fermi surface of a material, for example as function of pressure, is accessible with this method and without the need of any visual interpretation.
We have applied our method in bcc-Li, bcc-Na, fcc-Cu, fcc-Pb, bcc-W and hex-MgB 2 , demonstrating that a systematic procedure is easily implemented and that the method is robust, even in systems with complex band structures and/or several Fermi sheets.
We have expanded the cartesian components of the electron velocity and their product in terms of HFSH, and we have found that a relatively small number (∼10 3 ) of HFSH modes is enough for a very accurate description (error 1%). Thus, the HFSH mode representation enables to numerically tabulate any anisotropic magnitude defined on the Fermi surface, allowing not only qualitative and faithful comparisons of Figure A1. Tetrahedral subdivision of a cube in terms of cube vertices. We consider 6 tetrahedra in each cube with labeling corner (i, j, k).) these magnitudes calculated from different computational methods, but also allowing the efficient integration of any anisotropic function over the Fermi surface. The triangulation of the Fermi surface is calculated by checking, in a first step, that the Fermi level (e F ) is located between the minimum and the maximum of the energies corresponding to all the tetrahedral vertices. That is, if Min(e i ) < e F < Max(e i ), i = 1, 4 part of the Fermi surface is found in the interior of a given tetrahedra. Let us suppose for simplicity, that the energies are ordered in increasing order e 1 ≤e 2 ≤e 3 ≤e 4 . In this case, one can only find two nontrivial possibilities, (i) corresponding to the simple triangular section when e 1 ≤e F ≤e 2 ≤e 3 ≤e 4 or e 1 ≤e 2 ≤e 3 ≤e F ≤e 4 (left side of figure A2), and (ii) the case when the intersection is quadrilateral e 1 ≤e 2 ≤e F ≤e 3 ≤e 4 (right side figure A2), and the two triangle simplexes are generated. There are two possibilities to divide a quadrilateral by triangles, among them, our choice is the one in which both triangles have an area as similar as possible.
In all cases the vertices of the intersection are easily calculated by considering that the electron energy is a linear function inside each tetrahedra. The vertices of the -Fermi surface-intersection triangle corresponding to the case 1 ( e 1 ≤e F ≤e 2 ≤e 3 ≤e 4 ) are given by the vectors The case e 1 ≤e 2 ≤e 3 ≤e F ≤e 4 is equivalent to the one above by considering the opposite -descending-energy ordering. The vectors describing the quadrilateral intersection in case 2 (e 1 ≤e 2 ≤e F ≤e 3 ≤e 4 ) are similarly obtained, t 5 = t 1 + e F − e 1 e 3 − e 1 (t 3 − t 1 ) t 7 = t 2 + e F − e 1 e 4 − e 2 (t 4 − t 2 ) t 8 = t 2 + e F − e 1 e 3 − e 2 (t 3 − t 2 ) , and the triangle vertices would correspond to (5,8,6) and (6,8,7), or alternatively to (5,7,6) and (5,8,7). Our choice is the one in which both areas a maximally similar.
Appendix B. The maximally localized Wannier approach for calculating electron energies and velocities The Wannier method is used extensively in the present work as an efficient input tool for calculating the electron energy and velocities. In this section we introduce this method only succinctly for completeness, and refer the reader to the original works [10,18,19] for more details.
In the Wannier scheme as considered in this article, the first step is to perform an ordinary DFT electronic structure calculation on a regular and relatively coarse Monkhorst-Pack k mesh of density n 1 × n 2 × n 3 in order to obtain the energies k,i and eigenfunctions Φ k,i (r) for the given BZ division and band index i.
The maximaly localized Wannier functions are calculated (and defined) by considering a unitary transformation (mixing) among Bloch states and minimizing a spread functional. Considering, for simplicity, n isolated bands in each k point one considers a general lattice periodic function as ψ k,α (r) ≡ β U k,α,β Φ k,β (r), (B .1) where U k,α,β represents -in principle-any unitary matrix of dimension n. In this method, however, one is interested in making the functions ψ k,α as maximally flat as possible, such that any interpolation procedure works out effectively. Since maximally flat in reciprocal space is equivalent to maximally localized in real space, the U α,β unitary matrices are fixed by minimizing the spread functional where the Wannier functions are defined by with the unitary matrices in (B.1) fixed by minimizing the spread Ω in (B.2). In the equation above the summation in k goes over the initial Monkhorst-Pack coarse mesh of density n 1 × n 2 × n 3 considered for the self consistent DFT calculation.