A Practical Method to Detect, Analyze, and Engineer Higher Order Van Hove Singularities in Multi‐band Hamiltonians

Abstract

A practical method is presented to detect, diagnose, and engineer higher order Van Hove singularities in multiband systems, with no restrictions on the number of bands or hopping terms. The method allows us to directly compute the Taylor expansion of the dispersion of any band at arbitrary points in momentum space, using a generalized extension of the Feynman–Hellmann theorem, which is stated and proved. Being fairly, in general scope, it also allows to incorporate and analyze the effect of tuning parameters on the low energy dispersions, which can greatly aid the engineering of higher order Van Hove singularities. A certain class of degenerate bands can be handled within this framework. The use of this method is demonstrated, by applying it to the Haldane model.

Introduction

Nearly a century after its introduction, modeling of the band structure of crystalline solids continues to play a significant role in our understanding of real materials. Much progress has been made, both in the realm of techniques for the calculation of band structures, and the analysis of non-trivial band features. Progress in the latter is particularly epitomized by the rapid advances made in recent times, in the areas of band topology and geometry. The investigation of unusual band topologies and their physical importance was initiated by Lifshitz and Van Hove.^[1,2] Their work led to the notions of Fermi surface topological transitions, where the Fermi surface geometry changes suddenly when some parameters in the system are tuned and the associated density of states (DOS) divergences. In his pioneering early work, Lishitz identified two types of Fermi surface topological transitions, namely the pocket appearing/disappearing type and the neck formation/collapse type.^[1] Usually in Fermi surface topological transitions, the Fermi surface hosts one or more critical points of the dispersion.

Critical points are accompanied by a vanishing gradient of the energy dispersion, and typically correspond to either a maximum, a minimum, or a saddle point of the dispersion. For a 2D system, the canonical form of the Taylor expansion of the dispersion at a critical point to quadratic order, is normally given by $$\pm k_x^2 \pm k_y^2$$. Such a form of the dispersion leads to a logarithmic divergence in the DOS^[2] at the critical energy. Sometimes, however, we need to expand the dispersion beyond quadratic order to adequately capture the low energy features. These are situations where the determinant of the Hessian matrix of the dispersion, evaluated at the critical point, also vanishes. In such a scenario, the critical point is referred to as a higher order critical point, and the dispersion is said to have a higher order Van Hove singularity (HOVHS). A simple and commonly reported example of HOVHS is the cusp singularity with dispersion $$k_x^4 - k_y^{2}$$.^[3–5]

In the case of HOVHS, the divergence of the DOS is enhanced to a much stronger power law, often with asymmetric prefactors above and below the critical energy.^[3,6] This is expected to strengthen electronic correlations in the vicinity of the HOVHS. The consequences of this could be potentially dramatic,^[7] and may lead to novel quantum phases. The fundamental reason is that due to the vanishing gradient of the energy dispersion in the vicinity of those points, the fermions get heavy and interaction effects become important.

Ordinary Van Hove singularities have been observed in the context of Lifshitz transitions in many materials including cuprates and iron based superconductors, cobaltates, Sr₂RuO_4, and heavy fermion systems.^[8–17] With the recent advances in the Experimental Section of tuning materials and improved resolutions of various spectroscopic probes (particularly, angle resolved photoemission spectroscopy (ARPES)), there is both an increasing interest and emerging avenues for the investigation of the presence and consequence of exotic band geometry in real materials. This is particularly true for HOVHS, which have been linked with the unusual magnetic and thermodynamic behavior of Sr₃Ru₂O₇,^[4] twisted bilayer graphene near magic angle, close to half filling,^[18–20] the recently proposed unusual quantum phase referred to as “supermetal”, which has diverging susceptibilities without the presence of any long range order,^[21] the distinct Landau level spectrum of biased bilayer graphene^[5] and the “extended” Van Hove singularity in both doped graphene^[22] and highly overdoped graphene.^[23] More recently, HOVHS have been reported in the kagome superconductors,^[24–27] and may be relevant to the observed phases of doped Bernal bilayer graphene.^[28]

It is important to note that while a regular VHS is simply an ordinary saddle point singularity (having canonical dispersion $$k_x^2 - k_y^2$$), HOVHS are, in contrast, a large class of distinct singularities. Although in principle there are infinitely many of them, we need to consider only a finite subset of singularities when studying 2D and layered materials (see, in baseline. Ref. [3] for a justification, and Refs. [3, 6] for a thorough classification of these relevant singularities). Apart from this, another complication that arises in the context of HOVHS is that they are more elusive compared to regular VHS, that is, they typically require a fine tuning of parameters in the system (by the application of strain, pressure, twist, bias voltage, doping, etc). In addition, HOVHS may be sensitive to impurities and changes in the parameters, upon which the higher order critical point could split into a number of ordinary critical points. As a recent study of the effect of impurities on HOVHS^[29] shows, for a low concentration of impurities, the distinct quantitative signatures of the HOVHS do survive.

As already noted, there has been a renewed interest and improved ability to experimentally tune real materials by various probes. This renders a sound theoretical analysis of HOVHS a highly relevant and necessary task. To this end, in this work, we address the question of how we can detect a HOVHS given a k-space Hamiltonian, even when the description is detailed (in the sense of number of terms, hopping parameters, etc.) and involves multiple bands. The HOVHS can be diagnosed using the Taylor expansion of the band in question, at the chosen k point. Low energy theories determined either by polynomial fits or purely symmetry considerations are often also a stepping stone for further analysis, such as diagrammatic many-body treatment of electronic correlations and transport studies.^[7,29,30] But the full and correct Taylor expansion to a given order is needed for the unambiguous diagnosis of a HOVHS (see Ref. [3] for an illustration of how relying purely on symmetry considerations can lead to a misdiagnosis of the HOVHS, the DOS divergence exponent and the ratio of prefactors).

While HOVHS have been increasingly reported in real materials and theoretical models, detecting their presence and diagnosing their precise nature is not a straightforward task. Even when the Taylor expansion of the band at the critical point is available, naive inspection or even scaling arguments may not help us identify the correct singularity, and its associated power law exponents and ratio of prefactors (see Ref[3] for a discussion of the issue with an illustrative example). Nevertheless, a straightforward algorithm exists, which, given the Taylor expansion, can be used to exactly pinpoint the singularity that is present.^[3] To this end it becomes necessary for us to compute the Taylor expansion of the bands in tight-binding and Wannierised models, particularly when these models are able to account for experiments on materials (such as angle resolved photoemission spectroscopy and quasiparticle interference measurements).

What is widely used in this context is polynomial fits to the band along some directions in k-space. While this may indeed reveal some aspects of the Taylor expansion, it may at times altogether miss the presence of a higher order singularity if the directions of the fits do not coincide with the “principal” directions of the singularity. For example, the singularity with dispersion relation: 1 $$$\begin{equation} \varepsilon ({\bf k}) \propto (k_x^{\,} + k_y^{\,})^4 - (k_x^{\,} - k_y^{\,})^2 \end{equation}$$$will not be revealed as a higher order saddle by polynomial fits of the dispersion along the $$k_x^{\,}$$ and $$k_y^{\,}$$ directions. Even in the cases where the polynomial fit reveals a higher order singularity, by virtue of absence or smallness of the quadratic part of the fit, the full Taylor expansion to a specified order may not be revealed by polynomial fits along some chosen directions. In fact, the fitted polynomial may even have very different coefficients compared to the Taylor expansion arising from a tight-binding model that adequately describes the same dispersion. Therefore, given a tight-binding model, we need a general technique for extracting the Taylor expansion of bands at arbitrary k-points unambiguously.

In this paper, we provide a technique based on an extension of the Feynman–Hellmann theorem, to directly compute the Taylor expansion to any order, of band dispersion given the k-space Hamiltonian, with no restriction on number of bands and hopping terms. The method can be used to treat a case of degenerate bands, although it can be extended to any case. An explicit proof of the method, to any degree of Taylor expansion is provided in the appendix. We first deal with multi-band models in Section 2 where we provide the method and some general formulae. It is worth emphasizing that the method can be applied to any Hamiltonian H(k), not just ones arising from tight binding models. In Section 3, we apply the method and the generic results to the case of the Haldane model tuned to a higher order singularity. Finally, we conclude in Section 4. For clarity and to facilitate the flow, details of the calculations and the full general proof are presented in the appendix.

Multi-band Tight Binding Models

As described in the introduction, the active bands near the Fermi level of many real materials can be captured by complicated tight binding and k · p models that might, in practice, require a large number of Wannier functions or orbitals and several hopping terms to be able to adequately reproduce DFT and ARPES bands. This can often make the analysis of low energy theory around some k point quite difficult, especially if we want to avoid resorting to polynomial fits or if we are interested in tracking the effect of changing tuning parameters on the dispersion. Since the characteristic equation has a degree equal to the dimension of the Hamiltonian matrix (i.e., the total number of bands), any analytic solution to the eigenvalue problem for general k, is possible only up to the four band case, since the general quintic polynomial (of degree 5) and higher order polynomials can not be solved in terms of radicals.^[31] Already for cubic and quartic polynomials, the general solutions are quite complicated. We, therefore, need a more practical method to obtain the Taylor expansion around a particular k point in a specified band. In the next section, we describe such a method that can handle large and detailed k-space Hamiltonians (including those based on TBMs) and can be readily implemented numerically.

Preliminaries

To compute the Taylor expansion to order N at $$\bm{k}_0^{\,}$$, in principle we need to compute all the mixed partial derivatives $$\partial _1^{l_1} \partial _2^{l_2} \cdots \partial _d^{l_d} \varepsilon _n (\bm{k})|_{\bm{k}_0}$$ for $$1 \leqslant l_1^{\,} + l_2^{\,} + \cdots + l_d^{\,} \leqslant N$$ in a d-dimensional system. However, computing all these mixed partial derivatives within the Feynman–Hellmann method will be rather involved. We can simplify this procedure considerably by computing instead, the derivatives with respect to a single variable λ, taking the form $$\partial ^M \varepsilon _n^{\,} (\bm{k}_0^{\,} + \lambda \, \bm{k}) / \partial \lambda ^M$$ and setting λ to zero at the end. By repeated application of the chain rule for differentiation, it can be shown that 2 $$$\begin{eqnarray} &&\frac{\partial ^M \varepsilon _n^{\,} (\bm{k}_0^{\,} + \lambda \, \bm{k})}{\partial \lambda ^M} \bigg |_{\lambda = 0} = \sum _{\substack{1 \leqslant l_1, l_2, \ldots , l_d < M \\ l_1 + l_2 + \cdots + l_d = M}} \frac{M!}{l_1^{\,}! \, l_2^{\,}! \cdots l_d^{\,}!} \, \frac{\partial ^M \varepsilon _n^{\,} (\bm{k})}{\partial k_1^{l_1} \, \partial k_2^{l_2} \cdots \partial k_d^{l_d}} \bigg |_{\bm{k} = \bm{k}_0} \nonumber\\ &&\quad k_1^{l_1} \, k_2^{l_2} \cdots k_d^{l_d} \end{eqnarray}$$$Thus, $$\frac{1}{M!} \frac{\partial ^M}{\partial \lambda ^M} \varepsilon _n^{\,} (\bm{k}_0^{\,} + \lambda \, \bm{k}) |_{\lambda = 0}$$ is the sum of all the monomials of degree M in the Taylor expansion of $$\varepsilon _n^{\,} (\bm{k})$$ at $$\bm{k}_0^{\,}$$. The full Taylor expansion to order N is then given by 3 $$$\begin{equation} \varepsilon _n^{\,} (\bm{k}_0^{\,} + \bm{k}) = \sum _{M = 1}^{N} \frac{1}{M!} \, \frac{\partial ^M \varepsilon _n^{\,} (\bm{k}_0^{\,} + \lambda \, \bm{k})}{\partial \lambda ^M} \bigg |_{\lambda = 0} + \mathcal {O} {\left(k^{N+1} \right)} \end{equation}$$$In order to compute the Taylor expansion of the band dispersion, it is sufficient to use this trick along with the Feynman–Hellmann formula, whose extension to higher orders will be stated and derived below.

The Method

Before proceeding to the method, we state the important assumptions. The Hamiltonian is a finite dimensional matrix that is a smooth function of k. This is always the case with tight-binding models, however large and complicated they may be. We allow for the band in question to be degenerate with some other bands but demand that its Taylor expansion to order N also coincides with the other bands. This is a restrictive assumption which will be relaxed in a future work when the method will be extended to treat the more complicated scenarios of bands touching and intersecting each other in complicated ways (like a linear band crossing for example. Such multi-band intersections, sometimes collectively referred to as multi-fold fermions, can involve several bands and may be stabilized by symmetries.^[32]).

We are interested in computing the Taylor expansion at $$\bm{k}_0^{\,}$$ to order N, of the n^th band (denoted by $$\varepsilon _n^{\,}(\bm{k})$$) arising from the Hamiltonian matrix $$\widehat{H}(\bm{k})$$. The method proceeds to compute the terms order by order.

• 1.

  Compute numerically, all the eigenvalues and eigenvectors of $$\widehat{H}(\bm{k})$$ at the chosen point $$\bm{k}_0^{\,}$$. Let us order these eigenvalues and denote them by $$\varepsilon _1^{\,}$$, $$\varepsilon _2^{\,}$$,…, $$\varepsilon _{\mathcal {N}}^{\,}$$. The corresponding eigenvectors are arranged into the columns of a matrix $$\widehat{E}$$.

• 2.

  Perform a change of basis from the orbital basis to the eigenbasis at $$\bm{k}_0^{\,}$$ by transforming the Hamiltonian: $$\widehat{H}(\bm{k}) \rightarrow \widehat{E}^{\dag } \, \widehat{H}(\bm{k}) \, \widehat{E}$$. (Note that we are transforming the Hamiltonian as a function of k, not just the Hamiltonian at the chosen $$\bm{k}_0^{\,}$$. The eigenvectors that constitute the new basis are however the numerically computed ones at $$\bm{k}_0^{\,}$$. Technically this change of basis is not needed. But it makes the algorithm more efficient since we will not have to repeatedly compute matrix elements of the form $$\langle m_1^{\,} | \partial ^{M} \widehat{H} | m_2^{\,} \rangle$$ with respect to the eigenvectors $$|m_1^{\,} \rangle$$ and $$|m_2^{\,} \rangle$$ by matrix multiplication, and can instead directly read off the corresponding matrix elements).

• 3.

  Calculate the following matrices analytically by directly differentiating the basis transformed Hamiltonian: 4 $$$\begin{eqnarray} &&\partial \widehat{H} \coloneqq \frac{\partial \widehat{H}(\bm{k}_0^{\,} + \lambda \, \bm{k})}{\partial \lambda } \bigg |_{\lambda = 0}, \partial ^2 \widehat{H} \coloneqq \frac{\partial ^2 \widehat{H}(\bm{k}_0^{\,} + \lambda \, \bm{k})}{\partial \lambda ^2} \bigg |_{\lambda = 0}, \ldots ,\nonumber\\ &&\partial ^N \widehat{H} \coloneqq \frac{\partial ^N \widehat{H}(\bm{k}_0^{\,} + \lambda \bm{k})}{\partial \lambda ^N} \bigg |_{\lambda = 0} \end{eqnarray}$$$These matrices when divided by the appropriate factorial weights and assembled together make up the Taylor expansion of the Hamiltonian at $$\bm{k}_0^{\,}$$ to order N.

• 4.

  To compute the Taylor expansion, we shall now outline the procedure to compute the following polynomials and assemble them into the Taylor series as in Equation (3): 5 $$$\begin{eqnarray} &&\partial \varepsilon _n^{\,} \coloneqq \frac{\partial \varepsilon _n^{\,}(\bm{k}_0^{\,} + \lambda \, \bm{k})}{\partial \lambda } \bigg |_{\lambda = 0}, \partial ^2 \varepsilon _n^{\,} \coloneqq \frac{\partial ^2 \varepsilon _n^{\,}(\bm{k}_0^{\,} + \lambda \, \bm{k})}{\partial \lambda ^2} \bigg |_{\lambda = 0}, \ldots ,\nonumber\\ &&\partial ^N \varepsilon _n^{\,} \coloneqq \frac{\partial ^N \varepsilon _n^{\,}(\bm{k}_0^{\,} + \lambda \, \bm{k})}{\partial \lambda ^N} \bigg |_{\lambda = 0} \end{eqnarray}$$$

