# A practical method to detect, analyse and engineer higher order Van Hove
singularities in multi-band Hamiltonians

Anirudh Chandrasekaran  Department of Physics and Centre for the Science of
Materials, Loughborough University, Loughborough LE11 3TU, UK.  Joseph J.
Betouras  Department of Physics and Centre for the Science of Materials,
Loughborough University, Loughborough LE11 3TU, UK.

###### Abstract

We present a practical method to detect, diagnose and engineer higher order
Van Hove singularities in multiband systems, with no restrictions on the
number of bands and 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 generalised extension of the Feynman Hellmann theorem,
which we state and prove. Being fairly general in scope, it also allows us to
incorporate and analyse 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. We demonstrate the use of the method, by applying it to the Haldane
model.

##  I Introduction

Nearly a century after its introduction, modelling 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 two dimensional 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, 4, 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, $\text{Sr}_{2}\text{RuO}_{4}$ and heavy fermion
systems [8, 9, 10, 11, 12, 13, 14, 15, 16, 17]. With the recent advances in
the experimental methods 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 behaviour of
$\mathrm{Sr_{3}Ru_{2}O_{7}}$ [4], twisted bilayer graphene near magic angle,
close to half filling [18, 19, 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, 25,
26, 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 two dimensional and layered materials
(see 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, i.e, 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 $\mathbf{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, 30, 29]. 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:

$\varepsilon({\bf
k})\propto(k_{x}^{\,}+k_{y}^{\,})^{4}-(k_{x}^{\,}-k_{y}^{\,})^{2}$ (1)

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 II where we provide the method
and some general formulae. It is worth emphasizing that the method can be
applied to any Hamiltonian $H(\mathbf{k})$, not just ones arising from tight
binding models. In Sec. III 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 Sec. IV. For clarity and to facilitate the flow, details of the
calculations and the full general proof are presented in the appendix.

##  II 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\cdot
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 ${\bf 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
$\mathbf{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
$\mathbf{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.

###  II.1 Preliminaries

To compute the Taylor expansion to order $N$ at $\mathbf{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}(\mathbf{k})|_{\mathbf{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 $\lambda$, taking the form
$\partial^{M}\varepsilon_{n}^{\,}(\mathbf{k}_{0}^{\,}+\lambda\,\mathbf{k})/\partial\lambda^{M}$
and setting $\lambda$ to zero at the end. By repeated application of the chain
rule for differentiation, it can be shown that

$\frac{\partial^{M}\varepsilon_{n}^{\,}(\mathbf{k}_{0}^{\,}+\lambda\,\mathbf{k})}{\partial\lambda^{M}}\bigg{|}_{\lambda=0}=\sum_{\begin{subarray}{c}1\leqslant
l_{1},l_{2},\cdots,l_{d}<M\\\
l_{1}+l_{2}+\cdots+l_{d}=M\end{subarray}}\frac{M!}{l_{1}^{\,}!\,l_{2}^{\,}!\cdots
l_{d}^{\,}!}\,\frac{\partial^{M}\varepsilon_{n}^{\,}(\mathbf{k})}{\partial
k_{1}^{l_{1}}\,\partial k_{2}^{l_{2}}\cdots\partial
k_{d}^{l_{d}}}\bigg{|}_{\mathbf{k}=\mathbf{k}_{0}}\,k_{1}^{l_{1}}\,k_{2}^{l_{2}}\cdots
k_{d}^{l_{d}}.$ (2)

Thus,
$\frac{1}{M!}\frac{\partial^{M}}{\partial\lambda^{M}}\varepsilon_{n}^{\,}(\mathbf{k}_{0}^{\,}+\lambda\,\mathbf{k})\bigg{|}_{\lambda=0}$
is the sum of _all_ the monomials of degree $M$ in the Taylor expansion of
$\varepsilon_{n}^{\,}(\mathbf{k})$ at $\mathbf{k}_{0}^{\,}$. The full Taylor
expansion to order $N$ is then given by

$\varepsilon_{n}^{\,}(\mathbf{k}_{0}^{\,}+\mathbf{k})=\sum_{M=1}^{N}\frac{1}{M!}\,\frac{\partial^{M}\varepsilon_{n}^{\,}(\mathbf{k}_{0}^{\,}+\lambda\,\mathbf{k})}{\partial\lambda^{M}}\bigg{|}_{\lambda=0}+\mathcal{O}\left(k^{N+1}\right).$
(3)

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.

###  II.2 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
$\mathbf{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 $\mathbf{k}_{0}^{\,}$
to order $N$, of the $n^{\text{th}}$ band (denoted by
$\varepsilon_{n}^{\,}(\mathbf{k})$) arising from the Hamiltonian matrix
$\widehat{H}(\mathbf{k})$. The method proceeds to compute the terms order by
order.

  1. 1.

Compute numerically, all the eigenvalues and eigenvectors of
$\widehat{H}(\mathbf{k})$ at the chosen point $\mathbf{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. 2.

Perform a change of basis from the orbital basis to the eigenbasis _at_
$\mathbf{k}_{0}^{\,}$ by transforming the Hamiltonian:
$\widehat{H}(\mathbf{k})\rightarrow\widehat{E}^{{\dagger}}\,\widehat{H}(\mathbf{k})\,\widehat{E}$.
(Note that we are transforming the Hamiltonian as a function of $\mathbf{k}$,
not just the Hamiltonian at the chosen $\mathbf{k}_{0}^{\,}$. The eigenvectors
that constitute the new basis are however the numerically computed ones at
$\mathbf{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. 3.

Calculate the following matrices analytically by directly differentiating the
basis transformed Hamiltonian:

$\partial\widehat{H}\coloneqq\frac{\partial\widehat{H}(\mathbf{k}_{0}^{\,}+\lambda\,\mathbf{k})}{\partial\lambda}\bigg{|}_{\lambda=0},\partial^{2}\widehat{H}\coloneqq\frac{\partial^{2}\widehat{H}(\mathbf{k}_{0}^{\,}+\lambda\,\mathbf{k})}{\partial\lambda^{2}}\bigg{|}_{\lambda=0},\cdots,\partial^{N}\widehat{H}\coloneqq\frac{\partial^{N}\widehat{H}(\mathbf{k}_{0}^{\,}+\lambda\mathbf{k})}{\partial\lambda^{N}}\bigg{|}_{\lambda=0}.$
(4)

These matrices when divided by the appropriate factorial weights and assembled
together make up the Taylor expansion of the Hamiltonian at
$\mathbf{k}_{0}^{\,}$ to order $N$.

  4. 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 Eq 3:

$\partial\varepsilon_{n}^{\,}\coloneqq\frac{\partial\varepsilon_{n}^{\,}(\mathbf{k}_{0}^{\,}+\lambda\,\mathbf{k})}{\partial\lambda}\bigg{|}_{\lambda=0},\partial^{2}\varepsilon_{n}^{\,}\coloneqq\frac{\partial^{2}\varepsilon_{n}^{\,}(\mathbf{k}_{0}^{\,}+\lambda\,\mathbf{k})}{\partial\lambda^{2}}\bigg{|}_{\lambda=0},\cdots,\partial^{N}\varepsilon_{n}^{\,}\coloneqq\frac{\partial^{N}\varepsilon_{n}^{\,}(\mathbf{k}_{0}^{\,}+\lambda\,\mathbf{k})}{\partial\lambda^{N}}\bigg{|}_{\lambda=0}.$
(5)

  5. 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}^{\,}\neq\varepsilon_{n}^{\,}$ and $\kappa\geqslant 1$:

$f(\kappa,m,n)=\frac{\left[\partial^{\kappa}\widehat{H}\right]_{m,n}^{\,}}{\varepsilon_{n}^{\,}-\varepsilon_{m}^{\,}}+\sum\limits_{\widetilde{\kappa}=1}^{\kappa-1}\begin{pmatrix}\kappa\\\
\widetilde{\kappa}\end{pmatrix}\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),$
(6)

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}$, $\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\rangle$ (including
itself), so that the $|\widetilde{m}\rangle$ sum runs over all the
eigenvectors at $\mathbf{k}_{0}^{\,}$ that are _not_ degenerate to
$|n\rangle$.

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(\kappa,m,n)$ for $\kappa=1,2$

$\displaystyle f(1,m,n)=$
$\displaystyle\frac{[\partial\widehat{H}]_{m,n}}{\varepsilon_{n}^{\,}-\varepsilon_{m}^{\,}},$
(7a) $\displaystyle f(2,m,n)=$
$\displaystyle\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}}.$
(7b)

The sequence of product of $f$’s always terminates when
$\widetilde{\kappa}=1$, and we have in that case
$f\left(\widetilde{\kappa}=1,\widetilde{m},n\right)=[\partial\widehat{H}]_{\widetilde{m},n}^{\,}/(\varepsilon_{n}^{\,}-\varepsilon_{\widetilde{m}}^{\,})$.
For any $\kappa>1$ and given $m$ and $n$, to evaluate $f(\kappa,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}^{\,},\cdots,\partial^{\kappa-1}\varepsilon_{n}^{\,}$.
But since we know the first derivative from Feynman - Hellmann lemma,
$\partial\varepsilon_{n}^{\,}=\left[\partial\widehat{H}\right]_{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):

$\partial^{M}\varepsilon_{n}^{\,}=\left[\partial^{M}\widehat{H}\right]_{n,n}^{\,}+\sum_{\kappa=1}^{M-1}\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).$ (8)

  6. 6.

Having computed $\partial\varepsilon_{n}^{\,}$ through
$\partial^{N}\varepsilon_{n}^{\,}$, we assemble them into the Taylor
polynomial of degree $N$ at $\mathbf{k}_{0}^{\,}$:

$j^{N}(\varepsilon_{n}^{\,})|_{\mathbf{k}_{0}}=\varepsilon_{n}^{\,}+\frac{\partial\varepsilon_{n}^{\,}}{1!}+\frac{\partial^{2}\varepsilon_{n}^{\,}}{2!}+\cdots+\frac{\partial^{N}\varepsilon_{n}^{\,}}{N!}.$
(9)

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 Eq 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 Eq 4 as follows and apply
the rest of the algorithm

$\partial\widehat{H}\coloneqq\frac{\partial\widehat{H}(\mathbf{k}_{0}^{\,}+\lambda\,\mathbf{k},t_{0}^{\,}+\lambda\,\delta
t)}{\partial\lambda}\bigg{|}_{\lambda=0},\partial^{2}\widehat{H}\coloneqq\frac{\partial^{2}\widehat{H}(\mathbf{k}_{0}^{\,}+\lambda\,\mathbf{k},t_{0}^{\,}+\lambda\,\delta
t)}{\partial\lambda^{2}}\bigg{|}_{\lambda=0},\cdots\\\
\cdots,\partial^{N}\widehat{H}\coloneqq\frac{\partial^{N}\widehat{H}(\mathbf{k}_{0}^{\,}+\lambda\mathbf{k},t_{0}^{\,}+\lambda\,\delta
t)}{\partial\lambda^{N}}\bigg{|}_{\lambda=0}.$ (10)

This will then yield terms in the Taylor expansion that contain $\delta
t,\delta t^{2},\cdots$ 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 III.

###  II.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 Eqs 6 and 8 which
will be used later. The detailed derivations are included in the Appendix.

$\displaystyle\partial\varepsilon_{n}^{\,}=$ $\displaystyle\langle
n|\partial\widehat{H}|n\rangle,$ (11a)
$\displaystyle\partial^{2}\varepsilon_{n}^{\,}=$ $\displaystyle\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}}^{\,}},$
(11b) $\displaystyle\partial^{3}\varepsilon_{n}^{\,}=$ $\displaystyle\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}}^{\,}}+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}}^{\,}}$
$\displaystyle\quad+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}}^{\,})}$
$\displaystyle\quad-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}}.$
(11c)

Note that the formulae listed above only require us to be able to numerically
diagonalize the Hamiltonian $\widehat{H}(\mathbf{k})$ at a chosen
$\mathbf{k}_{0}^{\,}$, and be able to differentiate $\widehat{H}(\mathbf{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 $\mathbf{k}_{0}^{\,}$ and the matrix elements of
various derivatives of $\widehat{H}(\mathbf{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 (Eq 11a) is just the familiar Feynman
- Hellmann result.

##  III 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

$\mathbf{a}_{1}^{\phantom{{\dagger}}}=\begin{pmatrix}1\\\ \\\
0\end{pmatrix},\mathbf{a}_{2}^{\phantom{{\dagger}}}=\begin{pmatrix}-\frac{1}{2}\\\
\\\
\frac{\sqrt{3}}{2}\end{pmatrix},\mathbf{a}_{3}^{\phantom{{\dagger}}}=\begin{pmatrix}-\frac{1}{2}\\\
\\\ -\frac{\sqrt{3}}{2}\end{pmatrix}.$ (12)

Figure 1: The constant energy contours of the Haldane model near critical
tuning, in a Brillouin zone centred at the $K_{-}^{\,}$ point. In (a), the
staggered chemical potential is tuned slightly below the critical tuning
$-M_{0}^{\,}$. This causes the higher order critical point to split into a
number of ordinary critical points: a central minimum surrounded by three
regular saddle points. At the critical tuning (i.e when $\delta M=0$ or
equivalently $M=-M_{0}^{\,}$), we get a three-fold rotation symmetric higher
order singularity, known as the monkey saddle. As the staggered chemical
potential $M$ is tuned slightly above $-M_{0}^{\,}$ once again the higher
order critical point is split into a number of ordinary critical points: a
central maximum surrounded by three regular saddle points. This is depicted in
(c).

The next nearest neighbor vectors are given by
$\mathbf{b}_{1}^{\phantom{{\dagger}}}=\mathbf{a}_{2}^{\phantom{{\dagger}}}-\mathbf{a}_{3}^{\phantom{{\dagger}}}$,
$\mathbf{b}_{2}^{\phantom{{\dagger}}}=\mathbf{a}_{3}^{\phantom{{\dagger}}}-\mathbf{a}_{1}^{\phantom{{\dagger}}}$
and
$\mathbf{b}_{3}^{\phantom{{\dagger}}}=\mathbf{a}_{1}^{\phantom{{\dagger}}}-\mathbf{a}_{2}^{\phantom{{\dagger}}}$.
The full $\mathbf{k}$-space Hamiltonian is then the sum of nearest neighbor,
staggered chemical potential and next nearest neighbor terms:
$\widehat{H}(\mathbf{k})=\widehat{H}_{0}^{\phantom{{\dagger}}}(\mathbf{k})+M\sigma_{z}^{\phantom{{\dagger}}}+2t_{2}^{\phantom{{\dagger}}}\sum_{i}\sin(\mathbf{k}\cdot\mathbf{b}_{i}^{\phantom{{\dagger}}})\,\sigma_{z}^{\phantom{{\dagger}}}$,
where $\sigma_{z}^{\phantom{{\dagger}}}$ is the familiar Pauli matrix and

$\widehat{H}_{0}^{\phantom{{\dagger}}}(\mathbf{k})=\begin{pmatrix}0&t_{1}^{\phantom{{\dagger}}}\sum_{i}e^{i\mathbf{k}\cdot\mathbf{a}_{i}}\\\
t_{1}^{\phantom{{\dagger}}}\sum_{i}e^{-i\mathbf{k}\cdot\mathbf{a}_{i}}&0\end{pmatrix}.$
(13)

Since the Hamiltonian is two-dimensional, it can be diagonalized exactly to
give two bands indexed by $n=1,2$

$\varepsilon_{n}^{\phantom{{\dagger}}}(\mathbf{k})=(-1)^{n}\bigg{[}\bigg{(}M+2t_{2}\sum_{i}\sin(\mathbf{k}\cdot\mathbf{b}_{i}^{\phantom{{\dagger}}})\bigg{)}^{2}+t_{1}^{2}\bigg{|}\sum_{i}e^{i\mathbf{k}\cdot\mathbf{a}_{i}^{\phantom{{\dagger}}}}\bigg{|}^{2}\bigg{]}^{\frac{1}{2}}.$
(14)

Consider the following point of three-fold rotation symmetry

$\mathbf{K}_{-}^{\phantom{{\dagger}}}=\frac{4\pi}{3\sqrt{3}}\begin{pmatrix}0\\\
\\\ 1\end{pmatrix}.$ (15)

Let us define a critical value for the staggered chemical potential

$M_{0}^{\,}=\frac{t_{1}^{2}-18\,t_{2}^{2}}{2\,\sqrt{3}\,t_{2}^{\phantom{{\dagger}}}}.$
(16)

We can analytically derive the Taylor expansion for the upper band $(n=2)$
dispersion, which around the $K_{-}^{\phantom{{\dagger}}}$ point reads

$\displaystyle\varepsilon_{2}^{\phantom{{\dagger}}}(\mathbf{k}+\mathbf{K}_{-}^{\phantom{{\dagger}}})\approx$
$\displaystyle\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)$
$\displaystyle\quad+\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).$
(17)

When we tune $M\rightarrow-M_{0}^{\,}$

$\displaystyle\varepsilon_{2}^{\phantom{{\dagger}}}(\mathbf{k}+\mathbf{K}_{-}^{\phantom{{\dagger}}})\approx$
$\displaystyle\frac{t_{1}^{2}}{2\,\sqrt{3}\,|t_{2}^{\phantom{{\dagger}}}|}-\frac{3\,\sqrt{3}\,|t_{2}^{\phantom{{\dagger}}}|}{2}\left(k_{y}^{3}-3\,k_{y}^{\phantom{{\dagger}}}\,k_{x}^{2}\right)+\mathcal{O}\left(k^{4}\right).$
(18)

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 $\delta M$, the Taylor expansion gets
modified as

$\displaystyle\varepsilon_{2}^{\phantom{{\dagger}}}(\mathbf{k}+\mathbf{K}_{-}^{\phantom{{\dagger}}})\approx$
$\displaystyle\left|-\frac{t_{1}^{2}}{2\,\sqrt{3}\,t_{2}^{\,}}+\delta
M\right|+\frac{27\,t_{2}^{2}\,\,\delta
M}{2\,t_{1}^{2}}\left(k_{x}^{2}+k_{y}^{2}\right)$
$\displaystyle\quad-\left(\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{{\dagger}}}\,k_{x}^{2}\right)+\mathcal{O}\left(k^{4},\delta
M^{2}\right).$ (19)

As we see, a small quadratic term modulated by $\delta 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
Fig 1). Using the numerical values: $t=1$ and $t=1/2$:

$\displaystyle\varepsilon_{2}^{\phantom{{\dagger}}}(\mathbf{k}+\mathbf{K}_{-}^{\phantom{{\dagger}}})\approx$
$\displaystyle\,\frac{1}{\sqrt{3}}-\delta M+\frac{27}{8}\,\delta
M\left(k_{x}^{2}+k_{y}^{2}\right)$
$\displaystyle\quad-\left(\frac{3\sqrt{3}}{4}+\frac{27}{16}\,\delta
M\right)\left(k_{y}^{3}-3\,k_{y}^{\phantom{{\dagger}}}\,k_{x}^{2}\right)+\mathcal{O}\left(k^{4},\delta
M^{2}\right).$ (20)