• 5.

  Before computing each of the above polynomials, we formally define the following recursive function for $$1 \leqslant m, n \leqslant \dim \widehat{H}$$, with $$\varepsilon _m^{\,} \ne \varepsilon _n^{\,}$$ and κ ⩾ 1: 6 $$$\begin{eqnarray} &&f(\kappa , m, n) = \frac{{\left[ \partial ^{\kappa } \widehat{H} \right]}_{m,n}^{\,}}{\varepsilon _n^{\,} - \varepsilon _m^{\,}} + \sum \limits _{\widetilde{\kappa } = 1}^{\kappa - 1} \def\eqcellsep{&}\begin{pmatrix} \kappa \\ \widetilde{\kappa } \end{pmatrix}\nonumber\\ &&\quad \sum \limits _{| \widetilde{m} \rangle \notin D(n)} \frac{ {\left[ \partial ^{\kappa -\widetilde{\kappa }} \widehat{H} \right]}_{m,\widetilde{m}}^{\,} - {\left(\partial ^{\kappa -\widetilde{\kappa }}\varepsilon _n^{\,} \right)} \, \delta _{m,\widetilde{m}}^{\,} }{\varepsilon _n^{\,} - \varepsilon _m^{\,}} \, f {\left(\widetilde{\kappa }, \widetilde{m}, n \right)} \end{eqnarray}$$$where $$[\partial ^{\kappa } \widehat{H}]_{m,n}^{\,}$$ is a shorthand for $$\langle m | \partial ^{\kappa } \widehat{H} |n \rangle$$, a matrix element of the matrix $$\partial ^{\kappa } \widehat{H}$$, $$\def\eqcellsep{&}\begin{pmatrix} \kappa \\ \widetilde{\kappa } \end{pmatrix}$$ is the binomial coefficient and D(n) is the set of all the eigenstates that are degenerate to |n⟩ (including itself), so that the $$|\widetilde{m}\rangle$$ sum runs over all the eigenvectors at $$\bm{k}_0^{\,}$$ that are not degenerate to |n⟩.

  Given that this is a recursive definition, we will obtain a sequence of product of the function f evaluated at different $$\widetilde{m}$$ and decreasing $$\widetilde{\kappa }$$ in the full expansion of the right hand side. We list the explicit expression for f(κ, m, n) for κ = 1, 2 7a $$$\begin{align} f(1,m,n) = \frac{[\partial \widehat{H}]_{m,n}}{\varepsilon _n^{\,} - \varepsilon _m^{\,}} \end{align}$$$7b $$$\begin{eqnarray} f(2,m,n) = \frac{[\partial ^2 \widehat{H}]_{m,n}}{\varepsilon _n^{\,} - \varepsilon _m^{\,}} + 2 \sum _{|\widetilde{m}\rangle \notin D(n)} \frac{[\partial \widehat{H}]_{n,\widetilde{m}}}{\varepsilon _n^{\,} - \varepsilon _m^{\,}} \, \frac{[\partial \widehat{H}]_{\widetilde{m},n}}{\varepsilon _n^{\,} - \varepsilon _{\widetilde{m}}^{\,}} - 2 \, \frac{\partial \varepsilon _n^{\,} \, [\partial \widehat{H}]_{m,n}}{{\left(\varepsilon _n^{\,} - \varepsilon _{m}^{\,} \right)}^2} \nonumber\hspace*{-10pt}\\ \end{eqnarray}$$$

  The sequence of product of f's always terminates when $$\widetilde{\kappa }=1$$, and we have in that case $$f (\widetilde{\kappa } = 1, \widetilde{m}, n ) = [\partial \widehat{H}]_{\widetilde{m},n}^{\,} / (\varepsilon _n^{\,} - \varepsilon _{\widetilde{m}}^{\,})$$. For any κ > 1 and given m and n, to evaluate f(κ, m, n), we need to know all the lower derivatives of the dispersion of band n, denoted by $$\partial \varepsilon _n^{\,}, \partial ^2 \varepsilon _n^{\,}, \ldots , \partial ^{\kappa -1} \varepsilon _n^{\,}$$. But since we know the first derivative from Feynman–Hellmann lemma, $$\partial \varepsilon _n^{\,} = [ \partial \widehat{H} ]_{n,n}^{\,}$$, we can explicitly evaluate all the higher derivatives from $$\partial ^2 \varepsilon _n^{\,}$$ onwards order by order using the extended Feynman–Hellmann formula (see Appendix for proof): 8 $$$\begin{eqnarray} \partial ^M \varepsilon _n^{\,} = {\left[ \partial ^M \widehat{H} \right]}_{n,n}^{\,} + \sum _{\kappa = 1}^{M - 1} \def\eqcellsep{&}\begin{pmatrix} M \\ \kappa \end{pmatrix} \sum _{|m\rangle \notin D(n)} {\left[ \partial ^{M-\kappa } \widehat{H} \right]}_{n,m} \, f(\kappa , m, n)\nonumber\\ \end{eqnarray}$$$

• 6.

  Having computed $$\partial \varepsilon _n^{\,}$$ through $$\partial ^N \varepsilon _n^{\,}$$, we assemble them into the Taylor polynomial of degree N at $$\bm{k}_0^{\,}$$: 9 $$$\begin{equation} j^N (\varepsilon _n^{\,}) |_{\bm{k}_0} = \varepsilon _n^{\,} + \frac{\partial \varepsilon _n^{\,}}{1!} + \frac{\partial ^2 \varepsilon _n^{\,}}{2!} + \cdots + \frac{\partial ^N \varepsilon _n^{\,}}{N!} \end{equation}$$$

Before we proceed further, it is important to emphasize that the method can be also used to track the changes in the dispersion due to small changes in the tuning parameters. This can be achieved by simply redefining $$\partial ^N \widehat{H}$$ in Equation (4). For example, consider the case where we have the Hamiltonian of the system parametrized for changes in some hopping term t around a critical value $$t_0^{\,}$$ (this could be the result of twist, strain, pressure, etc). To obtain the effect of the tuning directly in the Taylor expansion of the dispersion, we redefine Equation (4) as follows and apply the rest of the algorithm 10 $$$\begin{eqnarray} {{ \def\eqcellsep{&}\begin{array}{l} \partial \widehat{H} \coloneqq \displaystyle\frac{\partial \widehat{H}(\bm{k}_0^{\,} + \lambda \, \bm{k}, t_0^{\,} + \lambda \, \delta t)}{\partial \lambda } \bigg |_{\lambda = 0}, \partial ^2 \widehat{H} \coloneqq \displaystyle\frac{\partial ^2 \widehat{H}(\bm{k}_0^{\,} + \lambda \, \bm{k}, t_0^{\,} + \lambda \, \delta t)}{\partial \lambda ^2} \bigg |_{\lambda = 0}, \ldots \\[12pt] \quad \ldots , \partial ^N \widehat{H} \coloneqq \displaystyle\frac{\partial ^N \widehat{H}(\bm{k}_0^{\,} + \lambda \bm{k}, t_0^{\,} + \lambda \, \delta t)}{\partial \lambda ^N} \bigg |_{\lambda = 0} \\[-32pt] \end{array} }}\nonumber\hspace*{-10pt}\\ \end{eqnarray}$$$This will then yield terms in the Taylor expansion that contain δt, δt², … which will track the changes in the dispersion under the tuning of parameters. This procedure is, by no means, restricted to a single tuning parameter. If we search for unusual band geometries under tuning, this procedure can be very helpful. All this is illustrated using an example in Section 3.

Unpacking the Formula

Before proceeding with the application, we need to state the explicit form of the linear, quadratic and cubic derivatives arising from Equations 6 and 8 which will be used later. The detailed derivations are included in the Appendix. 11a $$$\begin{align} \partial \varepsilon _n^{\,} =& \langle n | \partial \widehat{H} | n \rangle \end{align}$$$11b $$$\begin{align} \partial ^2 \varepsilon _n^{\,} =& \langle n | \partial ^2 \widehat{H} | n \rangle + 2 \sum _{|m_1 \rangle \notin D(n)} \frac{\langle n | \partial \widehat{H} | m_1^{\,} \rangle \langle m_1^{\,} | \partial \widehat{H} | n \rangle }{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \end{align}$$$11c $$$\begin{eqnarray} \partial ^3 \varepsilon _n^{\,} &=& \langle n | \partial ^3 \widehat{H} | n \rangle + 3 \sum _{|m_1 \rangle \notin D(n)} \frac{\langle n | \partial ^2 \widehat{H} | m_1^{\,} \rangle \langle m_1^{\,} | \partial \widehat{H} | n \rangle }{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}}\nonumber \\ && +\ 3 \sum _{|m_1 \rangle \notin D(n)} \frac{\langle n | \partial \widehat{H} | m_1^{\,} \rangle \langle m_1^{\,} | \partial ^2 \widehat{H} | n \rangle }{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \nonumber \\ && +\ 6 \sum _{|m_1 \rangle \notin D(n)} \sum _{|m_2 \rangle \notin D(n)} \frac{\langle n | \partial \widehat{H} | m_1^{\,} \rangle \langle m_1^{\,} | \partial \widehat{H} | m_2^{\,} \rangle \langle m_2^{\,} | \partial \widehat{H} | n \rangle }{(\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,})(\varepsilon _n^{\,} - \varepsilon _{m_2}^{\,})} \nonumber \\ && -\ 6 \, \partial \varepsilon _n^{\,} \sum _{|m_1 \rangle \notin D(n)} \frac{\langle n | \partial \widehat{H} | m_1^{\,} \rangle \langle m_1^{\,} | \partial \widehat{H} | n \rangle }{(\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,})^2} \end{eqnarray}$$$Note that the formulae listed above only require us to be able to numerically diagonalize the Hamiltonian $$\widehat{H}(\bm{k})$$ at a chosen $$\bm{k}_0^{\,}$$, and be able to differentiate $$\widehat{H}(\bm{k})$$ analytically a finite number of times, all of which can be readily and quickly done using a computer algebra system (such as Mathematica^[33] or Maple^[34]). Knowing the eigenvalues at $$\bm{k}_0^{\,}$$ and the matrix elements of various derivatives of $$\widehat{H}(\bm{k})$$, we shall be able to compute the Taylor expansion order by order, with each higher order depending upon the lower order terms. The first order term (Equation (11a)) is just the familiar Feynman –Hellmann result.

Haldane Model as an Example

In this section, we benchmark the method against the widely studied Haldane model.^[35] An extension of the Haldane model can be tuned to host the monkey saddle singularity having dispersion $$k_y^3 - 3 \, k_x^2 k_y$$.^[29] Although the method can more generally be applied to larger systems that require numerical diagonalization, we choose the Haldane model for its simplicity and analytic tractability, so that it serves as a pedagogical example of how to apply the method. The underlying lattice for the Haldane model is a honeycomb lattice similar to graphene, with two triangular sublattices A and B. The nearest neighbor unit vectors originating from the A sublattice on to the B sublattice are given by 12 $$$\begin{equation} \bm {a}_1^{\phantom{\dag }} = \def\eqcellsep{&}\begin{pmatrix} 1 \\ \\ 0 \end{pmatrix}, \bm {a}_2^{\phantom{\dag }} = \def\eqcellsep{&}\begin{pmatrix} -\frac{1}{2} \\ \\ \frac{\sqrt {3}}{2} \end{pmatrix} , \bm {a}_3^{\phantom{\dag }} = \def\eqcellsep{&}\begin{pmatrix} -\frac{1}{2} \\ \\ -\frac{\sqrt {3}}{2} \end{pmatrix} \end{equation}$$$

The next nearest neighbor vectors are given by $$\bm {b}_1^{\phantom{\dag }} = \bm {a}_2^{\phantom{\dag }} - \bm {a}_3^{\phantom{\dag }}$$, $$\bm {b}_2^{\phantom{\dag }} = \bm {a}_3^{\phantom{\dag }} - \bm {a}_1^{\phantom{\dag }}$$ and $$\bm{b}_3^{\phantom{\dag }} = \bm{a}_1^{\phantom{\dag }} - \bm{a}_2^{\phantom{\dag }}$$. The full k-space Hamiltonian is then the sum of nearest neighbor, staggered chemical potential and next nearest neighbor terms: $$\widehat{H}(\bm{k}) = \widehat{H}_0^{\phantom{\dag }} (\bm{k}) + M \sigma _z^{\phantom{\dag }} + 2 t_2^{\phantom{\dag }} \sum _i \sin (\bm{k} \cdot \bm{b}_i^{\phantom{\dag }}) \, \sigma _z^{\phantom{\dag }}$$, where $$\sigma _z^{\phantom{\dag }}$$ is the familiar Pauli matrix and 13 $$$\begin{equation} \widehat{H}_0^{\phantom{\dag }} (\bm{k}) = \def\eqcellsep{&}\begin{pmatrix} 0 & t_1^{\phantom{\dag }} \sum _i e^{i \bm{k} \cdot \bm{a}_i} \\ t_1^{\phantom{\dag }} \sum _i e^{-i \bm{k} \cdot \bm{a}_i} & 0 \end{pmatrix} \end{equation}$$$Since the Hamiltonian is 2D, it can be diagonalized exactly to give two bands indexed by n = 1, 2 14 $$$\begin{equation} \varepsilon _{n}^{\phantom{\dag }} (\bm{k}) = (-1)^n {\left[ {\left(M + 2 t_2 \sum _i \sin (\bm{k} \cdot \bm{b}_i^{\phantom{\dag }}) \right)}^2 + t_1^2 \bigg | \sum _i e^{i \bm{k} \cdot \bm{a}_i^{\phantom{\dag }}} \bigg |^2 \right]}^{\frac{1}{2}} \end{equation}$$$Consider the following point of three-fold rotation symmetry 15 $$$\begin{equation} \bm{k}_-^{\phantom{\dag }} = \frac{4\pi }{3\sqrt {3}} \def\eqcellsep{&}\begin{pmatrix} 0 \\ \\ 1 \end{pmatrix} \end{equation}$$$Let us define a critical value for the staggered chemical potential 16 $$$\begin{equation} M_0^{\,} = \frac{t_1^2 - 18 \, t_2^2}{2 \, \sqrt {3} \, t_2^{\phantom{\dag }}} \end{equation}$$$We can analytically derive the Taylor expansion for the upper band (n = 2) dispersion, which around the $$K_-^{\phantom{\dag }}$$ point reads 17 $$$\begin{eqnarray} \varepsilon _2^{\phantom{\dag }} (\bm{k} + \bm{k}_-^{\phantom{\dag }})& \approx & {\left| M - 3\sqrt {3} \, t_2^{\,} \right|} + \frac{9 \sqrt {3} \, t_2^{\,} (M + M_0^{\,})}{4 {\left| M - 3\sqrt {3} \, t_2^{\,} \right|}} {\left(k_x^2 + k_y^2 \right)} \nonumber \\ && +\ \frac{6 \, t_2^{\,} {\left(\sqrt {3} \, M - 9 \, t_2^{\,} \right)} - 9 \, t_1^2}{16 {\left| M - 3\sqrt {3} \, t_2^{\,} \right|}} {\left(k_y^3 - 3 \, k_y^{\,} \, k_x^2 \right)} + \mathcal {O} {\left(k^4 \right)} \nonumber\\ \end{eqnarray}$$$When we tune $$M \rightarrow - M_0^{\,}$$18 $$$\begin{align} \varepsilon _2^{\phantom{\dag }} (\bm{k} + \bm{k}_-^{\phantom{\dag }}) \approx & \frac{t_1^2}{2 \, \sqrt {3} \, |t_2^{\phantom{\dag }}|} - \frac{3 \, \sqrt {3} \, |t_2^{\phantom{\dag }}|}{2} {\left(k_y^3 - 3 \, k_y^{\phantom{\dag }} \, k_x^2 \right)} + \mathcal {O} {\left(k^4 \right)} \end{align}$$$Thus, we have obtained a monkey saddle. Let us examine what happens when M is close to but not equal to $$-M_0^{\,}$$, i.e $$M = -M_0^{\,} + \delta M$$. Assuming $$t_2^{\,} > 0$$, for very small δM, the Taylor expansion gets modified as 19 $$$\begin{eqnarray} {{ \def\eqcellsep{&}\begin{array}{rcl} \varepsilon _2^{\phantom{\dag }} (\bm{k} + \bm{k}_-^{\phantom{\dag }}) &\approx & {\left| - \displaystyle\frac{t_1^2}{2 \, \sqrt {3} \, t_2^{\,}} + \delta M \right|} + \displaystyle\frac{27 \, t_2^2 \,\, \delta M}{2 \, t_1^2} {\left(k_x^2 + k_y^2 \right)} \\[12pt] && -\ {\left(\displaystyle\frac{3 \, \sqrt {3} \, t_2^{\,}}{2} + \frac{27 \, t_2^2 \,\, \delta M}{4 \, t_1^2} \right)} {\left(k_y^3 - 3 \, k_y^{\phantom{\dag }} \, k_x^2 \right)} + \mathcal {O} {\left(k^4, \delta M^2 \right)} \end{array} }}\nonumber\\ \end{eqnarray}$$$