We will now attempt to reproduce this Taylor expansion using the method
described above. Firstly we calculate the eigenvalues and eigenvectors at
$\mathbf{K}_{-}^{\,}$ for the _critical_ tuning (since $\delta M$ is treated
here as a perturbation, we are computing eigenvalues and eigenvectors of
$H(\mathbf{K}_{-}^{\,}+\lambda\,\mathbf{k};-M_{0}^{\,}+\lambda\,\delta M)$ at
$\lambda=0$, which is the same as computing them for the critical tuning). The
Hamiltonian takes the value

$\widehat{H}(\mathbf{K}_{-}^{\,})=\begin{pmatrix}-\frac{1}{\sqrt{3}}&0\\\
0&\frac{1}{\sqrt{3}}\end{pmatrix}$ (21)

Thus it is already diagonal. The eigenvalues and eigenvectors are then given
by

$\displaystyle\varepsilon_{1}^{\,}(\mathbf{K}_{-}^{\,})$
$\displaystyle=-\frac{1}{\sqrt{3}},\quad|E_{1}^{\,}\rangle=\begin{pmatrix}1\\\
0\end{pmatrix},$ $\displaystyle\varepsilon_{2}^{\,}(\mathbf{K}_{-}^{\,})$
$\displaystyle=\frac{1}{\sqrt{3}},\quad|E_{2}^{\,}\rangle=\begin{pmatrix}0\\\
1\end{pmatrix}.$ (22)

We now compute the matrices $\partial\widehat{H}$, $\partial^{2}\widehat{H}$
and $\partial^{3}\widehat{H}$, defined in Eq 4, in the eigenbasis:

$\displaystyle\partial\widehat{H}=$
$\displaystyle\left(\begin{array}[]{cc}\delta
M&\frac{3}{2}i\left(k_{x}+ik_{y}\right)\\\
-\frac{3}{2}i\left(k_{x}-ik_{y}\right)&-\delta M\\\ \end{array}\right)$ (23c)
$\displaystyle\partial^{2}\widehat{H}=$
$\displaystyle\left(\begin{array}[]{cc}\frac{9}{4}\sqrt{3}\left(k_{x}^{2}+k_{y}^{2}\right)&-\frac{3}{4}\left(k_{x}-ik_{y}\right){}^{2}\\\
-\frac{3}{4}\left(k_{x}+ik_{y}\right){}^{2}&-\frac{9}{4}\sqrt{3}\left(k_{x}^{2}+k_{y}^{2}\right)\\\
\end{array}\right)$ (23f) $\displaystyle\partial^{3}\widehat{H}=$
$\displaystyle\left(\begin{array}[]{cc}\frac{9}{8}\sqrt{3}k_{y}\left(k_{y}^{2}-3k_{x}^{2}\right)&-\frac{9}{8}i\left(k_{x}-ik_{y}\right)\left(k_{x}+ik_{y}\right){}^{2}\\\
\frac{9}{8}i\left(k_{x}-ik_{y}\right){}^{2}\left(k_{x}+ik_{y}\right)&-\frac{9}{8}\sqrt{3}k_{y}\left(k_{y}^{2}-3k_{x}^{2}\right)\\\
\end{array}\right)$ (23i)

The derivatives $\partial^{N}\varepsilon_{2}^{\,}$ for $N=1,2,3$ are then
given by

$\displaystyle\partial\varepsilon_{2}^{\,}=$ $\displaystyle\langle
E_{2}^{\,}|\partial\widehat{H}|E_{2}^{\,}\rangle$ $\displaystyle=$
$\displaystyle-\delta M,$ $\displaystyle\partial^{2}\varepsilon_{2}^{\,}=$
$\displaystyle\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}^{\,}(\mathbf{K}_{-}^{\,})-\varepsilon_{1}^{\,}(\mathbf{K}_{-}^{\,})}$
$\displaystyle=$
$\displaystyle-\frac{9\sqrt{3}}{4}\left(k_{x}^{2}+k_{y}^{2}\right)+2\,\frac{-\frac{3i}{2}\left(k_{x}^{\,}-ik_{y}^{\,}\right)\times\frac{3i}{2}\left(k_{x}^{\,}+ik_{y}^{\,}\right)}{\frac{1}{\sqrt{3}}-\left(-\frac{1}{\sqrt{3}}\right)}$
$\displaystyle=$ $\displaystyle 0,$
$\displaystyle\partial^{3}\varepsilon_{n}^{\,}=$ $\displaystyle\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}^{\,}(\mathbf{K}_{-}^{\,})-\varepsilon_{1}^{\,}(\mathbf{K}_{-}^{\,})}+3\,\frac{\langle
E_{2}^{\,}|\partial\widehat{H}|E_{1}^{\,}\rangle\langle
E_{1}^{\,}|\partial^{2}\widehat{H}|E_{2}^{\,}\rangle}{\varepsilon_{2}^{\,}(\mathbf{K}_{-}^{\,})-\varepsilon_{1}^{\,}(\mathbf{K}_{-}^{\,})}$
$\displaystyle\quad+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}^{\,}(\mathbf{K}_{-}^{\,})-\varepsilon_{1}^{\,}(\mathbf{K}_{-}^{\,})\right)\left(\varepsilon_{2}^{\,}(\mathbf{K}_{-}^{\,})-\varepsilon_{1}^{\,}(\mathbf{K}_{-}^{\,})\right)}$
$\displaystyle\quad-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}^{\,}(\mathbf{K}_{-}^{\,})-\varepsilon_{1}^{\,}(\mathbf{K}_{-}^{\,})\right)^{2}}$
$\displaystyle=$
$\displaystyle-\frac{9\sqrt{3}}{8}\left(k_{y}^{3}-3\,k_{x}^{2}k_{y}^{\,}\right)+3\,\frac{2\,\text{Re}\left[-\frac{3}{4}\left(k_{x}^{\,}+ik_{y}^{\,}\right)^{2}\times\frac{3i}{2}\left(k_{x}^{\,}+i\,k_{y}^{\,}\right)\right]}{\frac{1}{\sqrt{3}}-\left(-\frac{1}{\sqrt{3}}\right)}$
$\displaystyle\quad+6\,\frac{-\frac{3i}{2}\left(k_{x}^{\,}-ik_{y}^{\,}\right)\times\delta
M\times\frac{3i}{2}\left(k_{x}^{\,}+ik_{y}^{\,}\right)}{\left(\frac{1}{\sqrt{3}}-\left(-\frac{1}{\sqrt{3}}\right)\right)^{2}}$
$\displaystyle\quad-6(-\delta
M)\,\frac{-\frac{3i}{2}\left(k_{x}^{\,}-ik_{y}^{\,}\right)\times\frac{3i}{2}\left(k_{x}^{\,}+ik_{y}^{\,}\right)}{\left(\frac{1}{\sqrt{3}}-\left(-\frac{1}{\sqrt{3}}\right)\right)^{2}}$
$\displaystyle=$
$\displaystyle-\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).$ (24)

We can assemble these into the Taylor expansion using Eq 9:

$\varepsilon_{2}^{\phantom{{\dagger}}}(\mathbf{k}+\mathbf{K}_{-}^{\phantom{{\dagger}}})\approx\,\frac{1}{\sqrt{3}}-\delta
M+\frac{27}{8}\,\delta
M\left(k_{x}^{2}+k_{y}^{2}\right)-\frac{3\sqrt{3}}{4}\left(k_{y}^{3}-3\,k_{y}^{\phantom{{\dagger}}}\,k_{x}^{2}\right)+\cdots$
(25)

As it is evident, the Taylor expansion matches with the one we obtained by
directly Taylor expanding the band dispersion, except for the extra $\delta 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 $\lambda$ (basically from a term of the form
$\lambda^{4}\,\delta M\,\left(k_{y}^{3}-3k_{x}^{2}k_{y}\right)$), 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.

##  IV Conclusion