As we see, a small quadratic term modulated by δM appears in the dispersion, splitting the higher order critical point into an ordinary maximum/minimum at the centre, surrounded by three regular saddle points (see Figure 1). Using the numerical values: t = 1 and t = 1/2: 20 $$$\begin{align} \varepsilon _2^{\phantom{\dag }} (\bm{k} + \bm{k}_-^{\phantom{\dag }}) \approx & \, \frac{1}{\sqrt {3}} - \delta M + \frac{27}{8} \, \delta M {\left(k_x^2 + k_y^2 \right)} \nonumber \\ & \quad - {\left(\frac{3\sqrt {3}}{4} + \frac{27}{16} \, \delta M \right)} {\left(k_y^3 - 3 \, k_y^{\phantom{\dag }} \, k_x^2 \right)} + \mathcal {O} {\left(k^4, \delta M^2 \right)} \end{align}$$$We will now attempt to reproduce this Taylor expansion using the method described above. Firstly we calculate the eigenvalues and eigenvectors at $$\bm{k}_-^{\,}$$ for the critical tuning (since δM is treated here as a perturbation, we are computing eigenvalues and eigenvectors of $$H(\bm{k}_-^{\,} + \lambda \, \bm{k} ; -M_0^{\,} + \lambda \, \delta M)$$ at λ = 0, which is the same as computing them for the critical tuning). The Hamiltonian takes the value 21 $$$\begin{equation} \widehat{H}(\bm{k}_-^{\,}) = \def\eqcellsep{&}\begin{pmatrix} -\frac{1}{\sqrt {3}} & 0 \\ 0 & \frac{1}{\sqrt {3}} \end{pmatrix} \end{equation}$$$Thus, it is already diagonal. The eigenvalues and eigenvectors are then given by 22 $$$\begin{align} \varepsilon _{1}^{\,} (\bm{k}_-^{\,}) &= - \frac{1}{\sqrt {3}}, \quad |E_1^{\,} \rangle = \def\eqcellsep{&}\begin{pmatrix} 1 \\ 0 \end{pmatrix} , \nonumber \\ \varepsilon _{2}^{\,} (\bm{k}_-^{\,}) &= \frac{1}{\sqrt {3}}, \quad |E_2^{\,} \rangle = \def\eqcellsep{&}\begin{pmatrix} 0 \\ 1 \end{pmatrix} \end{align}$$$We now compute the matrices $$\partial \widehat{H}$$, $$\partial ^2 \widehat{H\text{,}}$$ and $$\partial ^3 \widehat{H}$$, defined in Equation (4), in the eigenbasis: 23a $$$\begin{align} \partial \widehat{H} =& {\left( \def\eqcellsep{&}\begin{array}{cc}\delta M & \frac{3}{2} i {\left(k_x+i k_y\right)} \\[3pt] -\frac{3}{2} i {\left(k_x-i k_y\right)} & -\delta M \\[3pt] \end{array} \right)} \end{align}$$$23b $$$\begin{align} \partial ^2 \widehat{H} =& {\left( \def\eqcellsep{&}\begin{array}{cc}\frac{9}{4} \sqrt {3} {\left(k_x^2+k_y^2\right)} & -\frac{3}{4} {\left(k_x-i k_y\right)}{}^2 \\[3pt] -\frac{3}{4} {\left(k_x+i k_y\right)}{}^2 & -\frac{9}{4} \sqrt {3} {\left(k_x^2+k_y^2\right)} \\[3pt] \end{array} \right)} \end{align}$$$23c $$$\begin{align} \partial ^3 \widehat{H} =& {\left( \def\eqcellsep{&}\begin{array}{cc}\frac{9}{8} \sqrt {3} k_y {\left(k_y^2-3 k_x^2\right)} & -\frac{9}{8} i {\left(k_x-i k_y\right)} {\left(k_x+i k_y\right)}{}^2 \\[3pt] \frac{9}{8} i {\left(k_x-i k_y\right)}{}^2 {\left(k_x+i k_y\right)} & -\frac{9}{8} \sqrt {3} k_y {\left(k_y^2-3 k_x^2\right)} \\[3pt] \end{array} \right)} \end{align}$$$The derivatives $$\partial ^N \varepsilon _2^{\,}$$ for N = 1, 2, 3 are then given by 24 $$$\begin{eqnarray} \partial \varepsilon _2^{\,} &=& \langle E_2^{\,} | \partial \widehat{H} | E_2^{\,} \rangle \nonumber \\ &=& -\delta M, \nonumber \\ \partial ^2 \varepsilon _2^{\,} &=& \langle E_2^{\,} | \partial ^2 \widehat{H} | E_2^{\,} \rangle + 2 \, \frac{\langle E_2^{\,} | \partial \widehat{H} | E_1^{\,} \rangle \langle E_1^{\,} | \partial \widehat{H} | E_2^{\,} \rangle }{\varepsilon _2^{\,} (\bm{k}_-^{\,}) - \varepsilon _1^{\,} (\bm{k}_-^{\,})} \nonumber \\ &=& -\frac{9\sqrt {3}}{4} {\left(k_x^2 + k_y^2 \right)} + 2 \, \frac{-\frac{3 i}{2} {\left(k_x^{\,} - i k_y^{\,} \right)} \times \frac{3 i}{2} {\left(k_x^{\,} + i k_y^{\,} \right)}}{ \frac{1}{\sqrt {3}} - {\left(- \frac{1}{\sqrt {3}} \right)}} \nonumber \\ &=& 0, \nonumber \\ \partial ^3 \varepsilon _n^{\,} &=& \langle E_2^{\,} | \partial ^3 \widehat{H} | E_2^{\,} \rangle + 3 \, \frac{\langle E_2^{\,} | \partial ^2 \widehat{H} | E_1^{\,} \rangle \langle E_1^{\,} | \partial \widehat{H} | E_2^{\,} \rangle }{\varepsilon _2^{\,} (\bm{k}_-^{\,}) - \varepsilon _1^{\,} (\bm{k}_-^{\,})}\nonumber\\ && +\ 3 \, \frac{\langle E_2^{\,} | \partial \widehat{H} | E_1^{\,} \rangle \langle E_1^{\,} | \partial ^2 \widehat{H} | E_2^{\,} \rangle }{\varepsilon _2^{\,} (\bm{k}_-^{\,}) - \varepsilon _1^{\,} (\bm{k}_-^{\,})} \nonumber \\ && +\ 6 \, \frac{\langle E_2^{\,} | \partial \widehat{H} | E_1^{\,} \rangle \langle E_1^{\,} | \partial \widehat{H} | E_1^{\,} \rangle \langle E_1^{\,} | \partial \widehat{H} | E_2^{\,} \rangle }{{\left(\varepsilon _2^{\,} (\bm{k}_-^{\,}) - \varepsilon _1^{\,} (\bm{k}_-^{\,}) \right)} {\left(\varepsilon _2^{\,} (\bm{k}_-^{\,}) - \varepsilon _1^{\,} (\bm{k}_-^{\,}) \right)}} \nonumber \\ && -\ 6 \, \partial \varepsilon _2^{\,} \, \frac{\langle E_2^{\,} | \partial \widehat{H} | E_1^{\,} \rangle \langle E_1^{\,} | \partial \widehat{H} | E_2^{\,} \rangle }{{\left(\varepsilon _2^{\,} (\bm{k}_-^{\,}) - \varepsilon _1^{\,} (\bm{k}_-^{\,}) \right)}^2} \nonumber \\ &=& - \frac{9 \sqrt {3}}{8} {\left(k_y^3 - 3 \, k_x^2 k_y^{\,} \right)}\nonumber \\ && +\ 3 \, \frac{2 \, \text{Re} {\left[ -\frac{3}{4} {\left(k_x^{\,} + i k_y^{\,} \right)}^2 \times \frac{3 i}{2} {\left(k_x^{\,} + i \, k_y^{\,} \right)} \right]}}{\frac{1}{\sqrt {3}} - {\left(- \frac{1}{\sqrt {3}} \right)}} \nonumber \\ && +\ 6 \, \frac{-\frac{3 i}{2} {\left(k_x^{\,} - i k_y^{\,} \right)} \times \delta M \times \frac{3 i}{2} {\left(k_x^{\,} + i k_y^{\,} \right)}}{ {\left(\frac{1}{\sqrt {3}} - {\left(- \frac{1}{\sqrt {3}} \right)} \right)}^2} \nonumber \\ && -\ 6 (- \delta M) \, \frac{-\frac{3 i}{2} {\left(k_x^{\,} - i k_y^{\,} \right)} \times \frac{3 i}{2} {\left(k_x^{\,} + i k_y^{\,} \right)}}{ {\left(\frac{1}{\sqrt {3}} - {\left(- \frac{1}{\sqrt {3}} \right)} \right)}^2} \nonumber \\ &= & - \frac{9 \sqrt {3}}{2} {\left(k_y^3 - 3 \, k_x^2 k_y^{\,} \right)} + \frac{81}{4} \, \delta M {\left(k_x^2 + k_y^2 \right)} \end{eqnarray}$$$We can assemble these into the Taylor expansion using Equation (9): 25 $$$\begin{eqnarray} {{ \def\eqcellsep{&}\begin{array}{l} \varepsilon _2^{\phantom{\dag }} (\bm{k} + \bm{k}_-^{\phantom{\dag }}) \approx \, \displaystyle\frac{1}{\sqrt {3}} - \delta M + \frac{27}{8} \, \delta M {\left(k_x^2 + k_y^2 \right)} - \displaystyle\frac{3\sqrt {3}}{4} {\left(k_y^3 - 3 \, k_y^{\phantom{\dag }} \, k_x^2 \right)} + \cdots \end{array} }}\nonumber\\ \end{eqnarray}$$$

[IMAGE OMITTED. SEE PDF]

As it is evident, the Taylor expansion matches with the one we obtained by directly Taylor expanding the band dispersion, except for the extra δM modulation of the coefficient of the cubic term. The reason why we did not obtain that term is that it arises from the expansion to fourth order in the auxiliary parameter λ (basically from a term of the form $$\lambda ^4 \, \delta M \, (k_y^3 - 3 k_x^2 k_y )$$), while here we have expanded only up to cubic order for the sake of brevity. Nevertheless, the main feature, which is the appearance of the small quadratic part that perturbs the monkey saddle when we tune off the critical value of the staggered chemical potential, is quite manifest here.

Conclusion

In this paper, we set in detail the method to find HOVHS in practice, when the Hamiltonian H(k) is known by various means. The knowledge of H(k) can come from tight-binding models or ab-initio calculations or even motivated by experimental results. The fundamental idea outlined above was to use an extension the Feynman–Hellmann theorem to compute the various derivatives of the dispersion of any specified band using the technique explained in Section 2.2. These derivatives were then used to assemble the Taylor expansion of the band to any specified order, at arbitrary k points, with the added possibility of incorporating the influence of tuning of any parameters in the Hamiltonian. The latter result is particularly timely, given the increasing interest and the remarkable recent progress in experimental methods to engineer exotic phases in materials by tuning them (via stress, twist, pressure, bias voltage, etc). Widely studied metallic systems like Sr₂RuO₄ or semiconducting ones like two dimensional InSe and In₂Se₃ as well as a wealth of other quantum materials, can be readily explored for the possibility to host HOVHS by this method. This is exactly where the fundamental computational method developed here, will demonstrate its full potential.

In the case of a non-trivial degeneracy at a point in k-space (such as band crossing, band touching etc.), the process of the computation of the derivatives and terms in the series expansion becomes more complicated. This is the subject of a future work and beyond the scope of the present.

Note

A readily usable implementation of the algorithm in Wolfram Mathematica is available.^[36]

Appendix A: Derivation of the Expressions

In this section, we explicitly prove the formula in Equation (8) for linear, quadratic and cubic derivatives of the band. For the sake of convenience of notation, we drop the “hats” on the Hamiltonian $$\widehat{H}$$ and write it simply as H. We further assume that all the eigenvectors are normalised, that is, ⟨n|n⟩ = 1.

A.1 First Derivative

The eigenvalue equation reads A1 $$$\begin{equation} H | n \rangle = \varepsilon _n^{\,} | n \rangle \end{equation}$$$Differentiating this once and applying the product rule for derivatives, we get A2 $$$\begin{equation} \partial H | n \rangle + H | \partial n \rangle = \partial \varepsilon _n^{\,} | n \rangle + \varepsilon _n^{\,} | \partial n \rangle \end{equation}$$$(There is no problem with applying product rule since we are working with a finite dimensional vector space where all the terms in the eigenvalue equation are simply finite sums of product of vector and matrix coefficients in some basis). Now, we take inner product of this equation with |n⟩ A3 $$$\begin{align} \langle n | \partial H | n \rangle + \langle n | H | \partial n \rangle &= \partial \varepsilon _n^{\,} \langle n | n \rangle + \varepsilon _n^{\,} \langle n | \partial n \rangle , \nonumber \\ \Rightarrow \partial \varepsilon _n^{\,} &= \langle n | \partial H | n \rangle \end{align}$$$This is of course the familiar Feynman–Hellmann result. We now take inner product of Equation (A2) with another eigenstate |m⟩ that is not degenerate to |n⟩: A4 $$$\begin{align} \langle m | \partial H | n \rangle + \langle m | H | \partial n \rangle &= \partial \varepsilon _n^{\,} \langle m | n \rangle + \varepsilon _n^{\,} \langle m | \partial n \rangle , \nonumber \\ \Rightarrow \langle m | \partial n \rangle &= \frac{\langle m | \partial H | n \rangle }{\varepsilon _n^{\,} - \varepsilon _m^{\,}} \end{align}$$$Thus, we can find the overlap of the derivative of the eigenvector |n⟩, denoted by |∂n⟩ with the other eigenvectors not degenerate to it. It is not possible within this framework to compute ⟨n|∂n⟩ or even overlap with eigenvectors degenerate to it. We can, however, differentiate ⟨n|n⟩ = 1 to obtain Re ⟨n|∂n⟩ = 0. The imaginary part is of course the famous Berry connection that can not be gauged out in general (specifically when we have topological bands). The fact that we can not precisely compute ⟨n|∂n⟩ will not be an impediment to computing the higher derivatives of the band dispersion $$\varepsilon _n^{\,} (\bm{k})$$, as we shall show below. We also compute the overlap of Equation (A2) with respect to another eigenstate that is degenerate to |n⟩, denoted by $$| n_1^{\,} \rangle \in D(n)$$ (where D(n) is the set of eigenvectors degenerate to n) A5 $$$\begin{align} \langle n_1^{\,} | \partial H | n \rangle + \langle n_1^{\,} | H | \partial n \rangle &= \partial \varepsilon _n^{\,} \langle n_1^{\,} | n \rangle + \varepsilon _n^{\,} \langle n_1^{\,} | \partial n \rangle , \nonumber \\ \Rightarrow \langle n_1^{\,} | \partial H | n \rangle &= 0 , \,\, \text{for} \,\, \varepsilon _{n_1}^{\,} = \varepsilon _n^{\,} \,\, \text{but} \, |n_1^{\,} \rangle \ne | n \rangle \end{align}$$$In practice, however, we may find the projection of the first derivative of Hamiltonian matrix ∂H to the degenerate subspace to be non-diagonal. That is, we may find $$\langle n_1^{\,} | \partial H | n_2^{\,} \rangle \ne 0$$, for distinct but degenerate $$|n_1^{\,} \rangle$$ and $$|n_2^{\,} \rangle$$. In such situations, we diagonalize the projection of ∂H in this subspace, and replace the original set of $$|n_i^{\,} \rangle$$ with the new eigenbasis, so that Equation (A5) is automatically satisfied.

A.2 Second Derivative