In this paper, we set in detail the method to find HOVHS in practice, when the
Hamiltonian $H(\mathbf{k})$ is known by various means. The knowledge of
$H(\mathbf{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 Sec II.2. These derivatives were then used to assemble the Taylor
expansion of the band to any specified order, at arbitrary $\mathbf{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 Sr2RuO4 or semiconducting ones like two dimensional InSe and
In2Se3 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].

Supporting Information

Supporting Information is available from the Wiley Online Library or from the
author.

Acknowledgements

We would like to thank Claudio Chamon, Mark Greenaway, Nick Hine, Phil King,
Samuel Magorrian, Edgar Abarca Morales, Andreas Rost and Peter Wahl for useful
discussions. The work was supported by the EPSRC grant EP/T034351/1.

##  Appendix A Derivation of the expressions

In this section, we explicitly prove the formula in Eq 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, i.e $\langle
n|n\rangle=1$.

###  A.1 First derivative

The eigenvalue equation reads

$H|n\rangle=\varepsilon_{n}^{\,}|n\rangle.$ (26)

Differentiating this once and applying the product rule for derivatives, we
get

$\partial H|n\rangle+H|\partial
n\rangle=\partial\varepsilon_{n}^{\,}|n\rangle+\varepsilon_{n}^{\,}|\partial
n\rangle.$ (27)

(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\rangle$

$\displaystyle\langle n|\partial H|n\rangle+\langle n|H|\partial n\rangle$
$\displaystyle=\partial\varepsilon_{n}^{\,}\langle
n|n\rangle+\varepsilon_{n}^{\,}\langle n|\partial n\rangle,$
$\displaystyle\implies\partial\varepsilon_{n}^{\,}$ $\displaystyle=\langle
n|\partial H|n\rangle.$ (28)

This is of course the familiar Feynman Hellmann result. We now take inner
product of Eq 27 with another eigenstate $|m\rangle$ that is not degenerate to
$|n\rangle$:

$\displaystyle\langle m|\partial H|n\rangle+\langle m|H|\partial n\rangle$
$\displaystyle=\partial\varepsilon_{n}^{\,}\langle
m|n\rangle+\varepsilon_{n}^{\,}\langle m|\partial n\rangle,$
$\displaystyle\implies\langle m|\partial n\rangle$
$\displaystyle=\frac{\langle m|\partial
H|n\rangle}{\varepsilon_{n}^{\,}-\varepsilon_{m}^{\,}}.$ (29)

Thus, we can find the overlap of the derivative of the eigenvector
$|n\rangle$, denoted by $|\partial n\rangle$ with the other eigenvectors _not
degenerate_ to it. It is not possible within this framework to compute
$\langle n|\partial n\rangle$ or even overlap with eigenvectors degenerate to
it. We can however differentiate $\langle n|n\rangle=1$ to obtain
$\text{Re}\,\langle n|\partial n\rangle=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 $\langle n|\partial n\rangle$ will not be an impediment to
computing the higher derivatives of the band dispersion
$\varepsilon_{n}^{\,}(\mathbf{k})$, as we shall show below. We also compute
the overlap of Eq 27 with respect to another eigenstate that is degenerate to
$|n\rangle$, denoted by $|n_{1}^{\,}\rangle\in D(n)$ (where $D(n)$ is the set
of eigenvectors degenerate to $n$)

$\displaystyle\langle n_{1}^{\,}|\partial H|n\rangle+\langle
n_{1}^{\,}|H|\partial n\rangle$
$\displaystyle=\partial\varepsilon_{n}^{\,}\langle
n_{1}^{\,}|n\rangle+\varepsilon_{n}^{\,}\langle n_{1}^{\,}|\partial n\rangle,$
$\displaystyle\implies\langle n_{1}^{\,}|\partial H|n\rangle$
$\displaystyle=0,\,\,\text{for}\,\,\varepsilon_{n_{1}}^{\,}=\varepsilon_{n}^{\,}\,\,\text{but}\,|n_{1}^{\,}\rangle\neq|n\rangle.$
(30)

In practice however, we may find the projection of the first derivative of
Hamiltonian matrix $\partial H$ to the degenerate subspace to be non-diagonal.
That is, we may find $\langle n_{1}^{\,}|\partial H|n_{2}^{\,}\rangle\neq 0$,
for distinct but degenerate $|n_{1}^{\,}\rangle$ and $|n_{2}^{\,}\rangle$. In
such situations, we diagonalise the projection of $\partial H$ in this
subspace, and replace the original set of $|n_{i}^{\,}\rangle$ with the new
eigenbasis, so that Eq 30 is automatically satisfied.

###  A.2 Second derivative

Differentiating Eq 27 once again we get

$\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.$ (31)

Taking inner product with $|n\rangle$, we get

$\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+2\,\partial\varepsilon_{n}^{\,}\,\langle n|\partial
n\rangle+\varepsilon_{n}^{\,}\langle n|\partial^{2}n\rangle$ (32)

On rearranging terms, inserting a resolution of identity and splitting into
projectors over the degenerate and non-degenerate subspaces, this becomes

$\displaystyle\partial^{2}\varepsilon_{n}^{\,}=$ $\displaystyle\langle
n|\partial^{2}H|n\rangle+2\,\langle n|\partial H|\partial
n\rangle-2\,\partial\varepsilon_{n}^{\,}\,\langle n|\partial n\rangle$
$\displaystyle=$ $\displaystyle\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+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$
$\displaystyle=$ $\displaystyle\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}^{\,}},$ (33)

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\neq|n\rangle$ in going to the
last line. We now compute the overlap of the terms in Eq 31 with
$|m\rangle\notin D(n)$:

$\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+2\,\partial\varepsilon_{n}^{\,}\langle m_{1}^{\,}|\partial
n\rangle+\varepsilon_{n}^{\,}\langle m_{1}^{\,}|\partial^{2}n\rangle.$ (34)

Simplifying this further by rearranging terms and introducing another
resolution of identity:

$\displaystyle\langle m_{1}^{\,}|\partial^{2}n\rangle=$
$\displaystyle\frac{1}{\varepsilon_{n}^{\,}-\varepsilon_{m_{1}}^{\,}}\big{[}\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\big{]}$ $\displaystyle=$
$\displaystyle\frac{1}{\varepsilon_{n}^{\,}-\varepsilon_{m_{1}}^{\,}}\bigg{[}\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$
$\displaystyle\quad+2\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\bigg{]}$ $\displaystyle=$ $\displaystyle\frac{\langle
m_{1}^{\,}|\partial^{2}H|n\rangle}{\varepsilon_{n}^{\,}-\varepsilon_{m_{1}}^{\,}}+2\sum_{|m_{2}\rangle\notin
D(n)}\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}}^{\,})}-2\,\partial\varepsilon_{n}^{\,}\frac{\langle
m_{1}^{\,}|\partial
H|n\rangle}{(\varepsilon_{n}^{\,}-\varepsilon_{m_{1}}^{\,})^{2}}$
$\displaystyle\quad+2\sum_{n_{1}^{\,}\in D(n)}\frac{\langle
m_{1}^{\,}|\partial
H|n_{1}^{\,}\rangle}{\varepsilon_{n}^{\,}-\varepsilon_{m_{1}}^{\,}}\langle
n_{1}^{\,}|\partial n\rangle.$ (35)

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\rangle$. Lastly, we compute the overlap of Eq 33 with
$|n_{1}^{\,}\rangle\in D(n)$ and $|n_{1}^{\,}\rangle\neq|n\rangle$

$\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+2\,\partial\varepsilon_{n}^{\,}\langle n_{1}^{\,}|\partial
n\rangle+\varepsilon_{n}^{\,}\langle n_{1}^{\,}|\partial^{2}n\rangle.$ (36)

Using $\langle n_{1}^{\,}|H=\varepsilon_{n}^{\,}\langle n_{1}^{\,}|$,
introducing resolution of identity and simplifying

$\displaystyle\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+2\sum_{|n_{2}\rangle\notin D(n)}\langle n_{1}^{\,}|\partial
H|n_{2}^{\,}\rangle\langle n_{2}^{\,}|\partial n\rangle=$ $\displaystyle
2\,\partial\varepsilon_{n}^{\,}\langle n_{1}^{\,}|\partial n\rangle,$
$\displaystyle\implies\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}}^{\,}}=$ $\displaystyle
0.$ (37)

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\rangle$, 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

$\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+3\,\partial^{2}\varepsilon_{n}^{\,}|\partial
n\rangle+3\,\partial\varepsilon_{n}^{\,}|\partial^{2}n\rangle+\varepsilon_{n}^{\,}|\partial^{3}n\rangle.$
(38)

We compute the overlap of this equation with $|n\rangle$, rearrange terms,
introduce resolutions of identity (splitting them into projectors over the
degenerate and non-degenerate subspaces) and do some simplifications

$\displaystyle\partial^{3}\varepsilon_{n}^{\,}=$ $\displaystyle\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-3\,\partial^{2}\varepsilon_{n}^{\,}\,\langle
n|\partial n\rangle-3\,\partial\varepsilon_{n}^{\,}\,\langle
n|\partial^{2}n\rangle$ $\displaystyle=$ $\displaystyle\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+3\sum_{|m_{1}\rangle\notin D(n)}\langle n|\partial
H|m_{1}^{\,}\rangle\langle m_{1}^{\,}|\partial^{2}n\rangle$
$\displaystyle\quad+3\sum_{|n_{1}^{\,}\rangle\in D(n)}\langle
n|\partial^{2}H|n_{1}^{\,}\rangle\langle n_{1}^{\,}|\partial
n\rangle+3\sum_{|n_{1}^{\,}\rangle\in D(n)}\langle n|\partial
H|n_{1}^{\,}\rangle\langle n_{1}^{\,}|\partial^{2}n\rangle$
$\displaystyle\quad-3\,\,\partial^{2}\varepsilon_{n}^{\,}\,\langle n|\partial
n\rangle-3\,\,\partial\varepsilon_{n}^{\,}\,\langle n|\partial^{2}n\rangle.$
(39)

We can substitute for $\langle m_{1}^{\,}|\partial n\rangle$ and $\langle
m_{1}^{\,}|\partial^{2}n\rangle$ from Eqs 29 and 35 respectively to obtain

$\displaystyle\partial^{3}\varepsilon_{n}^{\,}=$ $\displaystyle\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}}^{\,}}$
$\displaystyle\quad+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}}^{\,}}$
$\displaystyle\quad+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}}^{\,})}$
$\displaystyle\quad-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}}$
$\displaystyle\quad+(*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$
$\displaystyle\quad+(*2)\,\,3\sum_{|n_{1}^{\,}\rangle\in D(n)}\langle
n|\partial^{2}H|n_{1}^{\,}\rangle\langle n_{1}^{\,}|\partial n\rangle$
$\displaystyle\quad+(*3)\,\,3\sum_{|n_{1}^{\,}\rangle\in D(n)}\langle
n|\partial H|n_{1}^{\,}\rangle\langle n_{1}^{\,}|\partial^{2}n\rangle$
$\displaystyle\quad-3\,\,\partial^{2}\varepsilon_{n}^{\,}\,\langle n|\partial
n\rangle-3\,\,\partial\varepsilon_{n}^{\,}\,\langle n|\partial^{2}n\rangle.$
(40)

Now let us look at some of the terms marked with a *. To simplify them, we
pull out the term containing $|n\rangle$ from the sum over
$|n_{1}^{\,}\rangle\in D(n)$ and substitute the expression for
$\partial^{2}\varepsilon_{n}^{\,}$:

$\displaystyle(*2)+($ $\displaystyle*1)$ $\displaystyle=$
$\displaystyle\,3\sum_{|n_{1}^{\,}\rangle\in D(n)}\langle
n|\partial^{2}H|n_{1}^{\,}\rangle\langle n_{1}^{\,}|\partial
n\rangle+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$ $\displaystyle=$
$\displaystyle\,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$
$\displaystyle+3\sum_{\begin{subarray}{c}|n_{1}^{\,}\rangle\in D(n)\\\
|n_{1}^{\,}\rangle\neq|n\rangle\end{subarray}}\bigg{[}\langle
n|\partial^{2}H|n_{1}^{\,}\rangle+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}}^{\,}}\bigg{]}\langle
n_{1}^{\,}|\partial n_{1}^{\,}\rangle$ $\displaystyle=$
$\displaystyle\,3\,\,\partial^{2}\varepsilon_{n}^{\,}\,\langle n|\partial
n\rangle,$ (41)