Differentiating Equation (A2) once again, we get A6 $$$\begin{eqnarray} \partial ^2 H | n \rangle + 2 \, \partial H | \partial n \rangle + H | \partial ^2 n \rangle = \partial ^2 \varepsilon _n^{\,} | n \rangle + 2 \, \partial \varepsilon _n^{\,} | \partial n \rangle + \varepsilon _n^{\,} | \partial ^2 n \rangle \nonumber\hspace*{-10pt}\\ \end{eqnarray}$$$Taking inner product with |n⟩, we get A7 $$$\begin{eqnarray} \langle n | \partial ^2 H | n \rangle &+& 2 \, \langle n | \partial H | \partial n \rangle + \langle n | H | \partial ^2 n \rangle = \partial ^2 \varepsilon _n^{\,} \langle n | n \rangle\nonumber\\ &+& 2 \, \partial \varepsilon _n^{\,} \, \langle n | \partial n \rangle + \varepsilon _n^{\,} \langle n | \partial ^2 n \rangle \end{eqnarray}$$$On rearranging terms, inserting a resolution of identity and splitting into projectors over the degenerate and non-degenerate subspaces, this becomes A8 $$$\begin{eqnarray} \partial ^2 \varepsilon _n^{\,} &=& \langle n | \partial ^2 H | n \rangle + 2 \, \langle n | \partial H | \partial n \rangle - 2 \, \partial \varepsilon _n^{\,} \, \langle n | \partial n \rangle \nonumber \\ &=& \langle n | \partial ^2 H | n \rangle + 2 \sum _{| m \rangle \notin D(n)} \langle n | \partial H | m \rangle \langle m | \partial n \rangle\nonumber\\ && +\ 2 \sum _{|n_1^{\,} \rangle \in D(n)} \langle n | \partial H | n_1^{\,} \rangle \langle n_1^{\,} | \partial n \rangle - 2 \, \partial \varepsilon _n^{\,} \, \langle n | \partial n \rangle \nonumber \\ &=& \langle n | \partial ^2 H | n \rangle + 2 \sum _{| m \rangle \notin D(n)} \frac{\langle n | \partial H | m \rangle \langle m | \partial H | n \rangle }{\varepsilon _n^{\,} - \varepsilon _m^{\,}} \end{eqnarray}$$$where we have used $$\langle n | \partial H | n \rangle = \partial \varepsilon _n^{\,}$$ and $$\langle n | \partial H | n_1^{\,} \rangle = 0$$ when $$| n_1^{\,} \rangle \ne | n \rangle$$ in going to the last line. We now compute the overlap of the terms in Equation (A6) with |m⟩∉D(n): A9 $$$\begin{eqnarray} \langle m_1^{\,} | \partial ^2 H | n \rangle &+& 2 \, \langle m_1^{\,} | \partial H | \partial n \rangle + \langle m_1^{\,} | H | \partial ^2 n \rangle = \partial ^2 \varepsilon _n^{\,} \langle m_1^{\,} | n \rangle \nonumber\\ &+& 2 \, \partial \varepsilon _n^{\,} \langle m_1^{\,} | \partial n \rangle + \varepsilon _n^{\,} \langle m_1^{\,} | \partial ^2 n \rangle \end{eqnarray}$$$Simplifying this further by rearranging terms and introducing another resolution of identity: A10 $$$\begin{eqnarray} {{ \def\eqcellsep{&}\begin{array}{rcl} \langle m_1^{\,} | \partial ^2 n \rangle &= & \displaystyle\frac{1}{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} {\left[ \langle m_1^{\,} | \partial ^2 H | n \rangle + 2 \, \langle m_1^{\,} | \partial H | \partial n \rangle - 2 \, \partial \varepsilon _n^{\,} \langle m_1^{\,} | \partial n \rangle \right]} \\[12pt] &=& \displaystyle\frac{1}{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} {\left[ \langle m_1^{\,} | \partial ^2 H | n \rangle + 2 \sum _{|m_2\rangle \notin D(n)} \langle m_1^{\,} | \partial H | m_2^{\,} \rangle \langle m_2^{\,} | \partial n \rangle \right.} \\[12pt] && {\left.\quad + 2 \displaystyle\sum _{| n_1^{\,} \rangle \in D(n)} \langle m_1^{\,} | \partial H | n_1^{\,} \rangle \langle n_1^{\,} | \partial n \rangle - 2 \, \partial \varepsilon _n^{\,} \langle m_1^{\,} | \partial n \rangle \right]} \\[12pt] &= & \displaystyle\frac{\langle m_1^{\,} | \partial ^2 H | n \rangle }{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} + 2 \sum _{|m_2\rangle \notin D(n)} \displaystyle\frac{\langle m_1^{\,} | \partial H | m_2^{\,} \rangle \langle m_2^{\,} | \partial H | n \rangle }{(\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,})(\varepsilon _n^{\,} - \varepsilon _{m_2}^{\,})}\\ [12pt] && - 2 \, \partial \varepsilon _n^{\,} \displaystyle\frac{\langle m_1^{\,} | \partial H | n \rangle }{(\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,})^2} \\[12pt] &&+ 2 \sum _{n_1^{\,} \in D(n)} \displaystyle\frac{\langle m_1^{\,} | \partial H | n_1^{\,} \rangle }{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \langle n_1^{\,} | \partial n \rangle \end{array} }}\nonumber\hspace*{-10pt}\\ \end{eqnarray}$$$This will be used in the computation of the third derivative below. Note that once again we are not able to explicitly evaluate the overlap with |n⟩. Lastly, we compute the overlap of Equation (A8) with $$| n_1^{\,} \rangle \in D(n)$$ and $$| n_1^{\,} \rangle \ne |n \rangle$$A11 $$$\begin{eqnarray} \langle n_1^{\,} | \partial ^2 H | n \rangle &+& 2 \, \langle n_1^{\,} | \partial H | \partial n \rangle + \langle n_1^{\,} | H | \partial ^2 n \rangle = \partial ^2 \varepsilon _n^{\,} \langle n_1^{\,} | n \rangle\nonumber\\ & +& 2 \, \partial \varepsilon _n^{\,} \langle n_1^{\,} | \partial n \rangle + \varepsilon _n^{\,} \langle n_1^{\,} | \partial ^2 n \rangle \end{eqnarray}$$$Using $$\langle n_1^{\,} | H = \varepsilon _n^{\,} \langle n_1^{\,} |$$, introducing resolution of identity and simplifying A12 $$$\begin{align} &\langle n_1^{\,} | \partial ^2 H | n \rangle + 2 \sum _{|m_1 \rangle \notin D(n)} \langle n_1^{\,} | \partial H | m_1^{\,} \rangle \langle m_1^{\,} | \partial n \rangle\nonumber\\ &\quad + 2 \sum _{|n_2 \rangle \notin D(n)} \langle n_1^{\,} | \partial H | n_2^{\,} \rangle \langle n_2^{\,} | \partial n \rangle = 2 \, \partial \varepsilon _n^{\,} \langle n_1^{\,} | \partial n \rangle , \nonumber \\ &\Rightarrow \langle n_1^{\,} | \partial ^2 H | n \rangle + 2 \sum _{|m_1 \rangle \notin D(n)} \frac{\langle n_1^{\,} | \partial H | m_1^{\,} \rangle \langle m_1^{\,} | \partial H | n \rangle }{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} = 0 \end{align}$$$The expression on the left hand side resembles the expression for $$\partial ^2 \varepsilon _n^{\,}$$, with the difference that the leftmost bra in each of the terms is $$|n_1^{\,}\rangle$$ instead of |n⟩, and the expression ultimately evaluates to 0 rather than $$\partial ^2 \varepsilon _n^{\,}$$.

A.3 Third Derivative

The third derivative of the eigenvalue equation reads A13 $$$\begin{eqnarray} \partial ^3 H | n \rangle &+& 3 \, \partial ^2 H | \partial n \rangle + 3 \, \partial H | \partial ^2 n \rangle + H | \partial ^3 n \rangle = \partial ^3 \varepsilon _n^{\,} | n \rangle\nonumber\\ & +& 3 \, \partial ^2 \varepsilon _n^{\,} | \partial n \rangle + 3 \, \partial \varepsilon _n^{\,} | \partial ^2 n \rangle + \varepsilon _n^{\,} | \partial ^3 n \rangle \end{eqnarray}$$$We compute the overlap of this equation with |n⟩, rearrange terms, introduce resolutions of identity (splitting them into projectors over the degenerate and non-degenerate subspaces) and do some simplifications A14 $$$\begin{eqnarray} \partial ^3 \varepsilon _n^{\,} &=& \langle n | \partial ^3 H | n \rangle + 3 \langle n | \partial ^2 H | \partial n \rangle + 3 \langle n | \partial H | \partial ^2 n \rangle \nonumber\\ && -\ 3 \, \partial ^2 \varepsilon _n^{\,} \, \langle n | \partial n \rangle - 3 \, \partial \varepsilon _n^{\,} \, \langle n | \partial ^2 n \rangle \nonumber \\ &=& \langle n | \partial ^3 H | n \rangle + 3 \sum _{|m_1 \rangle \notin D(n)} \langle n | \partial ^2 H | m_1^{\,} \rangle \langle m_1^{\,} | \partial n \rangle\nonumber\\ && +\ 3 \sum _{|m_1 \rangle \notin D(n)} \langle n | \partial H | m_1^{\,} \rangle \langle m_1^{\,} | \partial ^2 n \rangle \nonumber \\ && +\ 3 \sum _{| n_1^{\,} \rangle \in D(n)} \langle n | \partial ^2 H | n_1^{\,} \rangle \langle n_1^{\,} | \partial n \rangle\nonumber \\ && +\ 3 \sum _{| n_1^{\,} \rangle \in D(n)} \langle n | \partial H | n_1^{\,} \rangle \langle n_1^{\,} | \partial ^2 n \rangle \nonumber \\ && -\ 3 \,\, \partial ^2 \varepsilon _n^{\,} \, \langle n | \partial n \rangle - 3 \,\, \partial \varepsilon _n^{\,} \, \langle n | \partial ^2 n \rangle \end{eqnarray}$$$We can substitute for $$\langle m_1^{\,} | \partial n \rangle$$ and $$\langle m_1^{\,} | \partial ^2 n \rangle$$ from Equations A4 and A10, respectively, to obtain A15 $$$\begin{eqnarray} \partial ^3 \varepsilon _n^{\,} &=& \langle n | \partial ^3 H | n \rangle + 3 \sum _{|m_1 \rangle \notin D(n)} \frac{\langle n | \partial ^2 H | m_1^{\,} \rangle \langle m_1^{\,} | \partial H | n \rangle }{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \nonumber \\ && +\ 3 \sum _{|m_1 \rangle \notin D(n)} \frac{\langle n | \partial H | m_1^{\,} \rangle \langle m_1^{\,} | \partial ^2 H | n \rangle }{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \nonumber \\ && +\ 6 \sum _{|m_1 \rangle \notin D(n)} \sum _{|m_2 \rangle \notin D(n)} \frac{\langle n | \partial H | m_1^{\,} \rangle \langle m_1^{\,} | \partial H | m_2^{\,} \rangle \langle m_2^{\,} | \partial H | n \rangle }{(\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,})(\varepsilon _n^{\,} - \varepsilon _{m_2}^{\,})} \nonumber \\ && -\ 6 \,\, \partial \varepsilon _n^{\,} \sum _{|m_1 \rangle \notin D(n)} \frac{\langle n | \partial H | m_1^{\,} \rangle \langle m_1^{\,} | \partial H | n \rangle }{(\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,})^2} \nonumber \\ && +\ (*1) \,\, 6 \sum _{| n_1^{\,} \rangle \in D(n)} \sum _{| m_1 \rangle \notin D(n)} \frac{\langle n | \partial H | m_1^{\,} \rangle \langle m_1^{\,} | \partial H | n_1^{\,} \rangle }{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \langle n_1^{\,} | \partial n \rangle \nonumber \\ && +\ (*2) \,\, 3 \sum _{| n_1^{\,} \rangle \in D(n)} \langle n | \partial ^2 H | n_1^{\,} \rangle \langle n_1^{\,} | \partial n \rangle \nonumber \\ && +\ (*3) \,\, 3 \sum _{| n_1^{\,} \rangle \in D(n)} \langle n | \partial H | n_1^{\,} \rangle \langle n_1^{\,} | \partial ^2 n \rangle \nonumber \\ && -\ 3 \,\, \partial ^2 \varepsilon _n^{\,} \, \langle n | \partial n \rangle - 3 \,\, \partial \varepsilon _n^{\,} \, \langle n | \partial ^2 n \rangle \end{eqnarray}$$$Now, let us look at some of the terms marked with a *. To simplify them, we pull out the term containing |n⟩ from the sum over $$|n_1^{\,} \rangle \in D(n)$$ and substitute the expression for $$\partial ^2 \varepsilon _n^{\,}$$: A16 $$$\begin{eqnarray} (*2) &+& (*1) \nonumber \\ &=& 3 \sum _{| n_1^{\,} \rangle \in D(n)} \langle n | \partial ^2 H | n_1^{\,} \rangle \langle n_1^{\,} | \partial n \rangle\nonumber \\ &&+\ 6 \sum _{| n_1^{\,} \rangle \in D(n)} \sum _{| m_1 \rangle \notin D(n)} \frac{\langle n | \partial H | m_1^{\,} \rangle \langle m_1^{\,} | \partial H | n_1^{\,} \rangle }{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \langle n_1^{\,} | \partial n_1^{\,} \rangle \nonumber \\ &=& 3 {\left(\langle n | \partial ^2 H | n \rangle + 2 \sum _{| m_1 \rangle \notin D(n)} \frac{\langle n | \partial H | m_1^{\,} \rangle \langle m_1^{\,} | \partial H | n \rangle }{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \right)} \langle n | \partial n \rangle \nonumber \\ && +\ 3 \sum _{\substack{|n_1^{\,} \rangle \in D(n) \\ | n_1^{\,} \rangle \ne | n \rangle }} \left[ \langle n | \partial ^2 H | n_1^{\,} \rangle\right. \nonumber \\ &&\left.+ 2 \sum _{| m_1 \rangle \notin D(n)} \frac{\langle n | \partial H | m_1^{\,} \rangle \langle m_1^{\,} | \partial H | n_1^{\,} \rangle }{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \right] \langle n_1^{\,} | \partial n_1^{\,} \rangle \nonumber \\ &=& \, 3 \,\, \partial ^2 \varepsilon _n^{\,} \, \langle n | \partial n \rangle \end{eqnarray}$$$where we substitute from Equation (A8) for the terms within parenthesis, while the terms within [ ] vanish due to Equation (A12) (since $$|n_1^{\,} \rangle \ne |n\rangle$$). Lastly, since $$\langle n | \partial H | n \rangle = \partial \varepsilon _n^{\,}$$ (from Equation (A3)) and $$\langle n | \partial H | n_1^{\,} \rangle = 0$$ for $$| n \rangle \ne | n_1^{\,} \rangle$$ (from Equation (A5)), we have A17 $$$\begin{equation} (*3) = 3 \sum _{| n_1^{\,} \rangle \in D(n)} \langle n | \partial H | n_1^{\,} \rangle \langle n_1^{\,} | \partial ^2 n \rangle = 3 \,\, \partial \varepsilon _n^{\,} \, \langle n | \partial ^2 n \rangle \end{equation}$$$Therefore the sum (*1) + (*2) + (*3) cancels exactly with the last line of Equation (A15) so that the expression for the third derivative of the band becomes A18 $$$\begin{eqnarray} \partial ^3 \varepsilon _n^{\,} &=& \langle n | \partial ^3 H | n \rangle + 3 \sum _{|m_1 \rangle \notin D(n)} \frac{\langle n | \partial ^2 H | m_1^{\,} \rangle \langle m_1^{\,} | \partial H | n \rangle }{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}}\nonumber \\ && +\ 3 \sum _{|m_1 \rangle \notin D(n)} \frac{\langle n | \partial H | m_1^{\,} \rangle \langle m_1^{\,} | \partial ^2 H | n \rangle }{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \nonumber \\ && +\ 6 \sum _{|m_1 \rangle \notin D(n)} \sum _{|m_2 \rangle \notin D(n)} \frac{\langle n | \partial H | m_1^{\,} \rangle \langle m_1^{\,} | \partial H | m_2^{\,} \rangle \langle m_2^{\,} | \partial H | n \rangle }{(\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,})(\varepsilon _n^{\,} - \varepsilon _{m_2}^{\,})} \nonumber \\ && -\ 6 \,\, \partial \varepsilon _n^{\,} \sum _{|m_1 \rangle \notin D(n)} \frac{\langle n | \partial H | m_1^{\,} \rangle \langle m_1^{\,} | \partial H | n \rangle }{(\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,})^2} \end{eqnarray}$$$To obtain this, we got rid of the terms containing ⟨n|∂n⟩ and ⟨n|∂²n⟩ by a careful application of the results obtained for the lower derivatives $$\partial \varepsilon _n^{\,}$$ and $$\partial ^2 \varepsilon _n^{\,}$$, along with the expressions involving the overlaps $$\langle n_1^{\,} | \partial H | n \rangle$$ and $$\langle n_1^{\,} | \partial ^2 H | n \rangle$$ for $$|n_1^{\,} \rangle \in D(n)$$ and $$|n_1^{\,} \rangle \ne | n \rangle$$.

Appendix B: General Proof

In this section, we prove the general case by induction. Since the proof is somewhat involved, we provide a brief outline of it. While computing the expression for $$\partial ^N \varepsilon _n^{\,}$$, we will encounter terms like $$\langle m_1^{\,} | \partial ^M m_2^{\,} \rangle$$, which can not be evaluated in general. The strategy will be to expand such terms (particularly when $$|m_1^{\,} \rangle$$ and $$|m_2^{\,} \rangle$$ are non-degenerate), substitute the expansions in the expression for $$\partial ^N \varepsilon _n^{\,}$$, and get the pathological terms to cancel out. We first apply the Leibniz theorem for higher derivatives of products to differentiate the eigenvalue equation: B1a $$$\begin{align} \partial ^{N} (H \, |n \rangle ) &= \partial ^{N} (\varepsilon _n^{\,} \, |n \rangle ) \end{align}$$$B1b $$$\begin{align} \Rightarrow \sum _{K=0}^{N} \def\eqcellsep{&}\begin{pmatrix} N \\ K \end{pmatrix} \partial ^{N-K} H \, | \partial ^{K} n \rangle &= \sum _{K=0}^{N} \def\eqcellsep{&}\begin{pmatrix} N \\ K \end{pmatrix} \partial ^{N-K} \varepsilon _n^{\,} \, | \partial ^{K} n \rangle \end{align}$$$ Isolating the K = 0 and K = N terms on both sides of the equation, we obtain B2 $$$\begin{eqnarray} H \, |\partial ^{N} n \rangle &+& \partial ^N H \, |n \rangle + \sum _{K=1}^{N-1} \def\eqcellsep{&}\begin{pmatrix} N \\ K \end{pmatrix} \partial ^{N-K} H \, | \partial ^{K} n \rangle = \varepsilon _n^{\,} \, |\partial ^{N} n \rangle\nonumber\\ &+& \partial ^N \varepsilon _n^{\,} \, |n \rangle + \sum _{K=1}^{N-1} \def\eqcellsep{&}\begin{pmatrix} N \\ K \end{pmatrix} \partial ^{N-K} \varepsilon _n^{\,} \, | \partial ^{K} n \rangle \end{eqnarray}$$$After taking inner product on both sides with |m⟩∉D(n) (i.e., an eigenstate not degenerate to |n⟩), using $$\langle m | H = \varepsilon _m^{\,} \, \langle m|$$ and rearranging, we obtain: B3 $$$\begin{eqnarray} &&\langle m | \partial ^{N} n \rangle = \frac{1}{\varepsilon _n^{\,} - \varepsilon _m^{\,}} \Bigg[ \langle m | \partial ^N H | n \rangle + \sum _{K = 1}^{N-1} \def\eqcellsep{&}\begin{pmatrix} N \\ K \end{pmatrix}\nonumber\\ &&\quad \left(\langle m | \partial ^{N-K} H ** | \partial ^{K} n \rangle- \partial ^{N-K} \varepsilon _n^{\,} \, \langle m | \partial ^{K} n \rangle \right) \Bigg] \end{eqnarray}$$$

We now introduce a resolution of identity in terms of the eigenvectors of the Hamiltonian at the given $$\bm{k}_0^{\,}$$ at the location indicated by ** and separate out the projector to the degenerate subspace $$\sum _{|\tilde{n}\rangle \in D(n)} | \tilde{n} \rangle \langle \tilde{n} |$$. For convenience, we add subscripts to the variables and change N → K_i so that the expression becomes more readily usable in the forthcoming calculations: B4 $$$\begin{multline} \langle m_i^{\,} | \partial ^{K_i} n \rangle = \frac{1}{\varepsilon _n^{\,} - \varepsilon _{m_i}^{\,}} {\left[ \langle m_i^{\,} | \partial ^{K_i} H | n \rangle + \sum _{K_{i+1} = 1}^{K_i-1} \def\eqcellsep{&}\begin{pmatrix} K_i \\ K_{i+1} \end{pmatrix} \sum _{|m_{i+1} \rangle \notin D(n)} \right.}\\ \left\lbrace {\left(\langle m_i^{\,} | \partial ^{K_i - K_{i+1}} H | m_{i+1}^{\,} \rangle - \partial ^{K_i - K_{i+1}} \varepsilon _n^{\,} \, \delta _{m_i, m_{i+1}} \right)} \langle m_{i+1} | \partial ^{K_{i+1}} n \rangle \right \rbrace \\ {\left. + \sum _{K_{i+1} = 1}^{K_i - 1} \def\eqcellsep{&}\begin{pmatrix}K_i \\ K_{i+1} \end{pmatrix} \sum _{|\tilde{n} \rangle \in D(n)} \langle m_i^{\,} | \partial ^{K_i - K_{i+1}} H | \tilde{n} \rangle \langle \tilde{n} | \partial ^{K_{i+1}} n \rangle \right]}\end{multline}$$$

To make the notation less cumbersome, for a given |n⟩, we define B5a $$$\begin{align} & \widetilde{\mathcal {H}}_{m_i, m_{j}}^{K} \coloneqq \langle m_i^{\,} | \partial ^{K} H | m_{j}^{\,} \rangle - \partial ^{K} \varepsilon _n^{\,} \, \delta _{m_i, m_{j}} \end{align}$$$B5b $$$\begin{align} & \mathcal {H}_{m_i,m_j}^{K} \coloneqq \langle m_i^{\,} | \partial ^{K} H | m_j^{\,} \rangle \end{align}$$$B5c $$$\begin{equation} \def\eqcellsep{&}\begin{array}{cc}& {\mathcal{O}}_{{m}_{i},{m}_{j}}^{K}:= \langle {m}_{i}|{\partial}^{K}|{m}_{j} \rangle ,\text{for non-degenerate}|{m}_{i} \rangle ,|{m}_{j} \rangle\end{array} \end{equation}$$$B5d $$$\begin{equation} \def\eqcellsep{&}\begin{array}{cc} & {\mathcal{R}}_{\tilde{n},n}^{K}:= \langle \tilde{n}|{\partial}^{K}n \rangle ,\text{for degenerate}|\tilde{n} \rangle ,|n \rangle . \end{array} \end{equation}$$$ where $$\mathcal {O}$$ is the overlap which will be represented as a recursive function, while the matrix elements $$\widetilde{\mathcal {H}}$$ and $$\mathcal {H}$$ will occur respectively in the middle and at the ends in each of the terms in the final expression for the band derivatives. The reason for separately defining the overlap $$\mathcal {R}$$, which occurs at the end of the sequence of product of terms, is that we can not define a recursive relationship for it unlike $$\mathcal {O}$$, and all of the terms containing $$\mathcal {R}^{\prime }$$s will be shown to cancel out in the final expressions. We can now write the expression for the overlap in Equation (B4) in a compact form: B6 $$$\begin{multline} \mathcal {O}_{m_i,n}^{K_i} = \frac{1}{\varepsilon _n^{\,} - \varepsilon _{m_i}^{\,}} {\left[ \mathcal {H}_{m_i,n}^{K_i} + \sum _{K_{i+1} = 1}^{K_i-1} \def\eqcellsep{&}\begin{pmatrix} K_i \\ K_{i+1} \end{pmatrix} \sum _{|m_{i+1} \rangle \notin D(n)} \widetilde{\mathcal {H}}_{m_i, m_{i+1}}^{K_i - K_{i+1}} \mathcal {O}_{m_{i+1},n}^{K_{i+1}} \right.}\\ {\left. + \sum _{K_{i+1} = 1}^{K_i-1} \def\eqcellsep{&}\begin{pmatrix} K_i \\ K_{i+1} \end{pmatrix} \sum _{|\tilde{n} \rangle \in D(n)} \mathcal {H}_{m_i,\tilde{n}}^{K_i - K_{i+1}} \, \mathcal {R}_{\tilde{n},n}^{K_{i+1}} \right]}\end{multline}$$$This clearly defines a recursive relationship between the overlaps $$\mathcal {O}$$. We now expand this by applying the analogous expression for $$\mathcal {O}_{m_{i+1}, K_{i+1}}$$: B7 $$$\begin{eqnarray} \mathcal {O}_{m_i,n}^{K_i} &= & \frac{\mathcal {H}_{m_i,n}^{K_i}}{\varepsilon _n^{\,} - \varepsilon _{m_i}^{\,}} + \sum _{K_{i+1} = 1}^{K_i-1} \def\eqcellsep{&}\begin{pmatrix} K_i \\ K_{i+1} \end{pmatrix} \sum _{|m_{i+1} \rangle \notin D(n)} \frac{\widetilde{\mathcal {H}}_{m_i, m_{i+1}}^{K_i - K_{i+1}}}{\varepsilon _n^{\,} - \varepsilon _{m_i}^{\,}} {\left[ \frac{\mathcal {H}_{m_{i+1},n}^{K_{i+1}}}{\varepsilon _n^{\,} - \varepsilon _{m_{i+1}}^{\,}} \right.}\nonumber \\ && +\ \sum _{K_{i+2} = 1}^{K_{i+1}-1} \def\eqcellsep{&}\begin{pmatrix} K_{i+1} \\ K_{i+2} \end{pmatrix} \sum _{|m_{i+2} \rangle \notin D(n)} \frac{\widetilde{\mathcal {H}}_{m_{i+1}, m_{i+2}}^{K_{i+1} - K_{i+2}}}{\varepsilon _n^{\,} - \varepsilon _{m_{i+1}}^{\,}} \, \mathcal {O}_{m_{i+2},n}^{K_{i+2}} \nonumber \\ &&{\left. + \sum _{K_{i+2} = 1}^{K_{i+1}-1} \def\eqcellsep{&}\begin{pmatrix} K_{i+1} \\ K_{i+2} \end{pmatrix} \sum _{|\tilde{n} \rangle \in D(n)} \frac{\mathcal {H}_{m_{i+1},\tilde{n}}^{K_{i+1} - K_{i+2}}}{\varepsilon _n^{\,} - \varepsilon _{m_{i+1}}^{\,}} \, \mathcal {R}_{\tilde{n},n}^{K_{i+2}} \right]} \nonumber \\ && +\ (*) \sum _{K_{i+1} = 1}^{K_i-1} \def\eqcellsep{&}\begin{pmatrix} K_i \\ K_{i+1} \end{pmatrix} \sum _{|\tilde{n} \rangle \in D(n)} \frac{\mathcal {H}_{m_i, \tilde{n}}^{K_i - K_{i+1}}}{\varepsilon _n^{\,} - \varepsilon _{m_i}^{\,}} \, \mathcal {R}_{\tilde{n},n}^{K_{i+1}} \end{eqnarray}$$$

Noticeably, the $$\widetilde{\mathcal {H}}$$'s always occur at the beginning and in the middle of the products while both $$\mathcal {H}$$ and $$\mathcal {R}$$ terminate the sequence of products (since the $$\mathcal {R}$$'s can not be expanded further). We rearrange the above equation to collect the terms with $$\mathcal {R}$$'s together (after the (*)). B8 $$$\begin{eqnarray} \mathcal {O}_{m_i,n}^{K_i} &= & \frac{\mathcal {H}_{m_i,n}^{K_i}}{\varepsilon _n^{\,} - \varepsilon _{m_i}^{\,}} + \sum _{K_{i+1} = 1}^{K_i-1} \def\eqcellsep{&}\begin{pmatrix} K_i \\ K_{i+1} \end{pmatrix} \sum _{|m_{i+1} \rangle \notin D(n)} \frac{\widetilde{\mathcal {H}}_{m_i, m_{i+1}}^{K_i - K_{i+1}}}{\varepsilon _n^{\,} - \varepsilon _{m_i}^{\,}} \, \frac{\mathcal {H}_{m_{i+1},n}^{K_{i+1}}}{\varepsilon _n^{\,} - \varepsilon _{m_{i+1}}^{\,}} \nonumber \\ && + \sum _{K_{i+1} = 1}^{K_i-1} \def\eqcellsep{&}\begin{pmatrix} K_i \\ K_{i+1} \end{pmatrix} \sum _{K_{i+2} = 1}^{K_{i+1}-1} \def\eqcellsep{&}\begin{pmatrix} K_{i+1} \\ K_{i+2} \end{pmatrix} \sum _{|m_{i+1} \rangle \notin D(n)}\nonumber\\ && \sum _{|m_{i+2} \rangle \notin D(n)} \frac{\widetilde{\mathcal {H}}_{m_i, m_{i+1}}^{K_i - K_{i+1}}}{\varepsilon _n^{\,} - \varepsilon _{m_i}^{\,}} \, \frac{\widetilde{\mathcal {H}}_{m_{i+1}, m_{i+2}}^{K_{i+1} - K_{i+2}}}{\varepsilon _n^{\,} - \varepsilon _{m_{i+1}}^{\,}} \, \mathcal {O}_{m_{i+2},n}^{K_{i+2}} \nonumber \\ && +(*) \sum _{K_{i+1} = 1}^{K_i-1} \def\eqcellsep{&}\begin{pmatrix} K_i \\ K_{i+1} \end{pmatrix} \sum _{|\tilde{n} \rangle \in D(n)} \frac{\mathcal {H}_{m_i,\tilde{n}}^{K_i - K_{i+1}}}{\varepsilon _n^{\,} - \varepsilon _{m_i}^{\,}} \, \mathcal {R}_{\tilde{n},n}^{K_{i+1}} \nonumber \\ && + \sum _{K_{i+1} = 1}^{K_i-1} \def\eqcellsep{&}\begin{pmatrix} K_i \\ K_{i+1} \end{pmatrix} \sum _{K_{i+2} = 1}^{K_{i+1}-1} \def\eqcellsep{&}\begin{pmatrix} K_{i+1} \\ K_{i+2} \end{pmatrix} \sum _{|m_{i+1} \rangle \notin D(n)}\nonumber\\ &&\sum _{|\tilde{n} \rangle \in D(n)} \frac{\widetilde{\mathcal {H}}_{m_i, m_{i+1}}^{K_i - K_{i+1}}}{\varepsilon _n^{\,} - \varepsilon _{m_i}^{\,}} \, \frac{\mathcal {H}_{m_{i+1},\tilde{n}}^{K_{i+1} - K_{i+2}}}{\varepsilon _n^{\,} - \varepsilon _{m_{i+1}}^{\,}} \, \mathcal {R}_{\tilde{n},n}^{K_{i+2}} \end{eqnarray}$$$

We can expand this further by substituting the expression for $$\mathcal {O}_{m_{i+2},n}^{K_{i+2}}$$ based on Equation (B6). That can be further expanded by substituting the expression for $$\mathcal {O}_{m_{i+3},n}^{K_{i+3}}$$, once again based on Equation (B6). We cannot, however, continue indefinitely expanding like this. The reason is that we have 1 ⩽ K_{j + 1} ⩽ K_j − 1, which, on repeated application, starting from K_{i + l} up to K_i, gives 1 ⩽ K_{i + l} ⩽ K_{i + l − 1} − 1 ⩽ K_{i + l − 2} − 2 ⩽ … ⩽ K_i − l. This means that the maximum allowed l itself equals K_i − 1 (since 1 ⩽ K_{i + l} ⩽ K_i − l). This also gives the upper limit of the recursive expansion. The full expanded form then reads B9 $$$\begin{eqnarray} \mathcal {O}_{m_i,n}^{K_i} &= & \frac{\mathcal {H}_{m_i,n}^{K_i}}{\varepsilon _n^{\,} - \varepsilon _{m_i}^{\,}} + \sum _{K_{i+1} = 1}^{K_i-1} \sum _{|m_{i+1} \rangle \notin D(n)} \def\eqcellsep{&}\begin{pmatrix} K_i \\ K_{i+1} \end{pmatrix} \frac{\widetilde{\mathcal {H}}_{m_i, m_{i+1}}^{K_i - K_{i+1}}}{\varepsilon _n^{\,} - \varepsilon _{m_i}^{\,}} \, \frac{\mathcal {H}_{m_{i+1},n}^{K_{i+1}}}{\varepsilon _n^{\,} - \varepsilon _{m_{i+1}}^{\,}} \nonumber \\ && +\ \sum _{K_{i+1} = 1}^{K_i-1} \sum _{K_{i+2} = 1}^{K_{i+1}-1} \sum _{|m_{i+1} \rangle \notin D(n)} \sum _{|m_{i+2} \rangle \notin D(n)} \def\eqcellsep{&}\begin{pmatrix} K_i \\ K_{i+1} \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} K_{i+1} \\ K_{i+2} \end{pmatrix}\nonumber\\ && \frac{\widetilde{\mathcal {H}}_{m_i, m_{i+1}}^{K_i - K_{i+1}}}{\varepsilon _n^{\,} - \varepsilon _{m_i}^{\,}} \, \frac{\widetilde{\mathcal {H}}_{m_{i+1}, m_{i+2}}^{K_{i+1} - K_{i+2}}}{\varepsilon _n^{\,} - \varepsilon _{m_{i+1}}^{\,}} \, \frac{\mathcal {H}_{m_{i+2},n}^{K_{i+2}}}{\varepsilon _n^{\,} - \varepsilon _{m_{i+2}}^{\,}} \nonumber \\ && + \cdots \nonumber \\ && + \cdots \nonumber \\ && + \sum _{K_{i+1} = 1}^{K_i - 1} \sum _{K_{i+2} = 1}^{K_{i+1} - 1} \cdots \sum _{K_{i+l} = 1}^{K_{i+l-1} - 1} \sum _{|m_{i+1} \rangle \ne |n \rangle } \sum _{|m_{i+2} \rangle \ne |n \rangle } \cdots \sum _{|m_{i+l} \rangle \ne |n \rangle }\nonumber\\ && {\left[ \def\eqcellsep{&}\begin{pmatrix} K_i \\ K_{i+1} \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} K_{i+1} \\ K_{i+2} \end{pmatrix} \cdots \def\eqcellsep{&}\begin{pmatrix} K_{i+l-1} \\ K_{i+l} \end{pmatrix} \right.}\nonumber \\ && {\left.\times \frac{\widetilde{\mathcal {H}}_{m_i, m_{i+1}}^{K_i - K_{i+1}}}{\varepsilon _n^{\,} - \varepsilon _{m_i}^{\,}} \times \frac{\widetilde{\mathcal {H}}_{m_{i+1}, m_{i+2}}^{K_{i+1} - K_{i+2}}}{\varepsilon _n^{\,} - \varepsilon _{m_{i+1}}^{\,}} \times \cdots \times \frac{\widetilde{\mathcal {H}}_{m_{i+l-1}, m_{i+l}}^{K_{i+l-1} - K_{i+l}}}{\varepsilon _n^{\,} - \varepsilon _{m_{i+l-1}}^{\,}} \, \frac{\mathcal {H}_{m_{i+l},n}^{K_{i+l}}}{\varepsilon _n^{\,} - \varepsilon _{m_{i+l}}^{\,}} \right]} \nonumber \\ && + (*) \sum _{K_{i+1} = 1}^{K_i-1} \sum _{|\tilde{n} \rangle \in D(n)} \def\eqcellsep{&}\begin{pmatrix} K_i \\ K_{i+1} \end{pmatrix} \frac{\mathcal {H}_{m_i,\tilde{n}}^{K_i - K_{i+1}}}{\varepsilon _n^{\,} - \varepsilon _{m_i}^{\,}} \, \mathcal {R}_{\tilde{n},n}^{K_{i+1}} \nonumber \\ && +\ \sum _{K_{i+1} = 1}^{K_i-1} \sum _{K_{i+2} = 1}^{K_{i+1}-1} \sum _{|m_{i+1} \rangle \ne |n \rangle } \sum _{|\tilde{n} \rangle \in D(n)} \def\eqcellsep{&}\begin{pmatrix} K_i \\ K_{i+1} \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} K_{i+1} \\ K_{i+2} \end{pmatrix}\nonumber\\ && \frac{\widetilde{\mathcal {H}}_{m_i, m_{i+1}}^{K_i - K_{i+1}}}{\varepsilon _n^{\,} - \varepsilon _{m_i}^{\,}} \, \frac{\mathcal {H}_{m_{i+1},\tilde{n}}^{K_{i+1} - K_{i+2}}}{\varepsilon _n^{\,} - \varepsilon _{m_{i+1}}^{\,}} \, \mathcal {R}_{\tilde{n},n}^{K_{i+2}} \nonumber \\ && + \cdots \nonumber \\ && + \cdots \nonumber \\ && + \sum _{K_{i+1} = 1}^{K_i - 1} \cdots \sum _{K_{i+l} = 1}^{K_{i+l-1} - 1} \sum _{|m_{i+1} \rangle \ne |n \rangle } \cdots \sum _{|m_{i+l-1} \rangle \ne |n \rangle }\nonumber\\ && \sum _{|\tilde{n} \rangle \in D(n)} {\left[ \def\eqcellsep{&}\begin{pmatrix} K_i \\ K_{i+1} \end{pmatrix} \cdots \def\eqcellsep{&}\begin{pmatrix} K_{i+l-1} \\ K_{i+l} \end{pmatrix} \times \frac{\widetilde{\mathcal {H}}_{m_i, m_{i+1}}^{K_i - K_{i+1}}}{\varepsilon _n^{\,} - \varepsilon _{m_i}^{\,}}\right.}\nonumber \\ && {\left. \times \frac{\widetilde{\mathcal {H}}_{m_{i+1}, m_{i+2}}^{K_{i+1} - K_{i+2}}}{\varepsilon _n^{\,} - \varepsilon _{m_{i+1}}^{\,}} \times \cdots \times \frac{\widetilde{\mathcal {H}}_{m_{i+l-2}, m_{i+l-1}}^{K_{i+l-2} - K_{i+l-1}}}{\varepsilon _n^{\,} - \varepsilon _{m_{i+l-2}}^{\,}} \frac{\mathcal {H}_{m_{i+l-1},\tilde{n}}^{K_{i+l-1} - K_{i+l}}}{\varepsilon _n^{\,} - \varepsilon _{m_{i+l-1}}^{\,}} \right]} \mathcal {R}_{\tilde{n},n}^{K_{i+l}} \end{eqnarray}$$$where the sequence of terms ends with maximum l equalling K_i − 1.