where we substitute from Eq 33 for the terms within parenthesis, while the
terms within [ ] vanish due to Eq 37 (since
$|n_{1}^{\,}\rangle\neq|n\rangle$). Lastly, since $\langle n|\partial
H|n\rangle=\partial\varepsilon_{n}^{\,}$ (from Eq 28) and $\langle n|\partial
H|n_{1}^{\,}\rangle=0$ for $|n\rangle\neq|n_{1}^{\,}\rangle$ (from Eq 30), we
have

$(*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.$ (42)

Therefore the sum $(*1)+(*2)+(*3)$ cancels exactly with the last line of Eq 40
so that the expression for the third derivative of the band becomes

$\displaystyle\partial^{3}\varepsilon_{n}^{\,}=$ $\displaystyle\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}}^{\,}}+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}}^{\,}}$
$\displaystyle\quad+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}}^{\,})}$
$\displaystyle\quad-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}}.$ (43)

To obtain this, we got rid of the terms containing $\langle n|\partial
n\rangle$ and $\langle n|\partial^{2}n\rangle$ 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\neq|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:

$\displaystyle\partial^{N}(H\,|n\rangle)$
$\displaystyle=\partial^{N}(\varepsilon_{n}^{\,}\,|n\rangle)$ (44a)
$\displaystyle\implies\sum_{K=0}^{N}\begin{pmatrix}N\\\
K\end{pmatrix}\partial^{N-K}H\,|\partial^{K}n\rangle$
$\displaystyle=\sum_{K=0}^{N}\begin{pmatrix}N\\\
K\end{pmatrix}\partial^{N-K}\varepsilon_{n}^{\,}\,|\partial^{K}n\rangle$ (44b)

Isolating the $K=0$ and $K=N$ terms on both sides of the equation, we obtain

$H\,|\partial^{N}n\rangle+\partial^{N}H\,|n\rangle+\sum_{K=1}^{N-1}\begin{pmatrix}N\\\
K\end{pmatrix}\partial^{N-K}H\,|\partial^{K}n\rangle=\varepsilon_{n}^{\,}\,|\partial^{N}n\rangle+\partial^{N}\varepsilon_{n}^{\,}\,|n\rangle+\sum_{K=1}^{N-1}\begin{pmatrix}N\\\
K\end{pmatrix}\partial^{N-K}\varepsilon_{n}^{\,}\,|\partial^{K}n\rangle$ (45)

After taking inner product on both sides with $|m\rangle\notin D(n)$ (i.e an
eigenstate not degenerate to $|n\rangle$), using $\langle
m|H=\varepsilon_{m}^{\,}\,\langle m|$ and rearranging, we obtain:

$\langle
m|\partial^{N}n\rangle=\frac{1}{\varepsilon_{n}^{\,}-\varepsilon_{m}^{\,}}\left[\langle
m|\partial^{N}H|n\rangle+\sum_{K=1}^{N-1}\begin{pmatrix}N\\\
K\end{pmatrix}\left(\langle
m|\partial^{N-K}H**|\partial^{K}n\rangle-\partial^{N-K}\varepsilon_{n}^{\,}\,\langle
m|\partial^{K}n\rangle\right)\right]$ (46)

We now introduce a resolution of identity in terms of the eigenvectors of the
Hamiltonian at the given $\mathbf{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\rightarrow
K_{i}$ so that the expression becomes more readily usable in the forthcoming
calculations:

$\langle
m_{i}^{\,}|\partial^{K_{i}}n\rangle=\frac{1}{\varepsilon_{n}^{\,}-\varepsilon_{m_{i}}^{\,}}\bigg{[}\langle
m_{i}^{\,}|\partial^{K_{i}}H|n\rangle+\sum_{K_{i+1}=1}^{K_{i}-1}\begin{pmatrix}K_{i}\\\
K_{i+1}\end{pmatrix}\sum_{|m_{i+1}\rangle\notin D(n)}\bigg{\\{}\\\
\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\bigg{\\}}\\\
+\sum_{K_{i+1}=1}^{K_{i}-1}\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\bigg{]}$
(47)

To make the notation less cumbersome, for a given $|n\rangle$, we define

$\displaystyle\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}}$
(48a) $\displaystyle\mathcal{H}_{m_{i},m_{j}}^{K}\coloneqq\langle
m_{i}^{\,}|\partial^{K}H|m_{j}^{\,}\rangle$ (48b)
$\displaystyle\mathcal{O}_{m_{i},m_{j}}^{K}\coloneqq\langle
m_{i}^{\,}|\partial^{K}m_{j}^{\,}\rangle\text{, for non-degenerate
}|m_{i}^{\,}\rangle,|m_{j}^{\,}\rangle$ (48c)
$\displaystyle\mathcal{R}_{\tilde{n},n}^{K}\coloneqq\langle\tilde{n}|\partial^{K}n\rangle\text{,
for degenerate }|\tilde{n}\rangle,|n\rangle.$ (48d)

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 Eq 47 in a compact form:

$\mathcal{O}_{m_{i},n}^{K_{i}}=\frac{1}{\varepsilon_{n}^{\,}-\varepsilon_{m_{i}}^{\,}}\bigg{[}\mathcal{H}_{m_{i},n}^{K_{i}}+\sum_{K_{i+1}=1}^{K_{i}-1}\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}}\\\
+\sum_{K_{i+1}=1}^{K_{i}-1}\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}}\bigg{]}$
(49)

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}}$:

$\displaystyle\mathcal{O}_{m_{i},n}^{K_{i}}=$
$\displaystyle\frac{\mathcal{H}_{m_{i},n}^{K_{i}}}{\varepsilon_{n}^{\,}-\varepsilon_{m_{i}}^{\,}}+\sum_{K_{i+1}=1}^{K_{i}-1}\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}}^{\,}}\bigg{[}\frac{\mathcal{H}_{m_{i+1},n}^{K_{i+1}}}{\varepsilon_{n}^{\,}-\varepsilon_{m_{i+1}}^{\,}}$
$\displaystyle\qquad+\sum_{K_{i+2}=1}^{K_{i+1}-1}\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}}$
$\displaystyle\qquad+\sum_{K_{i+2}=1}^{K_{i+1}-1}\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}}\bigg{]}$
$\displaystyle\qquad+(*)\sum_{K_{i+1}=1}^{K_{i}-1}\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}}.$
(50)

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 $(*)$).

$\displaystyle\mathcal{O}_{m_{i},n}^{K_{i}}=$
$\displaystyle\frac{\mathcal{H}_{m_{i},n}^{K_{i}}}{\varepsilon_{n}^{\,}-\varepsilon_{m_{i}}^{\,}}+\sum_{K_{i+1}=1}^{K_{i}-1}\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}}^{\,}}$
$\displaystyle+\sum_{K_{i+1}=1}^{K_{i}-1}\begin{pmatrix}K_{i}\\\
K_{i+1}\end{pmatrix}\sum_{K_{i+2}=1}^{K_{i+1}-1}\begin{pmatrix}K_{i+1}\\\
K_{i+2}\end{pmatrix}\sum_{|m_{i+1}\rangle\notin
D(n)}\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}}$
$\displaystyle+(*)\sum_{K_{i+1}=1}^{K_{i}-1}\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}}$
$\displaystyle+\sum_{K_{i+1}=1}^{K_{i}-1}\begin{pmatrix}K_{i}\\\
K_{i+1}\end{pmatrix}\sum_{K_{i+2}=1}^{K_{i+1}-1}\begin{pmatrix}K_{i+1}\\\
K_{i+2}\end{pmatrix}\sum_{|m_{i+1}\rangle\notin
D(n)}\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}}.$
(51)

We can expand this further by substituting the expression for
$\mathcal{O}_{m_{i+2},n}^{K_{i+2}}$ based on Eq 49. That can be further
expanded by substituting the expression for
$\mathcal{O}_{m_{i+3},n}^{K_{i+3}}$, once again based on Eq 49. We cannot
however continue indefinitely expanding like this. The reason is that we have
$1\leqslant K_{j+1}\leqslant K_{j}-1$, which, on repeated application,
starting from $K_{i+l}$ up to $K_{i}$, gives $1\leqslant K_{i+l}\leqslant
K_{i+l-1}-1\leqslant K_{i+l-2}-2\leqslant\cdots\leqslant K_{i}-l$. This means
that the maximum allowed $l$ itself equals $K_{i}-1$ (since $1\leqslant
K_{i+l}\leqslant K_{i}-l$). This also gives the upper limit of the recursive
expansion. The full expanded form then reads