We now proceed to compute the N^th derivative of band n by taking inner product of Equation (B2) with |n⟩, rearranging terms, introducing resolution of identity and separating out the projector to the degenerate subspace to obtain B10a $$$\begin{eqnarray} \partial ^{N} \varepsilon _{n}^{\,} &= & \langle n | \partial ^N H | n \rangle + \sum _{K_1=1}^{N-1} \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \langle n | \partial ^{N - K_1} H | \partial ^{K_1} n \rangle\nonumber\\ && -\ \sum _{K_1=1}^{N-1} \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \partial ^{N - K_1} \varepsilon _n^{\,} \langle n | \partial ^{K_1} n \rangle \end{eqnarray}$$$B10b $$$\begin{eqnarray} &= & \langle n | \partial ^N H | n \rangle + \sum _{K_1=1}^{N-1} \sum _{|m_1 \rangle \notin D(n)} \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \langle n | \partial ^{N - K_1} H | m_1^{\,} \rangle \langle m_1^{\,} | \partial ^{K_1} n \rangle \nonumber \\ && -\ \sum _{K_1=1}^{N-1} \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \partial ^{N - K_1} \varepsilon _n^{\,} \langle n | \partial ^{K_1} n \rangle \nonumber \\ &&+\ \sum _{K_1=1}^{N-1} \sum _{|\tilde{n} \rangle \in D(n)} \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \langle n | \partial ^{N-K_1} H | \tilde{n} \rangle \langle \tilde{n} | \partial ^{K_1} n \rangle \end{eqnarray}$$$ In the symbolic notation, the above reads B11 $$$\begin{eqnarray} \partial ^{N} \varepsilon _{n}^{\,} &= & \mathcal {H}_{n,n}^{N} + \sum _{K_1=1}^{N-1} \sum _{|m_1 \rangle \notin D(n)} \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \mathcal {H}_{n,m_1}^{N- K_1} \, \mathcal {O}_{m_1,n}^{K_1}\nonumber \\ && -\ \sum _{K_1=1}^{N-1} \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \partial ^{N - K_1} \varepsilon _n^{\,} \, \mathcal {R}^{K_1}_{n,n} \nonumber \\ && +\ \sum _{K_1=1}^{N-1} \sum _{\tilde{n} \notin D(n)} \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \mathcal {H}^{N-K_1}_{n,\tilde{n}} \, \mathcal {R}^{K_1}_{\tilde{n},n} \end{eqnarray}$$$We now substitute for $$O_{m_1}^{K_1}$$: B12 $$$\begin{eqnarray} \partial ^{N} \varepsilon _{n}^{\,} &= & \mathcal {H}_{n,n}^{N} + \sum _{K_1=1}^{N-1} \sum _{|m_1 \rangle \notin D(n)} \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \mathcal {H}_{n,m_1}^{N - K_1} \frac{\mathcal {H}_{m_1,n}^{K_1}}{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \nonumber \\ && +\ \sum _{K_1=1}^{N-1} \sum _{K_{2} = 1}^{K_1-1} \sum _{|m_1 \rangle \ne |n \rangle } \sum _{|m_{2} \rangle \ne |n \rangle } \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} K_1 \\ K_{2} \end{pmatrix}\nonumber \\ && \mathcal {H}_{n,m_1}^{N - K_1} \frac{\widetilde{\mathcal {H}}_{m_1, m_{2}}^{K_1 - K_{2}}}{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \frac{\mathcal {H}_{m_{2},n}^{K_{2}}}{\varepsilon _n^{\,} - \varepsilon _{m_{2}}^{\,}} \nonumber \\ && + \cdots \nonumber \\ && + \cdots \nonumber \\ && + \sum _{K_1=1}^{N-1} \sum _{K_{2} = 1}^{K_1-1} \cdots \sum _{K_{N-1} = 1}^{K_{N-2} - 1} \sum _{|m_1 \rangle \ne |n \rangle } \sum _{|m_{2} \rangle \ne |n \rangle }\nonumber \\ && \cdots \sum _{|m_{N-1} \rangle \ne |n \rangle } \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} K_1 \\ K_{2} \end{pmatrix} \cdots \def\eqcellsep{&}\begin{pmatrix} K_{N-2} \\ K_{N-1} \end{pmatrix} \nonumber \\ && {\Bigg[}\times \mathcal {H}_{n,m_1}^{N - K_1} \left(\frac{\widetilde{\mathcal {H}}_{m_1, m_{2}}^{K_1 - K_{2}}}{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \times \frac{\widetilde{\mathcal {H}}_{m_{2}, m_{3}}^{K_{2} - K_{3}}}{\varepsilon _n^{\,} - \varepsilon _{m_{2}}^{\,}} \times \cdots\right.\nonumber \\ &&\left. \times \frac{\widetilde{\mathcal {H}}_{m_{N-2}, m_{N-1}}^{K_{N-2} - K_{N-1}}}{\varepsilon _n^{\,} - \varepsilon _{m_{N-2}}^{\,}} \right) \frac{\mathcal {H}_{m_{N-1},n}^{K_{N-1}}}{\varepsilon _n^{\,} - \varepsilon _{m_{N-1}}^{\,}} \Bigg] \nonumber \\ && -\ (*1) \sum _{K_1=1}^{N-1} \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \partial ^{N - K_1} \varepsilon _n^{\,} \, \mathcal {R}^{K_1}_{n,n} \nonumber \\ && +\ (*2) \sum _{K_1=1}^{N-1} \sum _{|\tilde{n} \rangle \in D(n)} \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \mathcal {H}^{N-K_1}_{n,\tilde{n}} \, \mathcal {R}^{K_1}_{\tilde{n},n} \nonumber \\ && +\ \sum _{K_1=1}^{N-1} \sum _{K_{2} = 1}^{K_1-1} \sum _{|m_1 \rangle \notin D(n)} \sum _{|\tilde{n} \rangle \in D(n)} \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} K_1 \\ K_{2} \end{pmatrix} \mathcal {H}_{n,m_1}^{N - K_1} \frac{\mathcal {H}_{m_1,\tilde{n}}^{K_1 - K_{2}}}{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \mathcal {R}^{K_2}_{\tilde{n},n} \nonumber \hspace*{-20pt}\\ && + \cdots \nonumber \\ && + \cdots \nonumber \\ && + \sum _{K_1=1}^{N-1} \sum _{K_{2} = 1}^{K_1-1} \cdots \sum _{K_{N-1} = 1}^{K_{N-2} - 1} \sum _{|m_1 \rangle \notin D(n)} \cdots \sum _{|m_{N-2} \rangle \notin D(n)} \sum _{|\tilde{n} \rangle \in D(n)}\nonumber \\ && \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} K_1 \\ K_{2} \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} K_2 \\ K_{3} \end{pmatrix} \cdots \def\eqcellsep{&}\begin{pmatrix} K_{N-2} \\ K_{N-1} \end{pmatrix} \nonumber \\ && {\Bigg[}\times\ \mathcal {H}_{n,m_1}^{N - K_1} {\left(\frac{\widetilde{\mathcal {H}}_{m_1, m_{2}}^{K_1 - K_{2}}}{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \times \frac{\widetilde{\mathcal {H}}_{m_{2}, m_{3}}^{K_{2} - K_{3}}}{\varepsilon _n^{\,} - \varepsilon _{m_{2}}^{\,}} \times \cdots\right.}\nonumber \\ && {\left.\qquad \times \frac{\widetilde{\mathcal {H}}_{m_{N-3}, m_{N-2}}^{K_{N-3} - K_{N-2}}}{\varepsilon _n^{\,} - \varepsilon _{m_{N-3}}^{\,}} \right)} \frac{\mathcal {H}_{m_{N-2},\tilde{n}}^{K_{N-2} - K_{N-1}}}{\varepsilon _n^{\,} - \varepsilon _{m_{N-2}}^{\,}} \Bigg] \mathcal {R}^{K_{N-1}}_{\tilde{n},n} \hspace*{-20pt} \end{eqnarray}$$$Theorem: The term (*1) cancels with terms from (2*) onwards so that the band derivative takes the explicitly calculable form: B13 $$$\begin{eqnarray} \partial ^{N} \varepsilon _{n}^{\,} &= & \mathcal {H}_{n,n}^{N} + \sum _{K_1=1}^{N-1} \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \sum _{|m_1 \rangle \notin D(n)} \mathcal {H}_{n,m_1}^{N - K_1} \frac{\mathcal {H}_{m_1,n}^{K_1}}{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \nonumber \\ && +\ \sum _{K_1=1}^{N-1} \sum _{K_{2} = 1}^{K_1-1} \sum _{|m_1 \rangle \notin D(n)} \sum _{|m_{2} \rangle \notin D(n)} \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} K_1 \\ K_{2} \end{pmatrix}\nonumber\\ && \mathcal {H}_{n,m_1}^{N - K_1} \frac{\widetilde{\mathcal {H}}_{m_1, m_{2}}^{K_1 - K_{2}}}{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \frac{\mathcal {H}_{m_{2},n}^{K_{2}}}{\varepsilon _n^{\,} - \varepsilon _{m_{2}}^{\,}} \nonumber \\ && + \cdots \nonumber \\ && + \cdots \nonumber \\ && + \sum _{K_1=1}^{N-1} \sum _{K_{2} = 1}^{K_1-1} \cdots \sum _{K_{N-1} = 1}^{K_{N-2} - 1} \sum _{|m_1 \rangle \notin D(n)} \sum _{|m_{2} \rangle \notin D(n)} \cdots \sum _{|m_{N-1} \rangle \notin D(n)} \nonumber\\ && \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} K_1 \\ K_{2} \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} K_2 \\ K_{3} \end{pmatrix} \cdots \def\eqcellsep{&}\begin{pmatrix} K_{N-2} \\ K_{N-1} \end{pmatrix} \nonumber \\ && {\Bigg[} \times\ \mathcal {H}_{n,m_1}^{N - K_1} {\left(\frac{\widetilde{\mathcal {H}}_{m_1, m_{2}}^{K_1 - K_{2}}}{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \times \frac{\widetilde{\mathcal {H}}_{m_{2}, m_{3}}^{K_{2} - K_{3}}}{\varepsilon _n^{\,} - \varepsilon _{m_{2}}^{\,}} \times \cdots\right.}\nonumber\\ && {\left. \times \frac{\widetilde{\mathcal {H}}_{m_{N-2}, m_{N-1}}^{K_{N-2} - K_{N-1}}}{\varepsilon _n^{\,} - \varepsilon _{m_{N-2}}^{\,}} \right)} \frac{\mathcal {H}_{m_{N-1},n}^{K_{N-1}}}{\varepsilon _n^{\,} - \varepsilon _{m_{N-1}}^{\,}} \Bigg] \end{eqnarray}$$$Furthermore, for any $$|\tilde{n} \rangle$$ distinct from, but degenerate to |n⟩, the following is satisfied B14 $$$\begin{eqnarray} 0 &= & \mathcal {H}_{\tilde{n},n}^{N} + \sum _{K_1=1}^{N-1} \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \sum _{|m_1 \rangle \notin D(n)} \mathcal {H}_{\tilde{n},m_1}^{N - K_1} \frac{\mathcal {H}_{m_1,n}^{K_1}}{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \nonumber \\ && +\ \sum _{K_1=1}^{N-1} \sum _{K_{2} = 1}^{K_1-1} \sum _{|m_1 \rangle \notin D(n)} \sum _{|m_{2} \rangle \notin D(n)} \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} K_1 \\ K_{2} \end{pmatrix}\nonumber\\ &&\ \mathcal {H}_{\tilde{n},m_1}^{N - K_1} \frac{\widetilde{\mathcal {H}}_{m_1, m_{2}}^{K_1 - K_{2}}}{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \frac{\mathcal {H}_{m_{2},n}^{K_{2}}}{\varepsilon _n^{\,} - \varepsilon _{m_{2}}^{\,}} \nonumber \\ && + \cdots \nonumber \\ && + \cdots \nonumber \\ && + \sum _{K_1=1}^{N-1} \sum _{K_{2} = 1}^{K_1-1} \cdots \sum _{K_{N-1} = 1}^{K_{N-2} - 1} \sum _{|m_1 \rangle \notin D(n)} \sum _{|m_{2} \rangle \notin D(n)} \cdots \sum _{|m_{N-1} \rangle \notin D(n)}\nonumber\\ && \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} K_1 \\ K_{2} \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} K_2 \\ K_{3} \end{pmatrix} \cdots \def\eqcellsep{&}\begin{pmatrix} K_{N-2} \\ K_{N-1} \end{pmatrix} \nonumber \\ && {\Bigg[} \times\ \mathcal {H}_{\tilde{n},m_1}^{N - K_1} {\left(\frac{\widetilde{\mathcal {H}}_{m_1, m_{2}}^{K_1 - K_{2}}}{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \times \frac{\widetilde{\mathcal {H}}_{m_{2}, m_{3}}^{K_{2} - K_{3}}}{\varepsilon _n^{\,} - \varepsilon _{m_{2}}^{\,}} \times \ldots\right.}\nonumber\\ &&{\left. \times \frac{\widetilde{\mathcal {H}}_{m_{N-2}, m_{N-1}}^{K_{N-2} - K_{N-1}}}{\varepsilon _n^{\,} - \varepsilon _{m_{N-2}}^{\,}} \right)} \frac{\mathcal {H}_{m_{N-1},n}^{K_{N-1}}}{\varepsilon _n^{\,} - \varepsilon _{m_{N-1}}^{\,}} \Bigg] \end{eqnarray}$$$

Proof: It is easy to see that this form for $$\partial ^N \varepsilon _n^{\,}$$ is satisfied for N = 1, 2 and 3, by explicit computation and comparison with Equation (11). As for Equation (B14), its RHS has the same form as that of $$\partial ^N \varepsilon _n^{\,}$$, with the difference that the leftmost bra is $$\langle \tilde{n}|$$ instead of ⟨n|. As Equations A5 and A12 show, this too holds for N = 1, 2. For higher values of N, we prove both the results by the principle of mathematical induction. Let us begin with $$\partial ^N \varepsilon _n^{\,}$$. Assume that the expressions above are satisfied for $$\partial ^{M} \varepsilon _n^{\,}$$ for 1 ⩽ M ⩽ N − 1. To prove that they are satisfied for M = N, we have to show that the terms from (2*) cancel with (1*) in Equation (B12). Now (2*) is a sum of several individual terms which take the form B15 $$$\begin{eqnarray} &&\hspace*{-6pt}\def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} K_1 \\ K_{2} \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} K_2 \\ K_{3} \end{pmatrix} \ldots \def\eqcellsep{&}\begin{pmatrix} K_{l-1} \\ K_{l} \end{pmatrix} \nonumber\\ &&\hspace*{-6pt} \times \mathcal {H}_{n,m_1}^{N - K_1} {\left(\frac{\widetilde{\mathcal {H}}_{m_1, m_{2}}^{K_1 - K_{2}}}{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \times \frac{\widetilde{\mathcal {H}}_{m_{2}, m_{3}}^{K_{2} - K_{3}}}{\varepsilon _n^{\,} - \varepsilon _{m_{2}}^{\,}} \times \ldots \times \frac{\widetilde{\mathcal {H}}_{m_{l-2}, m_{l-1}}^{K_{l-2} - K_{l-1}}}{\varepsilon _n^{\,} - \varepsilon _{m_{l-2}}^{\,}} \right)} \frac{\mathcal {H}_{m_{l-1},\tilde{n}}^{K_{l-1} - K_{l}}}{\varepsilon _n^{\,} - \varepsilon _{m_{l-1}}^{\,}} \bigg] \mathcal {R}^{K_{l}}_{\tilde{n},n}\nonumber\hspace*{-10pt}\\ \end{eqnarray}$$$for all possible ordered sequences of l natural numbers {K₁, K₂, …, K_l} such that N > K₁ > K₂ > … > K_l ⩾ 1 and 1 ⩽ l ⩽ N − 1. The fact that every such ordered sequence occurs in the sum, and that too only once is ensured by the cascaded K_i sums in Equation (B12) B16 $$$\begin{equation} \sum _{K_1=1}^{N-1} \sum _{K_{2} = 1}^{K_1-1} \ldots \sum _{K_{l} = 1}^{K_{l-1} - 1} \end{equation}$$$Let us denote by A_l, the set of all such sequences of length l (with 1 ⩽ l ⩽ N − 1). That is, A_l = Set of all ordered sequences {K₁, K₂, …, K_l} such that $$K_i \in \mathbb {N}$$ and N > K₁ > K₂ > … > K_l ⩾ 1. It is then easy to see that B17 $$$\begin{equation} \sum _{K_1=1}^{N-1} \sum _{K_{2} = 1}^{K_1-1} \ldots \sum _{K_{l} = 1}^{K_{l-1} - 1} \equiv \sum _{\lbrace K_1, \ldots , K_l\rbrace \in A_l} \end{equation}$$$

Let us define another set B_l which is the set of all ordered sequences of l natural numbers $$\lbrace K_l, \widetilde{K}_1, \widetilde{K}_2, \ldots \widetilde{K}_{l-1} \rbrace$$ such that N − l ⩾ K_l ⩾ 1 and $$N - K_l > \widetilde{K}_1 > \widetilde{K}_2 > \ldots > \widetilde{K}_{l-1} \geqslant 1$$. It is also easy to see that B18 $$$\begin{equation} \sum _{K_l=1}^{N-l} \sum _{\widetilde{K}_{1} = 1}^{N-K_1-1} \sum _{\widetilde{K}_{2} = 1}^{\widetilde{K}_1-1} \ldots \sum _{\widetilde{K}_{l-1} = 1}^{\widetilde{K}_{l-2} - 1} \equiv \sum _{\lbrace K_l, \widetilde{K}_1, \ldots \widetilde{K}_{l-1}\rbrace \in B_l} \end{equation}$$$

Lemma: There exists a bijection between A_l and B_l for each l.

Proof: Let us define F: A_l → B_l, with F(K₁, K₂, …, K_{l − 1}, K_l) = {K_l, K₁ − K_l, K₂ − K_l, …, K_{l − 1} − K_l}. Since N > K₁ > K₂ > … > K_l ⩾ 1, it is easy to see that the range is indeed contained in B_l since N − l ⩾ K_l ⩾ 1 and N − K_l > K₁ − K_l > K₂ − K_l > … > K_{l − 1} − K_l > 0. (The last inequality K_{l − 1} − K_l > 0 is equivalent to K_{l − 1} − K_l ⩾ 1 since we are dealing with natural numbers).

Let us define the function G: B_l → A_l with $$G(K_l, \widetilde{K}_1, \widetilde{K}_2, \ldots \widetilde{K}_{l-1}) = \lbrace \widetilde{K}_1 + K_l, \widetilde{K}_2 + K_l, \ldots ,\widetilde{K}_{l-1} + K_l, K_l \rbrace$$. Since N − l ⩾ K_l ⩾ 1 and $$N - K_l > \widetilde{K}_1 > \widetilde{K}_2 > \ldots > \widetilde{K}_{l-1} \geqslant 1$$, it is easy to see that $$N > \widetilde{K}_1 + K_l > \widetilde{K}_2 + K_l > \ldots > \widetilde{K}_{l-1} + K_l > K_l \geqslant 1$$ so that the range of this function is contained in A_l. Since F and G are clearly inverses of each other, they are both bijections. In particular, this means that iterating over each member of A_l and evaluating a function $$\psi _l^{\,}$$ defined on A_l is equivalent to iterating over each member of B_l: B19 $$$\begin{align} {{ \def\eqcellsep{&}\begin{array}{rcl} && \displaystyle\sum _{l=1}^{N-1} \sum _{\lbrace K_1, \ldots , K_l\rbrace \in A_l} \psi _l^{\,}(K_1, \ldots , K_l) \equiv \sum _{l=1}^{N-1} \sum _{\lbrace K_l, \widetilde{K}_1, \ldots , \widetilde{K}_{l-1}\rbrace \in B_l } \psi _l^{\,} (G(K_l, \widetilde{K}_1, \ldots , \widetilde{K}_{l-1})) \nonumber \\[12pt] &\Rightarrow & \displaystyle\sum _{l=1}^{N-1} \sum _{K_1=1}^{N-1} \sum _{K_{2} = 1}^{K_1-1} \ldots \sum _{K_{l} = 1}^{K_{l-1} - 1} \psi _l^{\,}(K_1, K_2, \ldots , K_l) \nonumber \\[12pt] & =& \displaystyle\sum _{l=1}^{N-1} \sum _{K_l=1}^{N-l} \sum _{\widetilde{K}_{1} = 1}^{N - K_1 - 1} \sum _{\widetilde{K}_{2} = 1}^{\widetilde{K}_1 - 1} \ldots \sum _{\widetilde{K}_{l-1} = 1}^{\widetilde{K}_{l-2} - 1} \psi _l^{\,}(G(K_l, \widetilde{K}_1, \widetilde{K}_2, \ldots , \widetilde{K}_{l-1})) \\[-28pt] \end{array} }}\nonumber\\ \end{align}$$$The sums on LHS and RHS exhaustively go through all the elements in A_l and B_l respectively, as seen from Equations B17 and B18. The function g is needed on the RHS since, for each element of B_l, we need to feed its unique counterpart in A_l as an input to the function $$\psi _l^{\,}$$, to sum over.