$\displaystyle\mathcal{O}_{m_{i},n}^{K_{i}}=$
$\displaystyle\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)}\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}}^{\,}}$
$\displaystyle+\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)}\begin{pmatrix}K_{i}\\\
K_{i+1}\end{pmatrix}\begin{pmatrix}K_{i+1}\\\
K_{i+2}\end{pmatrix}\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}}^{\,}}$
$\displaystyle+\cdots$ $\displaystyle+\cdots$
$\displaystyle+\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\neq|n\rangle}\sum_{|m_{i+2}\rangle\neq|n\rangle}\cdots\sum_{|m_{i+l}\rangle\neq|n\rangle}\bigg{[}\begin{pmatrix}K_{i}\\\
K_{i+1}\end{pmatrix}\begin{pmatrix}K_{i+1}\\\
K_{i+2}\end{pmatrix}\cdots\begin{pmatrix}K_{i+l-1}\\\ K_{i+l}\end{pmatrix}$
$\displaystyle\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}}^{\,}}\bigg{]}$
$\displaystyle+(*)\sum_{K_{i+1}=1}^{K_{i}-1}\sum_{|\tilde{n}\rangle\in
D(n)}\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}}$
$\displaystyle+\sum_{K_{i+1}=1}^{K_{i}-1}\sum_{K_{i+2}=1}^{K_{i+1}-1}\sum_{|m_{i+1}\rangle\neq|n\rangle}\sum_{|\tilde{n}\rangle\in
D(n)}\begin{pmatrix}K_{i}\\\ K_{i+1}\end{pmatrix}\begin{pmatrix}K_{i+1}\\\
K_{i+2}\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},\tilde{n}}^{K_{i+1}-K_{i+2}}}{\varepsilon_{n}^{\,}-\varepsilon_{m_{i+1}}^{\,}}\,\mathcal{R}_{\tilde{n},n}^{K_{i+2}}$
$\displaystyle+\cdots$ $\displaystyle+\cdots$
$\displaystyle+\sum_{K_{i+1}=1}^{K_{i}-1}\cdots\sum_{K_{i+l}=1}^{K_{i+l-1}-1}\sum_{|m_{i+1}\rangle\neq|n\rangle}\cdots\sum_{|m_{i+l-1}\rangle\neq|n\rangle}\sum_{|\tilde{n}\rangle\in
D(n)}\bigg{[}\begin{pmatrix}K_{i}\\\
K_{i+1}\end{pmatrix}\cdots\begin{pmatrix}K_{i+l-1}\\\ K_{i+l}\end{pmatrix}$
$\displaystyle\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-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}}^{\,}}\bigg{]}\mathcal{R}_{\tilde{n},n}^{K_{i+l}},$
(52)

where the sequence of terms ends with maximum $l$ equalling $K_{i}-1$.

We now proceed to compute the $N^{\text{th}}$ derivative of band $n$ by taking
inner product of Eq 45 with $|n\rangle$, rearranging terms, introducing
resolution of identity and separating out the projector to the degenerate
subspace to obtain

$\displaystyle\partial^{N}\varepsilon_{n}^{\,}=$ $\displaystyle\langle
n|\partial^{N}H|n\rangle+\sum_{K_{1}=1}^{N-1}\begin{pmatrix}N\\\
K_{1}\end{pmatrix}\langle
n|\partial^{N-K_{1}}H|\partial^{K_{1}}n\rangle-\sum_{K_{1}=1}^{N-1}\begin{pmatrix}N\\\
K_{1}\end{pmatrix}\partial^{N-K_{1}}\varepsilon_{n}^{\,}\langle
n|\partial^{K_{1}}n\rangle$ (53a) $\displaystyle=$ $\displaystyle\langle
n|\partial^{N}H|n\rangle+\sum_{K_{1}=1}^{N-1}\sum_{|m_{1}\rangle\notin
D(n)}\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$
$\displaystyle-\sum_{K_{1}=1}^{N-1}\begin{pmatrix}N\\\
K_{1}\end{pmatrix}\partial^{N-K_{1}}\varepsilon_{n}^{\,}\langle
n|\partial^{K_{1}}n\rangle+\sum_{K_{1}=1}^{N-1}\sum_{|\tilde{n}\rangle\in
D(n)}\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.$
(53b)

In the symbolic notation, the above reads

$\displaystyle\partial^{N}\varepsilon_{n}^{\,}=$
$\displaystyle\mathcal{H}_{n,n}^{N}+\sum_{K_{1}=1}^{N-1}\sum_{|m_{1}\rangle\notin
D(n)}\begin{pmatrix}N\\\
K_{1}\end{pmatrix}\mathcal{H}_{n,m_{1}}^{N-K_{1}}\,\mathcal{O}_{m_{1},n}^{K_{1}}-\sum_{K_{1}=1}^{N-1}\begin{pmatrix}N\\\
K_{1}\end{pmatrix}\partial^{N-K_{1}}\varepsilon_{n}^{\,}\,\mathcal{R}^{K_{1}}_{n,n}$
$\displaystyle+\sum_{K_{1}=1}^{N-1}\sum_{\tilde{n}\notin
D(n)}\begin{pmatrix}N\\\
K_{1}\end{pmatrix}\mathcal{H}^{N-K_{1}}_{n,\tilde{n}}\,\mathcal{R}^{K_{1}}_{\tilde{n},n}$
(54)

We now substitute for $O_{m_{1}}^{K_{1}}$:

$\displaystyle\partial^{N}\varepsilon_{n}^{\,}=$
$\displaystyle\mathcal{H}_{n,n}^{N}+\sum_{K_{1}=1}^{N-1}\sum_{|m_{1}\rangle\notin
D(n)}\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}}^{\,}}$
$\displaystyle+\sum_{K_{1}=1}^{N-1}\sum_{K_{2}=1}^{K_{1}-1}\sum_{|m_{1}\rangle\neq|n\rangle}\sum_{|m_{2}\rangle\neq|n\rangle}\begin{pmatrix}N\\\
K_{1}\end{pmatrix}\begin{pmatrix}K_{1}\\\
K_{2}\end{pmatrix}\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}}^{\,}}$
$\displaystyle+\cdots$ $\displaystyle+\cdots$
$\displaystyle+\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\neq|n\rangle}\sum_{|m_{2}\rangle\neq|n\rangle}\cdots\sum_{|m_{N-1}\rangle\neq|n\rangle}\begin{pmatrix}N\\\
K_{1}\end{pmatrix}\begin{pmatrix}K_{1}\\\
K_{2}\end{pmatrix}\cdots\begin{pmatrix}K_{N-2}\\\
K_{N-1}\end{pmatrix}\bigg{[}$
$\displaystyle\qquad\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\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{]}$
$\displaystyle-(*1)\sum_{K_{1}=1}^{N-1}\begin{pmatrix}N\\\
K_{1}\end{pmatrix}\partial^{N-K_{1}}\varepsilon_{n}^{\,}\,\mathcal{R}^{K_{1}}_{n,n}$
$\displaystyle+(*2)\sum_{K_{1}=1}^{N-1}\sum_{|\tilde{n}\rangle\in
D(n)}\begin{pmatrix}N\\\
K_{1}\end{pmatrix}\mathcal{H}^{N-K_{1}}_{n,\tilde{n}}\,\mathcal{R}^{K_{1}}_{\tilde{n},n}$
$\displaystyle+\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)}\begin{pmatrix}N\\\
K_{1}\end{pmatrix}\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}$
$\displaystyle+\cdots$ $\displaystyle+\cdots$
$\displaystyle+\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)}\begin{pmatrix}N\\\ K_{1}\end{pmatrix}\begin{pmatrix}K_{1}\\\
K_{2}\end{pmatrix}\begin{pmatrix}K_{2}\\\
K_{3}\end{pmatrix}\cdots\begin{pmatrix}K_{N-2}\\\
K_{N-1}\end{pmatrix}\bigg{[}$
$\displaystyle\qquad\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\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}.$
(55)

Theorem: The term $(*1)$ cancels with terms from $(2*)$ onwards so that the
band derivative takes the explicitly calculable form:

$\displaystyle\partial^{N}\varepsilon_{n}^{\,}=$
$\displaystyle\mathcal{H}_{n,n}^{N}+\sum_{K_{1}=1}^{N-1}\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}}^{\,}}$
$\displaystyle+\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)}\begin{pmatrix}N\\\
K_{1}\end{pmatrix}\begin{pmatrix}K_{1}\\\
K_{2}\end{pmatrix}\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}}^{\,}}$
$\displaystyle+\cdots$ $\displaystyle+\cdots$
$\displaystyle+\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)}\begin{pmatrix}N\\\ K_{1}\end{pmatrix}\begin{pmatrix}K_{1}\\\
K_{2}\end{pmatrix}\begin{pmatrix}K_{2}\\\
K_{3}\end{pmatrix}\cdots\begin{pmatrix}K_{N-2}\\\
K_{N-1}\end{pmatrix}\bigg{[}$
$\displaystyle\qquad\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\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{]}.$
(56)

Furthermore, for any $|\tilde{n}\rangle$ distinct from, but degenerate to
$|n\rangle$, the following is satisfied