To be able to interchange the K_l and l sums in the RHS, let us define C to be the set of all ordered pairs {l, K_l} with 1 ⩽ l ⩽ N − 1 and 1 ⩽ K_l ⩽ N − l and D to be the set of all ordered pairs $$\lbrace \mathfrak {K}_1, \mathfrak {l} \rbrace$$ with $$1 \leqslant \mathfrak {K}_1 \leqslant N - 1$$ and $$1 \leqslant \mathfrak {l} \leqslant N - \mathfrak {K}_1$$. As before, we define functions $$\mathcal {F}: C \rightarrow D$$ with $$\mathcal {F} (l, K_l) = \lbrace K_l, l \rbrace$$ and $$\mathcal {G}: D \rightarrow C$$ with $$\mathcal {G} (\mathfrak {K}_1, \mathfrak {l}) = \lbrace \mathfrak {l}, \mathfrak {K}_1 \rbrace$$. It is easy to check that $$\text{Range} [\mathcal {F}] \subseteq D$$, $$\text{Range} [\mathcal {G}] \subseteq C$$, and the functions are inverses of each other, allowing us to interchange the original sums: B20 $$$\begin{equation} \sum _{l=1}^{N-1} \sum _{K_l=1}^{N-l} \equiv \sum _{\mathfrak {K}_1=1}^{N-1} \sum _{\mathfrak {l}=1}^{N-\mathfrak {K}_1} \end{equation}$$$This is geometrically depicted in Figure B1, where the grid points in the left panel are being summed over exhaustively by the sum on LHS in Equation (B20). Clearly, the sum on RHS is just a different but equivalent approach to summing over the same set of points, depicted in the right panel. Thus, we have B21 $$$\begin{align} {{ \def\eqcellsep{&}\begin{array}{rcl} &&\displaystyle\sum _{l=1}^{N-1} \sum _{K_1=1}^{N-1} \sum _{K_{2} = 1}^{K_1-1} \ldots \sum _{K_{l} = 1}^{K_{l-1} - 1} \psi _l^{\,}(K_1, K_2, \ldots , K_l) \nonumber \\[12pt] & =& \displaystyle\sum _{\mathfrak {K}_1=1}^{N-1} \sum _{\mathfrak {l}=1}^{N-\mathfrak {K}_1} \sum _{\widetilde{K}_{1} = 1}^{N - \mathfrak {K}_1 - 1} \sum _{\widetilde{K}_{2} = 1}^{\widetilde{K}_1 - 1} \ldots \sum _{\widetilde{K}_{\mathfrak {l}-1} = 1}^{\widetilde{K}_{\mathfrak {l}-2} - 1} \psi _\mathfrak {l}^{\,}(G(\mathfrak {K}_1, \widetilde{K}_1, \widetilde{K}_2, \ldots , \widetilde{K}_{\mathfrak {l}-1})) \nonumber \\[12pt] & =& \displaystyle\sum _{\mathfrak {K}_1=1}^{N-1} \sum _{\mathfrak {l}=1}^{N-\mathfrak {K}_1} \sum _{\widetilde{K}_{1} = 1}^{N - \mathfrak {K}_1 - 1} \sum _{\widetilde{K}_{2} = 1}^{\widetilde{K}_1 - 1} \ldots \sum _{\widetilde{K}_{\mathfrak {l}-1} = 1}^{\widetilde{K}_{\mathfrak {l}-2} - 1} \psi _\mathfrak {l}^{\,}(\widetilde{K}_1 + \mathfrak {K}_1, \widetilde{K}_2 + \mathfrak {K}_1, \ldots , \widetilde{K}_{\mathfrak {l}-1} + \mathfrak {K}_1, \mathfrak {K}_1) \end{array} }}\nonumber\\ \end{align}$$$For convenience, we shall henceforth refer to $$\mathfrak {K}_1$$ as K₁ and $$\mathfrak {l}$$ as l. Applying the above results to (2*) in Equation (B12) (suppressing the |m_i⟩∉D(n) and $$|\tilde{n} \rangle \in D(n)$$ sums for the sake convenience), we get B22a $$$\begin{eqnarray} {{ \def\eqcellsep{&}\begin{array}{rcl} && (2*) \nonumber \\ &=& \displaystyle\sum _{l=1}^{N-1} \sum _{K_1=1}^{N-1} \sum _{K_{2} = 1}^{K_1-1} \ldots \sum _{K_{l} = 1}^{K_{l-1} - 1} \left[ \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} K_1 \\ K_{2} \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} K_2 \\ K_{3} \end{pmatrix} \ldots \def\eqcellsep{&}\begin{pmatrix} K_{l-1} \\ K_{l} \end{pmatrix} \right.\nonumber\\[12pt] &&{\left. \times \mathcal {H}_{n,m_1}^{N - K_1} {\left(\displaystyle\frac{\widetilde{\mathcal {H}}_{m_1, m_{2}}^{K_1 - K_{2}}}{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \times \frac{\widetilde{\mathcal {H}}_{m_{2}, m_{3}}^{K_{2} - K_{3}}}{\varepsilon _n^{\,} - \varepsilon _{m_{2}}^{\,}} \times \ldots \times \frac{\widetilde{\mathcal {H}}_{m_{l-2}, m_{l-1}}^{K_{l-2} - K_{l-1}}}{\varepsilon _n^{\,} - \varepsilon _{m_{l-2}}^{\,}} \right)} \displaystyle\frac{\mathcal {H}_{m_{l-1},\tilde{n}}^{K_{l-1} - K_{l}}}{\varepsilon _n^{\,} - \varepsilon _{m_{l-1}}^{\,}} \right]} \mathcal {R}^{K_{l}}_{\tilde{n},n} \end{array} }}\nonumber\hspace*{-10pt}\\ \end{eqnarray}$$$B22b $$$\begin{align} {{ \def\eqcellsep{&}\begin{array}{rcl} & =& \displaystyle\sum _{K_1=1}^{N-1} \sum _{l=1}^{N-K_1} \sum _{\widetilde{K}_{1} = 1}^{N - K_1 - 1} \sum _{\widetilde{K}_{2} = 1}^{\widetilde{K}_1 - 1} \ldots \sum _{\widetilde{K}_{l-1} = 1}^{\widetilde{K}_{l-2} - 1} \\[12pt] && \quad {\Bigg[}\def\eqcellsep{&}\begin{pmatrix} N \\ \widetilde{K}_1 + K_1 \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} \widetilde{K}_1 + K_1 \\ \widetilde{K}_{2} + K_1 \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} \widetilde{K}_2 + K_1 \\ \widetilde{K}_{3} + K_1 \end{pmatrix} \ldots \def\eqcellsep{&}\begin{pmatrix} \widetilde{K}_{l-2} + K_1 \\ \widetilde{K}_{l-1} + K_1 \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} \widetilde{K}_{l-1} + K_1 \\ K_1 \end{pmatrix} \\[12pt] &&{\left. \quad \times \mathcal {H}_{n,m_1}^{N - K_1 - \widetilde{K}_1} {\left(\displaystyle\frac{\widetilde{\mathcal {H}}_{m_1, m_{2}}^{\widetilde{K}_1 - \widetilde{K}_{2}}}{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \times \frac{\widetilde{\mathcal {H}}_{m_{2}, m_{3}}^{\widetilde{K}_{2} - \widetilde{K}_{3}}}{\varepsilon _n^{\,} - \varepsilon _{m_{2}}^{\,}} \times \ldots \times \displaystyle\frac{\widetilde{\mathcal {H}}_{m_{l-2}, m_{l-1}}^{\widetilde{K}_{l-2} - \widetilde{K}_{l-1}}}{\varepsilon _n^{\,} - \varepsilon _{m_{l-2}}^{\,}} \right)} \displaystyle\frac{\mathcal {H}_{m_{l-1},\tilde{n}}^{\widetilde{K}_{l-1}}}{\varepsilon _n^{\,} - \varepsilon _{m_{l-1}}^{\,}} \right]} \mathcal {R}^{K_{1}}_{\tilde{n},n} \end{array} }}\nonumber\hspace*{-10pt}\\ \end{align}$$$ Now, for the last ingredient of the proof, we use the equivalent expression of the product of the combinatorial factors: B23 $$$\begin{align} \def\eqcellsep{&}\begin{pmatrix} N \\ \widetilde{K}_1 + K_1 \end{pmatrix} & \def\eqcellsep{&}\begin{pmatrix} \widetilde{K}_1 + K_1 \\ \widetilde{K}_{2} + K_1 \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} \widetilde{K}_2 + K_1 \\ \widetilde{K}_{3} + K_1 \end{pmatrix} \ldots \def\eqcellsep{&}\begin{pmatrix} \widetilde{K}_{l-2} + K_1 \\ \widetilde{K}_{l-1} + K_1 \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} \widetilde{K}_{l-1} + K_1 \\ K_1 \end{pmatrix} \nonumber \\ & = \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} N - K_1 \\ \widetilde{K}_{1} \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} \widetilde{K}_1 \\ \widetilde{K}_{2} \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} \widetilde{K}_2 \\ \widetilde{K}_{3} \end{pmatrix} \ldots \def\eqcellsep{&}\begin{pmatrix} \widetilde{K}_{l-2} \\ \widetilde{K}_{l-1} \end{pmatrix} \end{align}$$$Collecting all of the results above, we obtain for (2*) B24 $$$\begin{eqnarray} (2*) & =& \sum _{K_1=1}^{N-1} \, \sum _{l=1}^{N-K_1} \, \sum _{\widetilde{K}_{1} = 1}^{N - K_1 - 1} \, \sum _{\widetilde{K}_{2} = 1}^{\widetilde{K}_1 - 1} \cdots \sum _{\widetilde{K}_{l-1} = 1}^{\widetilde{K}_{l-2} - 1}\nonumber\\ &&{\left[ \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} N - K_1 \\ \widetilde{K}_{1} \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} \widetilde{K}_1 \\ \widetilde{K}_{2} \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} \widetilde{K}_2 \\ \widetilde{K}_{3} \end{pmatrix} \cdots \def\eqcellsep{&}\begin{pmatrix} \widetilde{K}_{l-2} \\ \widetilde{K}_{l-1} \end{pmatrix} \right.}\nonumber \\ && \times\ \mathcal {H}_{n,m_1}^{N - K_1 - \widetilde{K}_1} \left(\frac{\widetilde{\mathcal {H}}_{m_1, m_{2}}^{\widetilde{K}_1 - \widetilde{K}_{2}}}{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \times \frac{\widetilde{\mathcal {H}}_{m_{2}, m_{3}}^{\widetilde{K}_{2} - \widetilde{K}_{3}}}{\varepsilon _n^{\,} - \varepsilon _{m_{2}}^{\,}} \times \cdots\right.\nonumber\\ && \left.\times \frac{\widetilde{\mathcal {H}}_{m_{l-2}, m_{l-1}}^{\widetilde{K}_{l-2} - \widetilde{K}_{l-1}}}{\varepsilon _n^{\,} - \varepsilon _{m_{l-2}}^{\,}} \right) \frac{\mathcal {H}_{m_{l-1},\tilde{n}}^{\widetilde{K}_{l-1}}}{\varepsilon _n^{\,} - \varepsilon _{m_{l-1}}^{\,}} \Bigg] \mathcal {R}^{K_{1}}_{\tilde{n},n} \nonumber \\ & =& \sum _{K_1=1}^{N-1} \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \mathcal {R}^{K_{1}}_{\tilde{n},n} {\left[ \mathcal {H}_{n, \tilde{n}}^{N-K_1} + \sum _{\widetilde{K}_1=1}^{N - K_1 -1} \def\eqcellsep{&}\begin{pmatrix} N - K_1 \\ \widetilde{K}_1 \end{pmatrix}\right.}\nonumber\\ && {\left.\sum _{|m_1 \rangle \notin D(n)} \mathcal {H}_{n,m_1}^{N - K_1 - \widetilde{K}_1} \frac{\mathcal {H}_{m_1, \tilde{n}}^{\widetilde{K}_1}}{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}}\right.} \nonumber \\ && + \sum _{\widetilde{K}_1=1}^{N - K_1-1} \, \sum _{\widetilde{K}_{2} = 1}^{\widetilde{K}_1-1} \, \sum _{|m_1 \rangle \notin D(n)} \, \sum _{|m_{2} \rangle \notin D(n)} \def\eqcellsep{&}\begin{pmatrix} N - K_1 \\ \widetilde{K}_1 \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} \widetilde{K}_1 \\ \widetilde{K}_{2} \end{pmatrix} \nonumber\\ && \mathcal {H}_{n,m_1}^{N - K_1 - \widetilde{K}_1} \frac{\widetilde{\mathcal {H}}_{m_1, m_{2}}^{\widetilde{K}_1 - \widetilde{K}_{2}}}{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \frac{\mathcal {H}_{m_{2},\tilde{n}}^{\widetilde{K}_{2}}}{\varepsilon _n^{\,} - \varepsilon _{m_{2}}^{\,}} \nonumber \\ && +\ \cdots \nonumber \\ && +\ \cdots \nonumber \\ && +\ \sum _{\widetilde{K}_1=1}^{N - K_1 -1} \, \sum _{\widetilde{K}_{2} = 1}^{\widetilde{K}_1-1} \cdots \sum _{\widetilde{K}_{N - K_1 -1} = 1}^{\widetilde{K}_{N - K_1 -2} - 1} \, \sum _{|m_1 \rangle \notin D(n)} \, \sum _{|m_{2} \rangle \notin D(n)} \cdots \sum _{|m_{N - K_1 - 1} \rangle \notin D(n)} \nonumber \\ && \Bigg \lbrace\def\eqcellsep{&}\begin{pmatrix} N - K_1 \\ \widetilde{K}_1 \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} \widetilde{K}_1 \\ \widetilde{K}_{2} \end{pmatrix} \def\eqcellsep{&}\begin{pmatrix} \widetilde{K}_2 \\ \widetilde{K}_{3} \end{pmatrix} \cdots \def\eqcellsep{&}\begin{pmatrix} \widetilde{K}_{N - K_1 - 2} \\ \widetilde{K}_{N - K_1 - 1} \end{pmatrix} \nonumber \\ && \times\ \mathcal {H}_{n,m_1}^{N - K_1 - \widetilde{K}_1} {\left(\frac{\widetilde{\mathcal {H}}_{m_1, m_{2}}^{\widetilde{K}_1 - \widetilde{K}_{2}}}{\varepsilon _n^{\,} - \varepsilon _{m_1}^{\,}} \times \frac{\widetilde{\mathcal {H}}_{m_{2}, m_{3}}^{\widetilde{K}_{2} - \widetilde{K}_{3}}}{\varepsilon _n^{\,} - \varepsilon _{m_{2}}^{\,}} \times \ldots \times \frac{\widetilde{\mathcal {H}}_{m_{N - K_1 - 2}, m_{N - K_1 - 1}}^{\widetilde{K}_{N - K_1 - 2} - \widetilde{K}_{N - K_1 - 1}}}{\varepsilon _n^{\,} - \varepsilon _{m_{N - K_1 - 2}}^{\,}} \right)} \nonumber \\ && {\left. \times \frac{\mathcal {H}_{m_{N - K_1 - 1}, \tilde{n}}^{K_{N - K_1 - 1}}}{\varepsilon _n^{\,} - \varepsilon _{m_{N - K_1 - 1}}^{\,}} \Bigg \rbrace \right]} \end{eqnarray}$$$On inspection, we observe that the terms within the square brackets [ ] precisely add up to give $$\partial ^{N - K_1} \varepsilon _n^{\,}$$ when $$|\tilde{n} \rangle = |n \rangle$$ (according to Equation (B13)) and vanish when $$|\tilde{n} \rangle \ne |n \rangle$$ (as per Equation (B14)). This is the case because N − K₁ < N or equivalently N − K₁ ⩽ N − 1, and we have assumed that the form of the band derivative given in Equation (B13) holds up to N − 1. Thus, some of the terms in (2*) vanish while the rest add up to B25 $$$\begin{equation} (2*) = \sum _{K_1 = 1}^{N - 1} \def\eqcellsep{&}\begin{pmatrix} N \\ K_1 \end{pmatrix} \partial ^{N - K_1} \varepsilon _n^{\,} \, \mathcal {R}_{n,n}^{K_1} \end{equation}$$$which exactly cancels out with (1*). As for the proof of Equation (B14) at M = N, we take overlap of Equation (B2) with $$|\tilde{n} \rangle$$ instead of |n⟩, which causes the term containing $$\partial ^{N} \varepsilon _n^{\,}$$ to vanish. The rest of the proof proceeds exactly the same way as for the proof for $$\partial ^N \varepsilon _n^{\,}$$, with the extra assumption that the projection of the derivatives of the Hamiltonian ∂^MH to the degenerate subspace are diagonal and proportional to identity matrix. This ensures that the derivatives for the degenerate bands are also equal. Thus we have proved that the forms given in Equation (B13) and B14 hold for N as long as they hold up from 1 to N − 1. Since they have already been shown to hold for N = 1, 2, it follows that they hold for all N.

[IMAGE OMITTED. SEE PDF]

Acknowledgements

The authors would like to thank Claudio Chamon, Mark Greenaway, Nick Hine, Phil King, Samuel Magorrian, Edgar Abarca Morales, Andreas Rost, and Peter Wahl for their useful discussions. The work was supported by the Engineering and Physical Sciences Research Council (EPSRC) grant EP/T034351/1.

Conflict Of Interest

The authors declare no conflict of interest.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

References

I. Lifshitz, Sov. Phys. JETP 1960, 11, 1130.

L. Van Hove, Phys. Rev. 1953, 89, 1189.

A. Chandrasekaran, A. Shtyk, J. J. Betouras, C. Chamon, Phys. Rev. Res. 2020, 2, [eLocator: 013355].

D. V. Efremov, A. Shtyk, A. W. Rost, C. Chamon, A. P. Mackenzie, J. J. Betouras, Phys. Rev. Lett. 2019, 123, [eLocator: 207202].

A. Shtyk, G. Goldstein, C. Chamon, Phys. Rev. B 2017, 95, [eLocator: 035137].

N. F. Q. Yuan, L. Fu, Phys. Rev. B 2020, 101, [eLocator: 125120].

A. Zervou, G. Goldstein, D. V. Efremov, J. J. Betouras, The fate of density waves in the presence of a higher order van Hove singularity, arXiv:2205.08828.

D. Aoki, G. Seyfarth, A. Pourret, A. Gourgout, A. McCollam, J. A. N. Bruin, Y. Krupko, I. Sheikin, Phys. Rev. Lett. 2016, 116, [eLocator: 037202].

M. E. Barber, F. Lechermann, S. V. Streltsov, S. L. Skornyakov, S. Ghosh, B. J. Ramshaw, N. Kikugawa, D. A. Sokolov, A. P. Mackenzie, C. W. Hicks, I. I. Mazin, Phys. Rev. B 2019, 100, [eLocator: 245139].

S. Benhabib, A. Sacuto, M. Civelli, I. Paul, M. Cazayous, Y. Gallais, M.‐A. Méasson, R. D. Zhong, J. Schneeloch, G. D. Gu, D. Colson, A. Forget, Phys. Rev. Lett. 2015, 114, [eLocator: 147001].

S. N. Khan, D. D. Johnson, Phys. Rev. Lett. 2014, 112, [eLocator: 156401].

A. I. Coldea, S. F. Blake, S. Kasahara, A. A. Haghighirad, M. D. Watson, W. Knafo, E. S. Choi, A. McCollam, P. Reiss, T. Yamashita, M. Bruma, S. C. Speller, Y. Matsuda, T. Wolf, T. Shibauchi, A. J. Schofield, npj Quant. Mater. 2019, 4, 1.

Y. Okamoto, A. Nishio, Z. Hiroi, Phys. Rev. B 2010, 81, [eLocator: 121102].

Y. Sherkunov, A. V. Chubukov, J. J. Betouras, Phys. Rev. Lett. 2018, 121, [eLocator: 097001].

S. Slizovskiy, A. V. Chubukov, J. J. Betouras, Phys. Rev. Lett. 2015, 114, [eLocator: 066403].

I. Stewart, Rep. Prog. Phys. 1982, 45, 185.

E. Yelland, J. Barraclough, W. Wang, K. Kamenev, A. Huxley, Nat. Phys. 2011, 7, 890.

N. F. Yuan, H. Isobe, L. Fu, Nat. Commun. 2019, 10, 1.

Y. Sherkunov, J. J. Betouras, Phys. Rev. B 2018, 98, [eLocator: 205151].

L. Classen, A. V. Chubukov, C. Honerkamp, M. M. Scherer, Phys. Rev. B 2020, 102, [eLocator: 125141].

H. Isobe, L. Fu, Phys. Rev. Res. 2019, 1, [eLocator: 033206].

J. L. McChesney, A. Bostwick, T. Ohta, T. Seyller, K. Horn, J. González, E. Rotenberg, Phys. Rev. Lett. 2010, 104, [eLocator: 136803].

P. Rosenzweig, H. Karakachian, D. Marchenko, K. Küster, U. Starke, Phys. Rev. Lett. 2020, 125, [eLocator: 176403].

M. Kang, S. Fang, J.‐K. Kim, B. R. Ortiz, S. H. Ryu, J. Kim, J. Yoo, G. Sangiovanni, D. Di Sante, B.‐G. Park, C. Jozwiak, A. Bostwick, E. Rotenberg, E. Kaxiras, S. D. Wilson, J.‐H. Park, R. Comin, Nat. Phys. 2022, 1.

Y. Hu, X. Wu, B. R. Ortiz, S. Ju, X. Han, J. Ma, N. C. Plumb, M. Radovic, R. Thomale, S. D. Wilson, A. P. Schnyder, M. Shi, Nat. Commun. 2022, 13, 1.

A. Consiglio, T. Schwemmer, X. Wu, W. Hanke, T. Neupert, R. Thomale, G. Sangiovanni, D. Di Sante, Phys. Rev. B 2022, 105, [eLocator: 165146].

T. Neupert, M. M. Denner, J.‐X. Yin, R. Thomale, M. Z. Hasan, Nat. Phys. 2022, 18, 137.

H. Zhou, L. Holleis, Y. Saito, L. Cohen, W. Huynh, C. L. Patterson, F. Yang, T. Taniguchi, K. Watanabe, A. F. Young, Science 2022, 375, 774.

A. Chandrasekaran, J. J. Betouras, Phys. Rev. B 2022, 105, [eLocator: 075144].

E. K. Kokkinis, G. Goldstein, D. V. Efremov, J. J. Betouras, Phys. Rev. B 2022, 105, [eLocator: 155123].

R. G. Ayoub, Archive for history of exact sciences 1980, 253.

B. Bradlyn, J. Cano, Z. Wang, M. Vergniory, C. Felser, R. J. Cava, B. A. Bernevig, Science 2016, 353, [eLocator: aaf5037].

W. R. Inc., Mathematica, Version 13.1, https://www.wolfram.com/mathematica, Champaign, IL, 2022.

Maplesoft, a division of Waterloo Maple Inc., Waterloo, Ontario, Maple.

F. D. M. Haldane, Phys. Rev. Lett. 1988, 61, 2015.

A. Chandrasekaran, Band utilities, https://github.com/anirudhc‐git/Series‐expansion‐of‐bandstructure, A Mathematica package for analysing bandstructures arising from k‐space models.