$\displaystyle 0=$
$\displaystyle\mathcal{H}_{\tilde{n},n}^{N}+\sum_{K_{1}=1}^{N-1}\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}}^{\,}}$
$\displaystyle+\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)}\begin{pmatrix}N\\\
K_{1}\end{pmatrix}\begin{pmatrix}K_{1}\\\
K_{2}\end{pmatrix}\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}}^{\,}}$
$\displaystyle+\cdots$ $\displaystyle+\cdots$
$\displaystyle+\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)}\begin{pmatrix}N\\\ K_{1}\end{pmatrix}\begin{pmatrix}K_{1}\\\
K_{2}\end{pmatrix}\begin{pmatrix}K_{2}\\\
K_{3}\end{pmatrix}\cdots\begin{pmatrix}K_{N-2}\\\
K_{N-1}\end{pmatrix}\bigg{[}$
$\displaystyle\qquad\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\cdots\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{]}.$
(57)

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
Eq 11. As for Eq 57, 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 $\langle n|$. As Eqs 30 and 37 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\leqslant M\leqslant
N-1$. To prove that they are satisfied for $M=N$, we have to show that the
terms from $(2*)$ cancel with $(1*)$ in Eq 55. Now $(2*)$ is a sum of several
individual terms which take the form

$\begin{pmatrix}N\\\ K_{1}\end{pmatrix}\begin{pmatrix}K_{1}\\\
K_{2}\end{pmatrix}\begin{pmatrix}K_{2}\\\
K_{3}\end{pmatrix}\cdots\begin{pmatrix}K_{l-1}\\\ K_{l}\end{pmatrix}\\\
\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\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},$
(58)

for _all possible_ ordered sequences of $l$ natural numbers
$\\{K_{1},K_{2},\cdots,K_{l}\\}$ such that
$N>K_{1}>K_{2}>\cdots>K_{l}\geqslant 1$ and $1\leqslant l\leqslant 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 Eq 55

$\sum_{K_{1}=1}^{N-1}\sum_{K_{2}=1}^{K_{1}-1}\cdots\sum_{K_{l}=1}^{K_{l-1}-1}$
(59)

Let us denote by $A_{l}$, the set of all such sequences of length $l$ (with
$1\leqslant l\leqslant N-1$). That is, $A_{l}$ = Set of all ordered sequences
$\\{K_{1},K_{2},\cdots,K_{l}\\}$ such that $K_{i}\in\mathbb{N}$ and
$N>K_{1}>K_{2}>\cdots>K_{l}\geqslant 1$. It is then easy to see that

$\sum_{K_{1}=1}^{N-1}\sum_{K_{2}=1}^{K_{1}-1}\cdots\sum_{K_{l}=1}^{K_{l-1}-1}\equiv\sum_{\\{K_{1},\cdots,K_{l}\\}\in
A_{l}}.$ (60)

Let us define another set $B_{l}$ which is the set of all ordered sequences of
$l$ natural numbers
$\\{K_{l},\widetilde{K}_{1},\widetilde{K}_{2},\cdots\widetilde{K}_{l-1}\\}$
such that $N-l\geqslant K_{l}\geqslant 1$ and
$N-K_{l}>\widetilde{K}_{1}>\widetilde{K}_{2}>\cdots>\widetilde{K}_{l-1}\geqslant
1$. It is also easy to see that

$\sum_{K_{l}=1}^{N-l}\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}\equiv\sum_{\\{K_{l},\widetilde{K}_{1},\cdots\widetilde{K}_{l-1}\\}\in
B_{l}}.$ (61)

Lemma: There exists a bijection between $A_{l}$ and $B_{l}$ for _each_ $l$.
Proof: Let us define $F:A_{l}\rightarrow B_{l}$, with
$F(K_{1},K_{2},\cdots,K_{l-1},K_{l})=\\{K_{l},K_{1}-K_{l},K_{2}-K_{l},\cdots,K_{l-1}-K_{l}\\}$.
Since $N>K_{1}>K_{2}>\cdots>K_{l}\geqslant 1$, it is easy to see that the
range is indeed contained in $B_{l}$ since $N-l\geqslant K_{l}\geqslant 1$ and
$N-K_{l}>K_{1}-K_{l}>K_{2}-K_{l}>\cdots>K_{l-1}-K_{l}>0$. (The last inequality
$K_{l-1}-K_{l}>0$ is equivalent to $K_{l-1}-K_{l}\geqslant 1$ since we are
dealing with natural numbers).

Let us define the function $G:B_{l}\rightarrow A_{l}$ with
$G(K_{l},\widetilde{K}_{1},\widetilde{K}_{2},\cdots\widetilde{K}_{l-1})=\\{\widetilde{K}_{1}+K_{l},\widetilde{K}_{2}+K_{l},\cdots,\widetilde{K}_{l-1}+K_{l},K_{l}\\}$.
Since $N-l\geqslant K_{l}\geqslant 1$ and
$N-K_{l}>\widetilde{K}_{1}>\widetilde{K}_{2}>\cdots>\widetilde{K}_{l-1}\geqslant
1$, it is easy to see that
$N>\widetilde{K}_{1}+K_{l}>\widetilde{K}_{2}+K_{l}>\cdots>\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}$:

$\displaystyle\sum_{l=1}^{N-1}\sum_{\\{K_{1},\cdots,K_{l}\\}\in
A_{l}}\psi_{l}^{\,}(K_{1},\cdots,K_{l})\equiv\sum_{l=1}^{N-1}\sum_{\\{K_{l},\widetilde{K}_{1},\cdots,\widetilde{K}_{l-1}\\}\in
B_{l}}\psi_{l}^{\,}(G(K_{l},\widetilde{K}_{1},\cdots,\widetilde{K}_{l-1}))$
$\displaystyle\implies$
$\displaystyle\sum_{l=1}^{N-1}\sum_{K_{1}=1}^{N-1}\sum_{K_{2}=1}^{K_{1}-1}\cdots\sum_{K_{l}=1}^{K_{l-1}-1}\psi_{l}^{\,}(K_{1},K_{2},\cdots,K_{l})$
$\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}\cdots\sum_{\widetilde{K}_{l-1}=1}^{\widetilde{K}_{l-2}-1}\psi_{l}^{\,}(G(K_{l},\widetilde{K}_{1},\widetilde{K}_{2},\cdots,\widetilde{K}_{l-1})).$
(62)

The sums on LHS and RHS exhaustively go through all the elements in $A_{l}$
and $B_{l}$ respectively, as seen from Eqs 60 and 61. 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.

Figure 2: Two equivalent ways of iterating over a grid of points. The left
panel represents the sum $\sum_{l=1}^{N-1}\sum_{K_{l}=1}^{N-l}$, while the
right panel denotes
$\sum_{\mathfrak{K}_{1}=1}^{N-1}\sum_{\mathfrak{l}=1}^{N-\mathfrak{K}_{1}}$.

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\leqslant
l\leqslant N-1$ and $1\leqslant K_{l}\leqslant N-l$ and $D$ to be the set of
all ordered pairs $\\{\mathfrak{K}_{1},\mathfrak{l}\\}$ 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})=\\{K_{l},l\\}$ and $\mathcal{G}:D\rightarrow C$
with
$\mathcal{G}(\mathfrak{K}_{1},\mathfrak{l})=\\{\mathfrak{l},\mathfrak{K}_{1}\\}$.
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:

$\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}}.$
(63)

This is geometrically depicted in Fig 2, where the grid points in the left
panel are being summed over exhaustively by the sum on LHS in Eq 63. 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

$\displaystyle\sum_{l=1}^{N-1}\sum_{K_{1}=1}^{N-1}\sum_{K_{2}=1}^{K_{1}-1}\cdots\sum_{K_{l}=1}^{K_{l-1}-1}\psi_{l}^{\,}(K_{1},K_{2},\cdots,K_{l})$
$\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}\cdots\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},\cdots,\widetilde{K}_{\mathfrak{l}-1}))$
$\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}\cdots\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},\cdots,\widetilde{K}_{\mathfrak{l}-1}+\mathfrak{K}_{1},\mathfrak{K}_{1}).$
(64)

For convenience, we shall henceforth refer to $\mathfrak{K}_{1}$ as $K_{1}$
and $\mathfrak{l}$ as $l$. Applying the above results to $(2*)$ in Eq 55
(suppressing the $|m_{i}\rangle\notin D(n)$ and $|\tilde{n}\rangle\in D(n)$
sums for the sake convenience), we get

$\displaystyle(2*)$
$\displaystyle=\sum_{l=1}^{N-1}\sum_{K_{1}=1}^{N-1}\sum_{K_{2}=1}^{K_{1}-1}\cdots\sum_{K_{l}=1}^{K_{l-1}-1}\bigg{[}\begin{pmatrix}N\\\
K_{1}\end{pmatrix}\begin{pmatrix}K_{1}\\\
K_{2}\end{pmatrix}\begin{pmatrix}K_{2}\\\
K_{3}\end{pmatrix}\cdots\begin{pmatrix}K_{l-1}\\\ K_{l}\end{pmatrix}$ (65a)
$\displaystyle\qquad\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\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}$
$\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}\cdots\sum_{\widetilde{K}_{l-1}=1}^{\widetilde{K}_{l-2}-1}\bigg{[}$
$\displaystyle\qquad\begin{pmatrix}N\\\
\widetilde{K}_{1}+K_{1}\end{pmatrix}\begin{pmatrix}\widetilde{K}_{1}+K_{1}\\\
\widetilde{K}_{2}+K_{1}\end{pmatrix}\begin{pmatrix}\widetilde{K}_{2}+K_{1}\\\
\widetilde{K}_{3}+K_{1}\end{pmatrix}\cdots\begin{pmatrix}\widetilde{K}_{l-2}+K_{1}\\\
\widetilde{K}_{l-1}+K_{1}\end{pmatrix}\begin{pmatrix}\widetilde{K}_{l-1}+K_{1}\\\
K_{1}\end{pmatrix}$
$\displaystyle\qquad\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\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}.$
(65b)

Now for the last ingredient of the proof, we use the equivalent expression of
the product of the combinatorial factors:

$\displaystyle\begin{pmatrix}N\\\ \widetilde{K}_{1}+K_{1}\end{pmatrix}$
$\displaystyle\begin{pmatrix}\widetilde{K}_{1}+K_{1}\\\
\widetilde{K}_{2}+K_{1}\end{pmatrix}\begin{pmatrix}\widetilde{K}_{2}+K_{1}\\\
\widetilde{K}_{3}+K_{1}\end{pmatrix}\cdots\begin{pmatrix}\widetilde{K}_{l-2}+K_{1}\\\
\widetilde{K}_{l-1}+K_{1}\end{pmatrix}\begin{pmatrix}\widetilde{K}_{l-1}+K_{1}\\\
K_{1}\end{pmatrix}$ $\displaystyle=\begin{pmatrix}N\\\
K_{1}\end{pmatrix}\begin{pmatrix}N-K_{1}\\\
\widetilde{K}_{1}\end{pmatrix}\begin{pmatrix}\widetilde{K}_{1}\\\
\widetilde{K}_{2}\end{pmatrix}\begin{pmatrix}\widetilde{K}_{2}\\\
\widetilde{K}_{3}\end{pmatrix}\cdots\begin{pmatrix}\widetilde{K}_{l-2}\\\
\widetilde{K}_{l-1}\end{pmatrix}.$ (66)

Collecting all of the results above, we obtain for $(2*)$

$\displaystyle(2*)$
$\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}\cdots\sum_{\widetilde{K}_{l-1}=1}^{\widetilde{K}_{l-2}-1}\bigg{[}\begin{pmatrix}N\\\
K_{1}\end{pmatrix}\begin{pmatrix}N-K_{1}\\\
\widetilde{K}_{1}\end{pmatrix}\begin{pmatrix}\widetilde{K}_{1}\\\
\widetilde{K}_{2}\end{pmatrix}\begin{pmatrix}\widetilde{K}_{2}\\\
\widetilde{K}_{3}\end{pmatrix}\cdots\begin{pmatrix}\widetilde{K}_{l-2}\\\
\widetilde{K}_{l-1}\end{pmatrix}$
$\displaystyle\qquad\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\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}$
$\displaystyle=\sum_{K_{1}=1}^{N-1}\begin{pmatrix}N\\\
K_{1}\end{pmatrix}\mathcal{R}^{K_{1}}_{\tilde{n},n}\bigg{[}\mathcal{H}_{n,\tilde{n}}^{N-K_{1}}+\sum_{\widetilde{K}_{1}=1}^{N-K_{1}-1}\begin{pmatrix}N-K_{1}\\\
\widetilde{K}_{1}\end{pmatrix}\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}}^{\,}}$
$\displaystyle\qquad+\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)}\begin{pmatrix}N-K_{1}\\\
\widetilde{K}_{1}\end{pmatrix}\begin{pmatrix}\widetilde{K}_{1}\\\
\widetilde{K}_{2}\end{pmatrix}\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}}^{\,}}$
$\displaystyle\qquad+\cdots$ $\displaystyle\qquad+\cdots$
$\displaystyle\qquad+\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)}\bigg{\\{}$ $\displaystyle\qquad\qquad\begin{pmatrix}N-K_{1}\\\
\widetilde{K}_{1}\end{pmatrix}\begin{pmatrix}\widetilde{K}_{1}\\\
\widetilde{K}_{2}\end{pmatrix}\begin{pmatrix}\widetilde{K}_{2}\\\
\widetilde{K}_{3}\end{pmatrix}\cdots\begin{pmatrix}\widetilde{K}_{N-K_{1}-2}\\\
\widetilde{K}_{N-K_{1}-1}\end{pmatrix}$
$\displaystyle\qquad\qquad\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\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)$
$\displaystyle\qquad\qquad\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{\\}}\bigg{]}.$
(67)

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 Eq 56) and vanish when
$|\tilde{n}\rangle\neq|n\rangle$ (as per Eq 57). This is the case because
$N-K_{1}<N$ or equivalently $N-K_{1}\leqslant N-1$, and we have assumed that
the form of the band derivative given in Eq 56 holds up to $N-1$. Thus some of
the terms in $(2*)$ vanish while the rest add up to

$(2*)=\sum_{K_{1}=1}^{N-1}\begin{pmatrix}N\\\
K_{1}\end{pmatrix}\partial^{N-K_{1}}\varepsilon_{n}^{\,}\,\mathcal{R}_{n,n}^{K_{1}},$
(68)

which exactly cancels out with $(1*)$. As for the proof of Eq 57 at $M=N$, we
take overlap of Eq 45 with $|\tilde{n}\rangle$ instead of $|n\rangle$, 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 $\partial^{M}H$ 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 Eq 56 and 57 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$.

## References

  * [1] I. Lifshitz,  _Sov. Phys. JETP_ 1960, _11_ , 5 1130.
  * [2] L. Van Hove,  _Phys. Rev._ 1953, _89_ 1189\.
  * [3] A. Chandrasekaran, A. Shtyk, J. J. Betouras, C. Chamon,  _Phys. Rev. Research_ 2020, _2_ 013355\.
  * [4] D. V. Efremov, A. Shtyk, A. W. Rost, C. Chamon, A. P. Mackenzie, J. J. Betouras,  _Phys. Rev. Lett._ 2019, _123_ 207202\.
  * [5] A. Shtyk, G. Goldstein, C. Chamon,  _Phys. Rev. B_ 2017, _95_ 035137\.
  * [6] N. F. Q. Yuan, L. Fu,  _Phys. Rev. B_ 2020, _101_ 125120\.
  * [7] 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, 2022,  URL https://arxiv.org/abs/2205.08828.
  * [8] D. Aoki, G. Seyfarth, A. Pourret, A. Gourgout, A. McCollam, J. A. N. Bruin, Y. Krupko, I. Sheikin,  _Phys. Rev. Lett._ 2016, _116_ 037202\.
  * [9] 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_ 245139\.
  * [10] 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_ 147001\.
  * [11] S. N. Khan, D. D. Johnson,  _Phys. Rev. Lett._ 2014, _112_ 156401\.
  * [12] 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, et al.,  _npj Quantum Materials_ 2019, _4_ , 1 1.
  * [13] Y. Okamoto, A. Nishio, Z. Hiroi,  _Phys. Rev. B_ 2010, _81_ 121102\.
  * [14] Y. Sherkunov, A. V. Chubukov, J. J. Betouras,  _Phys. Rev. Lett._ 2018, _121_ 097001\.
  * [15] S. Slizovskiy, A. V. Chubukov, J. J. Betouras,  _Phys. Rev. Lett._ 2015, _114_ 066403\.
  * [16] I. Stewart,  _Reports on Progress in Physics_ 1982, _45_ , 2 185\.
  * [17] E. Yelland, J. Barraclough, W. Wang, K. Kamenev, A. Huxley,  _Nature physics_ 2011, _7_ , 11 890.
  * [18] N. F. Yuan, H. Isobe, L. Fu,  _Nature communications_ 2019, _10_ , 1 1.
  * [19] Y. Sherkunov, J. J. Betouras,  _Phys. Rev. B_ 2018, _98_ 205151\.
  * [20] L. Classen, A. V. Chubukov, C. Honerkamp, M. M. Scherer,  _Phys. Rev. B_ 2020, _102_ 125141\.
  * [21] H. Isobe, L. Fu,  _Phys. Rev. Research_ 2019, _1_ 033206\.
  * [22] J. L. McChesney, A. Bostwick, T. Ohta, T. Seyller, K. Horn, J. González, E. Rotenberg,  _Phys. Rev. Lett._ 2010, _104_ 136803\.
  * [23] P. Rosenzweig, H. Karakachian, D. Marchenko, K. Küster, U. Starke,  _Phys. Rev. Lett._ 2020, _125_ 176403\.
  * [24] 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, et al.,  _Nature Physics_ 2022, 1–8.
  * [25] Y. Hu, X. Wu, B. R. Ortiz, S. Ju, X. Han, J. Ma, N. C. Plumb, M. Radovic, R. Thomale, S. D. Wilson, et al.,  _Nature Communications_ 2022, _13_ , 1 1.
  * [26] A. Consiglio, T. Schwemmer, X. Wu, W. Hanke, T. Neupert, R. Thomale, G. Sangiovanni, D. Di Sante,  _Phys. Rev. B_ 2022, _105_ 165146\.
  * [27] T. Neupert, M. M. Denner, J.-X. Yin, R. Thomale, M. Z. Hasan,  _Nature Physics_ 2022, _18_ , 2 137.
  * [28] 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_ , 6582 774.
  * [29] A. Chandrasekaran, J. J. Betouras,  _Phys. Rev. B_ 2022, _105_ 075144\.
  * [30] E. K. Kokkinis, G. Goldstein, D. V. Efremov, J. J. Betouras,  _Phys. Rev. B_ 2022, _105_ 155123\.
  * [31] R. G. Ayoub,  _Archive for history of exact sciences_ 1980, 253–277.
  * [32] B. Bradlyn, J. Cano, Z. Wang, M. Vergniory, C. Felser, R. J. Cava, B. A. Bernevig,  _Science_ 2016, _353_ , 6299 aaf5037.
  * [33] W. R. Inc.,  Mathematica, Version 13.1,  URL https://www.wolfram.com/mathematica,  Champaign, IL, 2022.
  * [34] Maplesoft, a division of Waterloo Maple Inc..,  Maple,  URL https://hadoop.apache.org.
  * [35] F. D. M. Haldane,  _Phys. Rev. Lett._ 1988, _61_ 2015\.
  * [36] A. Chandrasekaran,  Band utilities,  https://github.com/anirudhc-git/Series-expansion-of-bandstructure,  A Mathematica package for analysing bandstructures arising from k-space models.