# Effect of disorder on density of states and conductivity in higher order Van
Hove singularities in two dimensional bands

Anirudh Chandrasekaran  Department of Physics, Loughborough University,
Loughborough LE11 3TU, UK.  Joseph J. Betouras  Department of Physics,
Loughborough University, Loughborough LE11 3TU, UK.

###### Abstract

We study systems with energy bands in two dimensions, hosting higher order Van
Hove singularities (HOVHS) in the presence of disorder, using standard
diagrammatic techniques for impurity averaging. In the clean limit, such
singularities cause power-law divergence in the density of states (DOS), and
this is expected to strongly affect electronic correlation. In order to
analyse the signatures of these singularities in disordered systems, we employ
various Born approximations, culminating in the self-consistent (non) Born
approximation. Although the divergence of the DOS is smeared, we find that the
shape of the DOS, as characterized by the power law tail and the universal
ratio of prefactors, is retained slightly away from the singularity. This
could help us to understand current and future experiments on materials that
can be tuned to host HOVHS. The impurity induced smearing is calculated and
analysed for several test cases of singularities. We also study the effects of
impurities on electrical conductivity and determine the regimes where the
quantitative features of the power law DOS manifest in the conductivity.

##  I Introduction

Novel phases of matter driven by non-trivial topology and geometry of
electronic band structure, have been the subject of much interest in recent
times. While the various effects of topology have been well studied and
documented, new avenues continue to emerge in the investigation of band
structure geometries. Pioneering early work by Lifshitz and Van Hove Lifshitz
(1960); Van Hove (1953) laid the foundation for a rich path of subsequent
explorations uncovering the exotic effects of Fermi surface geometry and its
consequences, especially in the context of Fermi surface topological
transitions. In Fermi surface topological transitions, the Fermi surface
geometry undergoes a sudden change when some parameters in the system are
changed. Lifshitz initially studied two particular forms of these topological
changes (pocket appearing/disappearing or neck formation/collapse). Lifshitz
(1960); Abrikosov (2017) Quite often, such transitions happen when the Fermi
surface hosts critical points of the dispersion.

When the gradient of the energy dispersion
$\varepsilon_{n}^{\phantom{{\dagger}}}(\mathbf{k})$ of the $n^{\mathrm{th}}$
band at some point $\mathbf{k}_{0}^{\phantom{{\dagger}}}$, represented by the
Jacobian
$\nabla\varepsilon_{n}^{\phantom{{\dagger}}}(\mathbf{k_{0}^{\phantom{{\dagger}}}})$
Vanishes, we have a critical point of the dispersion, that could be a maximum,
minimum or a saddle. It is often adequate to describe the dispersion around
these points with a Taylor expansion to quadratic order, taking a canonical
form $\pm k_{x}^{2}\pm k_{y}^{2}$. Such extrema are accompanied by a
logarithmic divergence in the DOS Van Hove (1953) at the corresponding energy.
There are however, a class of critical points around which the dispersion
needs to be Taylor expanded beyond quadratic order. These are points where the
determinant of the Hessian of dispersion relation also Vanishes, and they are
known as higher order critical points, and the corressponding dispersion
relation is said to have a higher order Van Hove singularity (HOVHS). An
example of HOVHS that is frequently reported in the literature is the cusp
singularity having the canonical dispersion $k_{x}^{4}-k_{y}^{2}$.
Chandrasekaran  _et al._ (2020); Efremov  _et al._ (2019); Shtyk  _et al._
(2017)

At a higher order critical point, the Fermi surface becomes singular as
happens in a regular Van Hove singularity. HOVHS are also accompanied by power
law diverging DOS, often with asymmetric ratio of prefactors above and below
the singular energy. Chandrasekaran  _et al._ (2020); Yuan and Fu (2020)
Signatures of this can be observed experimentally in the tunnelling
conductivity. Yuan  _et al._ (2019); Kerelsky  _et al._ (2019) This property
is expected to affect other measurable quantities such as electrical
conductivity and charge susceptibility as well. Furthermore, electronic
correlation is expected to be enhanced by the large DOS in the vicinity of the
HOVHS potentially leading to novel phases driven by electronic interaction.
The kinetic energy of the fermions in the vicinity of those points is
comparable or less to the interaction potential energy.

Lifshitz transitions and associated logarithmically-divergent Van Hove
singularities have been reported in in a variety of materials including
cuprates, iron based superconductors, cobaltates,
$\text{Sr}_{2}\text{RuO}_{4}$ and heavy fermions. (Aoki  _et al._ , 2016;
Barber  _et al._ , 2019; Benhabib  _et al._ , 2015; Khan and Johnson, 2014;
Coldea  _et al._ , 2019; Okamoto  _et al._ , 2010; Sherkunov and Betouras,
2018; Slizovskiy  _et al._ , 2015; Stewart, 1982; Yelland  _et al._ , 2011)
There is an even more recent surge of interest in higher order Van Hove
singularities. Classen  _et al._ (2020); Lin and Nandkishore (2020); Bi and Fu
(2021); Guerci  _et al._ (2021); Oriekhov  _et al._ (2021) Some of the
materials where they have been discovered include
$\text{Sr}_{3}\text{Ru}_{2}\text{O}_{7}$ where a higher order ($X_{9}$ with
$n=4$) Van Hove singularity was shown to exist in the presence of an external
magnetic field, (Efremov  _et al._ , 2019; Chandrasekaran  _et al._ , 2020)
while a different higher order Van Hove saddle has been reported in highly
overdoped graphene (Rosenzweig  _et al._ , 2020) and may be quite releVant for
the recently observed phases of Bernal bilayer graphene. Zhou  _et al._ (2021)

Given the host of exotic physical phenomena that HOVHS promise, it is not a
surprise that that HOVHS typically require delicate tuning of parameters in
the system to obtain (effected through strain, pressure, twist, bias voltage,
etc). This has lately been made possible with the adVances in experimental
techniques. Although infinitely many such distinct singularities exist,
catastrophe theory Castrigiano and Hayes (2019); Poston and Stewart (2014);
Bruce  _et al._ (1992) guarantees that we are typically likely to obtain only
a finite subset of singularities in real systems. For a thorough
classification of such HOVHS in two dimensional systems refer Chandrasekaran
_et al._ , 2020; Yuan and Fu, 2020. Mathematically, these singularities are
unstable to certain classes of perturbations that cause the higher order
critical points to break into a number of ordinary critical points. Thus,
where it is theoretically possible to obtain a HOVHS, it becomes pertinent to
investigate the effect of such factors as impurities, that may impede their
experimental detection.

A random distribution of impurities in a crystalline system with a
concentration $n_{\text{imp}}^{\phantom{{\dagger}}}$, however low, clearly
distorts the lattice periodicity, rendering a simplistic application of the
Bloch-band idea questionable. Nevertheless, we expect that for low
concentrations, a perturbative, diagrammatic treatment over the clean band
theoretic system would be sufficient. There a few such diagrammatic techniques
available to treat the problem of randomly distributed impurities, such as
quenched averaging, the replica trick, supersymmetry and the Keyldish
technique. Keldysh  _et al._ (1965); Kamenev and Andreev (1999); Efetov
(1999); Mézard  _et al._ (1987) In trying to understand the impact of
disorder, one typically tries to compute how certain quantities like the
spectral function $A_{\mathbf{k}}^{\phantom{{\dagger}}}(\omega)$ and the
electrical conductivity $\sigma(\omega)$ are affected due to perturbative
corrections. The spectral function gives information about the density of
states which is a measurable quantity (alongside response functions like
various conductivities). For the clean system, these quantities show singular
behaviour due to the HOVHS. It is therefore important that we estimate both
the qualitative and quantitative changes in the presence of impurities.

Intuitively, we would expect the signatures of the HOVHS to survive to some
extent when the impurity scattering is ‘weak’. One way to characterize weak
scattering is a low concentration of impurities, as compared to concentration
of electrons or atoms. (Bruus and Flensberg, 2004) More precisely, consider a
system with area (or volume) $\mathcal{A}$ having
$N_{\text{el}}^{\phantom{{\dagger}}}$ electrons and
$N_{\text{imp}}^{\phantom{{\dagger}}}$ impurities so that the concentration of
impurities is
$n_{\text{imp}}^{\phantom{{\dagger}}}=N_{\text{imp}}^{\phantom{{\dagger}}}/\mathcal{A}$
and concentration of electrons is
$n_{\text{el}}^{\phantom{{\dagger}}}=N_{\text{el}}^{\phantom{{\dagger}}}/\mathcal{A}$.
We then need $n_{\text{imp}}^{\phantom{{\dagger}}}\ll
n_{\text{el}}^{\phantom{{\dagger}}}$.

Another important quantity needed to characterise the weakness of impurity
scattering is the strength of a single impurity potential denoted by
$V(\mathbf{x})$. Here we assume that it is short ranged and denote its average
value in a unit cell centred around the impurity as $\overline{V}$. Since the
zeroth Fourier transform $u_{\mathbf{q}=0}^{\phantom{{\dagger}}}$ is the
integral of this short ranged potential over all space, we can express the
unit cell average of the potential in terms of it as $\overline{V}\approx
u_{0}^{\phantom{{\dagger}}}/A_{\text{u.c}}^{\phantom{{\dagger}}}$, where
$A_{\text{u.c}}^{\phantom{{\dagger}}}$ is the area of the real lattice unit
cell. If $E_{s}^{\,}$ is some energy scale corresponding to the singular
dispersion (see Appendix A for a strategy to choose $E_{s}^{\,}$), then we
would expect the impact of the impurity scattering on the singularity to be
weak when $\overline{V}\ll E_{s}^{\,}$, or equivalently

$u_{0}^{\phantom{{\dagger}}}\ll
A_{\text{u.c}}^{\phantom{{\dagger}}}\,E_{s}^{\,}.$ (1)

This estimate is independent of the condition of validity for the diagrammatic
full Born approximation to be derived below. It is used to carefully choose
the numbers for the numerical calculations that will follow.

In this paper, we use the diagrammatic quenched averaging technique to
determine the smearing of the power law diverging DOS. To this end, we shall
calculate the self energy, and through it the scattering lifetime of electrons
in a disordered system hosting a HOVHS in the clean limit. We begin in Sec II,
by detailing the hierarchy of approximations that we will use to calculate the
impurity averaged self energy in a pedagogical manner. We are careful to apply
the continuum limit only at the end of calculations, in order to avoid
mathematical complications such as analytic continuation of integrals as
opposed to finite sums. This is important due to the delicate nature of the
problem being studied. After analysing the smearing of the DOS within the
various approximations, we proceed to briefly discuss the consequences for
electrical conductivity due to HOVHS and impurities in Sec III. We then
summarize and conclude in Sec IV.

##  II Various approximations

Figure 1: The various approximations to the renormalized one particle Greens
functions are represented above diagrammatically. The simplest is the tree
level approximation shown in (a). Going to one loop we get the first Born
approximation (1BA) in (b). The full Born approximation consists of all
$O(n_{\text{imp}})$ diagrams. With respect to the self-consistent Born
approximation, there two distinct approximations that are referred to by this
name. The simpler one involves diagram (c) while the inclusion of all the
succeeding diagrams (which we refer to as the self-consistent Born
approximation in this paper) is sometimes referred to as self-consistent non-
Born approximation.

Since the self energy can not be calculated to arbitrary orders, we employ a
series of tractable approximations that retain only particular diagrams or
some entire classes of diagrams. The various approximations used below are
depicted in Fig 1. We briefly note the Feynman rules: dashed impurity lines
carrying momentum $\mathbf{q}$ contribute an amplitude
$u_{\mathbf{q}}^{\phantom{{\dagger}}}$ (the Fourier coefficient of a single
impurity potential), the solid Fermion lines carrying
$(iq_{n}^{\phantom{{\dagger}}},\mathbf{q})$ are accompanied by the
corresponding bare Fermionic propagator
$\mathcal{G}_{\mathbf{q}}^{0}(iq_{n}^{\phantom{{\dagger}}})$, the unfilled
diamond shaped impurity vertices each contribute a factor of
$n_{\text{imp}}^{\phantom{{\dagger}}}$, and momentum is conserved at all
vertices. Only internal Fermion momenta are summed over. Note that since
impurities are elastic scatterers, no frequency flows through the impurity
lines. The precise origin and motivation for these rules can be found in Ref
Bruus and Flensberg, 2004.

Before we proceed further, we explain some aspects of the notation used in
this work. We use $\varepsilon_{n}^{\phantom{{\dagger}}}(\mathbf{k})$ or
simply $\varepsilon(\mathbf{k})$ to denote the full and actual band dispersion
without reference to the Fermi level. In contrast,
$\xi^{\phantom{{\dagger}}}_{\mathbf{k}}=\varepsilon_{n}^{\phantom{{\dagger}}}(\mathbf{k})-\mu$
is the full or series expanded dispersion, adjusted by the chemical potential
$\mu$, i.e. at zero temperature, states with negative
$\xi^{\phantom{{\dagger}}}_{\mathbf{k}}$ are occupied while states with
positive $\xi^{\phantom{{\dagger}}}_{\mathbf{k}}$ are unoccupied. For $\mu=0$,
the Fermi level lies exactly at the higher order singularity. The Fermionic
Matsubara frequencies are denoted by
$k_{n}^{\phantom{{\dagger}}}=2(n+1)\pi/\beta$.

###  II.1 Tree level

At tree level, we only have a single impurity vertex and a single line as in
Fig 1 (a). The tree level self energy is then simply

$\Sigma_{\mathbf{k}}^{\phantom{{\dagger}}}(ik_{n}^{\phantom{{\dagger}}})=n_{\text{imp}}^{\phantom{{\dagger}}}\,u_{0}^{\phantom{{\dagger}}}.$
(2)

This contribution is purely real and $\mathbf{k}$-independent and provides
only a constant shift to the dispersion.

###  II.2 First Born approximation

This is the first loop correction to the self energy, shown in Fig 1 (b). It
takes the value

$\Sigma_{\mathbf{k}}^{\text{1BA}}(ik_{n}^{\phantom{{\dagger}}})=\frac{n_{\text{imp}}^{\phantom{{\dagger}}}}{\mathcal{A}}\sum_{\mathbf{q}}\left|u_{\mathbf{k}-\mathbf{q}}^{\phantom{{\dagger}}}\right|^{2}\,\frac{1}{ik_{n}^{\phantom{{\dagger}}}-\xi_{\mathbf{q}}^{\phantom{{\dagger}}}}.$
(3)

It is the first non-trivial $O(n_{\text{imp}}^{\phantom{{\dagger}}})$ loop
correction. Since the sum over $\mathbf{q}$ is finite in a finite subsystem,
we analytically continue to obtain the retarded self energy as

$\Sigma_{\mathbf{k}}^{\text{1BA}}(\omega)=\frac{n_{\text{imp}}^{\phantom{{\dagger}}}}{\mathcal{A}}\sum_{\mathbf{q}}\left|u_{\mathbf{k}-\mathbf{q}}^{\phantom{{\dagger}}}\right|^{2}\,\frac{1}{\omega-\xi_{\mathbf{q}}^{\phantom{{\dagger}}}+i\delta}.$
(4)

Now, we assume that $u_{\mathbf{k}-\mathbf{q}}^{\phantom{{\dagger}}}\approx
u_{0}^{\phantom{{\dagger}}}$, convert the $\mathbf{q}$ sum to an energy
integral with density of states $g(\epsilon)$ and take imaginary part in order
to obtain the scattering lifetime

$\displaystyle\frac{1}{2\tau_{\mathbf{k}}^{\phantom{{\dagger}}}(\omega)}$
$\displaystyle=-\text{Im}\,\Sigma_{\mathbf{k}}^{\text{1BA}}(\omega),$ (5a)
$\displaystyle=\pi\,n_{\text{imp}}\,u_{0}^{2}\int\limits_{\infty}^{\infty}d\epsilon\,g(\epsilon)\,\delta(\omega-\epsilon+\mu),$
$\displaystyle=\pi\,n_{\text{imp}}\,u_{0}^{2}\,g(\omega+\mu),$ (5b)

where we have used the density of states per unit area (or volume) or simply
density of states (DOS), defined as

$\displaystyle g(\epsilon)=$
$\displaystyle\frac{1}{\mathcal{A}}\sum_{\mathbf{q}}\delta(\epsilon-\xi_{\mathbf{q}}^{\phantom{{\dagger}}}),$
(6a) $\displaystyle\approx$
$\displaystyle\int\frac{d^{2}q}{(2\pi)^{2}}\,\delta(\epsilon-\xi_{\mathbf{q}}^{\phantom{{\dagger}}}).$
(6b)

The density of states will play an important role in our calculations. It is
discussed in greater detail in Appendix B. When $\omega,\mu\rightarrow 0$, the
inverse lifetime diverges due to the singular density of states of higher
order singularities. The full Born approximation contains _all_ the
$O(n_{\text{imp}}^{\phantom{{\dagger}}})$ contributions to the self energy,
which, diagrammatically speaking, arise from scattering through a single
impurity (and therefore a single $n_{\text{imp}}^{\phantom{{\dagger}}}$
contribution).

###  II.3 Interlude: the $t$-matrix

We will need the $t$-matrix for the full Born approximation as it contains
information about the scattering of an electron across a _single_ impurity.
The $t$-matrix is defined diagrammatically in Fig 1 (within the curly braces
$\\{\,\\}$). Momentum conservation between ingoing and outgoing momenta is not
demanded. Therefore the $t$-matrix has two momenta indices:
$t_{\mathbf{k}_{1},\mathbf{k}_{2}}^{\phantom{{\dagger}}}(ik_{n}^{\phantom{{\dagger}}})$.
As before, only one frequency dependence is needed since the impurity
potentials scatter elastically. It is easy to show that the $t$-matrix
satisfies a self consistent equation

$\displaystyle
t_{\mathbf{k}_{1},\mathbf{k}_{2}}^{\phantom{{\dagger}}}(ik_{n}^{\phantom{{\dagger}}})=$
$\displaystyle\,n_{\text{imp}}^{\phantom{{\dagger}}}\,u_{0}^{\phantom{{\dagger}}}\,\delta_{\mathbf{k}_{1},\mathbf{k}_{2}}^{\phantom{{\dagger}}}$
$\displaystyle+\frac{1}{\mathcal{A}}\sum_{\mathbf{q}}u_{\mathbf{k}_{1}-\mathbf{q}}^{\phantom{{\dagger}}}\,\mathcal{G}^{0}_{\mathbf{q}}(ik_{n}^{\phantom{{\dagger}}})\,t_{\mathbf{q},\mathbf{k}_{2}}^{\phantom{{\dagger}}}(ik_{n}^{\phantom{{\dagger}}}).$
(7)

Now, if we assume that the impurity amplitudes
$u_{\mathbf{q}}^{\phantom{{\dagger}}}$ are effectively momentum independent,
i.e $u_{\mathbf{q}}^{\phantom{{\dagger}}}\approx u_{0}^{\phantom{{\dagger}}}$,
then we immediately see that the right hand side of the above equation does
not give any non-trivial dependence on $\mathbf{k}_{1}^{\phantom{{\dagger}}}$
since the $\mathbf{k}_{1}^{\phantom{{\dagger}}}$ dependence is expressed
through the factor $u_{\mathbf{k}_{1}-\mathbf{q}}^{\phantom{{\dagger}}}$ which
is now set to a constant. We can thus assume that the $t$ matrix is also
momentum independent, which greatly simplifies calculations.

Analytically continuing
$ik_{n}^{\phantom{{\dagger}}}\rightarrow\omega+i\delta$ and assuming
$u_{\mathbf{k}-\mathbf{q}}^{\phantom{{\dagger}}}\approx
u_{0}^{\phantom{{\dagger}}}$ so that the $t$-matrix depends only on $\omega$
and $\mu$ (with momentum conservation implied) we obtain:

$t(\omega,\mu)=\frac{n_{\text{imp}}^{\phantom{{\dagger}}}\,u_{0}^{\phantom{{\dagger}}}}{1-u_{0}^{\phantom{{\dagger}}}\frac{1}{\mathcal{A}}\sum_{\mathbf{q}}\mathcal{G}_{\mathbf{q}}^{0}(\omega)}.$
(8)

We have made the $\mu$ dependence of the $t$ matrix explicit above. Once
again, the inverse area weighted momentum sums can be converted into integrals
over energy levels weighted by the DOS. The integrals take the following form
that is releVant for both the full Born (FBA) and self consistent Born (SBA)
approximations

$I(\omega,\mu,z_{0}^{\phantom{{\dagger}}})=\int_{-\infty}^{\infty}d\epsilon\,g(\epsilon)\,\frac{1}{\omega-\epsilon+\mu+z_{0}^{\phantom{{\dagger}}}},$
(9)

where $\text{Im}\,z_{0}^{\phantom{{\dagger}}}\neq 0$ and the density of states
is given by

$g(\epsilon)=\begin{cases}D_{+}^{\phantom{{\dagger}}}|\epsilon|^{-\nu},&\epsilon>0\\\
D_{-}^{\phantom{{\dagger}}}|\epsilon|^{-\nu},&\epsilon<0\end{cases}\;.$ (10)

The integral is evaluated in Appendix D. The final result reads

$I(\omega,\mu,z_{0}^{\phantom{{\dagger}}})=\frac{2\pi
i(\omega+\mu+z_{0}^{\phantom{{\dagger}}})^{-\nu}}{1-e^{-i2\pi\nu}}\,[e^{-i\pi\nu}D_{-}^{\phantom{{\dagger}}}-D_{+}^{\phantom{{\dagger}}}].$
(11)

In evaluating the above integral, the principal branch is used for the
fractional power, i.e $(Re^{i\theta})^{-\nu}=R^{-\nu}e^{-i\nu\theta}$ for
$\theta\in(0,2\pi)$. Now for the $t$ matrix, we have
$z_{0}^{\phantom{{\dagger}}}=+i\delta$. This implies that

$(x+i\delta)^{-\nu}=\begin{cases}x^{-\nu}(1-i\delta)&x>0\\\
|x|^{-\nu}e^{-i\pi\nu}(1+i\delta)&x<0\end{cases},$ (12)

where for us $x=\omega+\mu$. Defining
$r=D_{+}^{\phantom{{\dagger}}}/D_{-}^{\phantom{{\dagger}}}$, we have

$\displaystyle I(\omega,\mu,i\delta)=$ $\displaystyle
2\pi\,g(\omega+\mu)\,\,\frac{i(e^{-i\pi\nu}-r)}{1-e^{-i2\pi\nu}}$
$\displaystyle\times\begin{cases}\frac{1}{r}(1-i\delta)&\omega+\mu>0\\\
e^{-i\pi\nu}(1+i\delta)&\omega+\mu<0\end{cases}.$ (13)

The $t$ matrix is then given by

$t(\omega,\mu)=\frac{n_{\text{imp}}^{\phantom{{\dagger}}}\,u_{0}^{\phantom{{\dagger}}}}{1-u_{0}^{\phantom{{\dagger}}}\,I(\omega,\mu,i\delta)}.$
(14)

For concreteness, we include the $t$ matrix values for a few of the
singularities that have been identified in lattice models in Table 1.

Singularity | $D_{+}^{\phantom{{\dagger}}}/D_{-}^{\phantom{{\dagger}}}$ | $\nu$ | $t(\omega,\mu)$
---|---|---|---
Fold | $1/\sqrt{3}$ | $1/6$ | $\dfrac{n_{\text{imp}}^{\phantom{{\dagger}}}\,u_{0}^{\phantom{{\dagger}}}}{1-\pi u_{0}^{\phantom{{\dagger}}}\,g(\omega+\mu)\,\,\left(\left(\sqrt{3}\right)^{\text{sgn}(\omega+\mu)}-i\right)}$
Cusp | $1/\sqrt{2}$ | $1/4$ | $\dfrac{n_{\text{imp}}^{\phantom{{\dagger}}}\,u_{0}^{\phantom{{\dagger}}}}{1-\pi u_{0}^{\phantom{{\dagger}}}\,g(\omega+\mu)\,\,\left(\Theta(\omega+\mu)-i\right)}$
| Monkey
---
saddle
$1$ | $1/3$ | $\dfrac{n_{\text{imp}}^{\phantom{{\dagger}}}\,u_{0}^{\phantom{{\dagger}}}}{1-\pi u_{0}^{\phantom{{\dagger}}}\,g(\omega+\mu)\,\,\left(\frac{1}{\sqrt{3}}\,\text{sgn}(\omega+\mu)-i\right)}$
$X_{9}^{\phantom{{\dagger}}}$ | $1$ | $1/2$ | $\dfrac{n_{\text{imp}}^{\phantom{{\dagger}}}\,u_{0}^{\phantom{{\dagger}}}}{1-\pi u_{0}^{\phantom{{\dagger}}}\,g(\omega+\mu)\,\,(1-i\,\text{sgn}(\omega+\mu))}$
Table 1: The values of the $t$-matrix for some of the commonly observed
singularities, computed assuming momentum independence of the Fourier
transform of the impurity potential (i.e $u_{\mathbf{q}}^{\,}\approx
u_{0}^{\,}$). This ensures that the $t$-matrix depends only on frequency. The
$t$-matrix values tabulated above give the self energy within the full Born
approximation, that in turn yields the lifetime that smears the DOS. For
precise values of $D_{+}^{\,}$ and $D_{-}^{\,}$ of these singularities, see
Table 2.

###  II.4 Full Born approximation

The full Born approximation (FBA) is the truncation of the perturbation
expansion of
$\Sigma_{\mathbf{k}}^{\phantom{{\dagger}}}(ik_{n}^{\phantom{{\dagger}}})$ to
order $n_{\text{imp}}$, i.e it includes all the scattering across a single
impurity. But this is the information contained in the $t$-matrix as we saw
above. More precisely

$\displaystyle\Sigma_{\mathbf{k}}^{\text{FBA}}(ik_{n}^{\phantom{{\dagger}}})$
$\displaystyle=$ $\displaystyle
n_{\text{imp}}^{\phantom{{\dagger}}}\,u_{0}^{\phantom{{\dagger}}}+\frac{1}{\mathcal{A}}\sum_{\mathbf{q}}u_{\mathbf{k}-\mathbf{q}}^{\phantom{{\dagger}}}\,\mathcal{G}_{\mathbf{q}}^{0}(ik_{n}^{\phantom{{\dagger}}})\,t_{\mathbf{q},\mathbf{k}}^{\phantom{{\dagger}}}(ik_{n}^{\phantom{{\dagger}}})$
(15) $\displaystyle=$ $\displaystyle
t_{\mathbf{k},\mathbf{k}}(ik_{n}^{\phantom{{\dagger}}}).$

Thus, in the Born approximation, the self energy is simply a diagonal element
of the $t$-matrix. By assuming that
$u_{\mathbf{q}}^{\phantom{{\dagger}}}\approx u_{0}^{\phantom{{\dagger}}}$, the
$t$-matrix becomes momentum independent and the FBA self energy coincides with
the $t$-matrix:

$\displaystyle\Sigma_{\mathbf{k}}^{\text{FBA}}(\omega,\mu)=$
$\displaystyle\,t(\omega,\mu),$ (16a) $\displaystyle=$
$\displaystyle\,\frac{n_{\text{imp}}^{\phantom{{\dagger}}}\,u_{0}^{\phantom{{\dagger}}}}{1-u_{0}^{\phantom{{\dagger}}}\,I(\omega,\mu,i\delta)}.$
(16b)

###  II.5 Self-consistent Born approximation

This is an improvement over Born approximation wherein the full Green’s
function
$\mathcal{G}_{\mathbf{q}}^{\phantom{0}}(ik_{n}^{\phantom{{\dagger}}})$ is used
in place of the bare one
$\mathcal{G}_{\mathbf{q}}^{0}(ik_{n}^{\phantom{{\dagger}}})$, in the diagrams
of the Born approximation leading to a self consistent relation

$\displaystyle\Sigma_{\mathbf{k}}^{\phantom{{\dagger}}}(i$ $\displaystyle
k_{n}^{\phantom{{\dagger}}})$ $\displaystyle=$
$\displaystyle\,n_{\text{imp}}^{\phantom{{\dagger}}}\,u_{0}^{\phantom{{\dagger}}}+\frac{1}{\mathcal{A}}\sum_{\mathbf{q}}u_{\mathbf{k}-\mathbf{q}}^{\phantom{{\dagger}}}\,\mathcal{G}_{\mathbf{q}}(ik_{n}^{\phantom{{\dagger}}})\,\tilde{t}_{\mathbf{q},\mathbf{k}}^{\phantom{{\dagger}}}(ik_{n}^{\phantom{{\dagger}}}),$
(17a) $\displaystyle=$
$\displaystyle\,n_{\text{imp}}^{\phantom{{\dagger}}}\,u_{0}^{\phantom{{\dagger}}}$
$\displaystyle+\frac{1}{\mathcal{A}}\sum_{\mathbf{q}}u_{\mathbf{k}-\mathbf{q}}^{\phantom{{\dagger}}}\,\frac{1}{ik_{n}^{\phantom{{\dagger}}}-\xi_{\mathbf{q}}^{\phantom{{\dagger}}}-\Sigma_{\mathbf{q}}^{\phantom{{\dagger}}}(ik_{n}^{\phantom{{\dagger}}})}\,\tilde{t}_{\mathbf{q},\mathbf{k}}^{\phantom{{\dagger}}}(ik_{n}^{\phantom{{\dagger}}}).$
(17b)

Here, we introduce the modified version of the $t$ matrix with the full
$\mathcal{G}_{\mathbf{q}}^{\phantom{0}}(ik_{n}^{\phantom{{\dagger}}})$ in the
place of $\mathcal{G}_{\mathbf{q}}^{0}(ik_{n}^{\phantom{{\dagger}}})$ in its
expansion. The self consistent relation for this $\tilde{t}$-matrix is

$\displaystyle\tilde{t}_{\mathbf{k}_{1},\mathbf{k}_{2}}^{\phantom{{\dagger}}}$
$\displaystyle(ik_{n}^{\phantom{{\dagger}}})$ $\displaystyle=$
$\displaystyle\,n_{\text{imp}}^{\phantom{{\dagger}}}\,u_{0}^{\phantom{{\dagger}}}\,\delta_{\mathbf{k}_{1},\mathbf{k}_{2}}^{\phantom{{\dagger}}}$
$\displaystyle+\frac{1}{\mathcal{A}}\sum_{\mathbf{q}}u_{\mathbf{k}_{1}-\mathbf{q}}^{\phantom{{\dagger}}}\,\mathcal{G}^{\phantom{0}}_{\mathbf{q}}(ik_{n}^{\phantom{{\dagger}}})\,\tilde{t}_{\mathbf{q},\mathbf{k}_{2}}^{\phantom{{\dagger}}}(ik_{n}^{\phantom{{\dagger}}}),$
(18a) $\displaystyle=$
$\displaystyle\,n_{\text{imp}}^{\phantom{{\dagger}}}\,u_{0}^{\phantom{{\dagger}}}\,\delta_{\mathbf{k}_{1},\mathbf{k}_{2}}^{\phantom{{\dagger}}}$
$\displaystyle+\frac{1}{\mathcal{A}}\sum_{\mathbf{q}}u_{\mathbf{k}_{1}-\mathbf{q}}^{\phantom{{\dagger}}}\,\,\frac{1}{ik_{n}^{\phantom{{\dagger}}}-\xi_{\mathbf{q}}^{\phantom{{\dagger}}}-\Sigma_{\mathbf{q}}^{\phantom{{\dagger}}}(ik_{n}^{\phantom{{\dagger}}})}\,\tilde{t}_{\mathbf{q},\mathbf{k}_{2}}^{\phantom{{\dagger}}}(ik_{n}^{\phantom{{\dagger}}}).$
(18b)

The self energy, again, is simply the diagonal element of the $\tilde{t}$
matrix:
$\Sigma_{\mathbf{k}}^{\text{SBA}}(ik_{n}^{\phantom{{\dagger}}})=\tilde{t}_{\mathbf{k},\mathbf{k}}(ik_{n}^{\phantom{{\dagger}}})$.
We now analytically continue
$ik_{n}^{\phantom{{\dagger}}}\rightarrow\omega+i\delta$ and assume that both
the $\tilde{t}$ matrix and the self energy depend only on frequency. We can
then write

$\displaystyle\Sigma^{\text{SBA}}(\omega,\mu)=$
$\displaystyle\,\tilde{t}(\omega,\mu),$ (19a) $\displaystyle=$
$\displaystyle\,\frac{n_{\text{imp}}^{\phantom{{\dagger}}}\,u_{0}^{\phantom{{\dagger}}}}{1-u_{0}^{\phantom{{\dagger}}}\frac{1}{\mathcal{A}}\sum_{\mathbf{q}}\mathcal{G}_{\mathbf{q}}^{\phantom{0}}(\omega)},$
(19b) $\displaystyle=$
$\displaystyle\,\frac{n_{\text{imp}}^{\phantom{{\dagger}}}\,u_{0}^{\phantom{{\dagger}}}}{1-u_{0}^{\phantom{{\dagger}}}\,I(\omega,\mu,-\Sigma^{\text{SBA}}(\omega,\mu))}.$
(19c)

After substituting for $I$, and writing $\nu=m/n$ for integer $m$ and $n$, we
can rearrange the resultant expression to obtain a polynomial equation for
$\Sigma^{\text{SBA}}(\omega,\mu)\equiv\Sigma$:

$\left(\omega+\mu-\Sigma\right)^{m}\Sigma^{n}-\frac{1}{(2\pi
iu_{0}^{\phantom{{\dagger}}})^{n}}\left(\frac{1-e^{-i2\pi\nu}}{e^{-i\pi\nu}D_{-}^{\phantom{{\dagger}}}-D_{+}^{\phantom{{\dagger}}}}\right)^{n}\\\
\times\left(\Sigma-
n_{\text{imp}}^{\phantom{{\dagger}}}u_{0}^{\phantom{{\dagger}}}\right)^{n}=0.$
(20)

However, it is hard to pick the right root for this equation, even if we were
able to solve it in particular instances, therefore we adopt an alternate
strategy. From the branch chosen to define the function $f(z)=z^{-\nu}$, we
have that $\text{arg}\,f(z)\in(-2\pi\nu,0)$. We can then define
$x=\left(\omega+\mu-\Sigma\right)^{-\nu}$. Clearly we need to choose $x$ such
that $\text{arg}\,x\in(-2\pi\nu,0)$. For simplicity let us treat the case
where $m=1$. We can get a polynomial equation for $x$:

$(\omega+\mu)\,\varphi\,x^{n+1}+(n_{\text{imp}}^{\phantom{{\dagger}}}\,u_{0}^{\phantom{{\dagger}}}-\omega-\mu)\,x^{n}-\varphi\,x+1=0,$
(21)

where we have defined a constant characteristic of the singularity

$\varphi=2\pi
i\,(u_{0}^{\phantom{{\dagger}}}\,D_{-}^{\phantom{{\dagger}}})\left(\frac{e^{-i\pi\nu}-r}{1-e^{-i2\pi\nu}}\right).$
(22)

In general, we will have to solve this equation numerically and pick the root
$x_{0}^{\phantom{{\dagger}}}$ with
$\text{arg}\,x_{0}^{\phantom{{\dagger}}}\in(-2\pi\nu,0)$.

To compute the smeared DOS, we express it in terms of the spectral function:

$\displaystyle\tilde{g}(\omega)=$
$\displaystyle\frac{1}{\mathcal{A}}\sum_{\mathbf{k}}A_{\mathbf{k}}^{\phantom{{\dagger}}}(\omega),$
(23a) $\displaystyle=$
$\displaystyle\frac{1}{2\pi\mathcal{A}}\sum_{\mathbf{k}}\frac{1/\tau_{\mathbf{k}}^{\phantom{{\dagger}}}(\omega)}{(\omega-\xi_{\mathbf{k}}^{\phantom{{\dagger}}})^{2}+1/4\tau_{\mathbf{k}}^{2}(\omega)},$
(23b) $\displaystyle\approx$
$\displaystyle\frac{1}{2\pi}\int\frac{d^{2}k}{(2\pi)^{2}}\,\frac{1/\tau_{\mathbf{k}}^{\phantom{{\dagger}}}(\omega)}{(\omega-\xi_{\mathbf{k}}^{\phantom{{\dagger}}})^{2}+1/4\tau_{\mathbf{k}}^{2}(\omega)},$
(23c) $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}d\epsilon\,\,g(\epsilon)\frac{1/\tau(\omega)}{(\omega-\epsilon)^{2}+1/4\tau^{2}(\omega)}.$
(23d)

Using Eq 5a, where $\tau_{\mathbf{k}}^{\phantom{{\dagger}}}(\omega)$ is
defined, we can numerically compute the smeared DOS within the various
approximations outlined above. We plot the smeared DOS for the cusp and monkey
saddle singularities in Fig 2.

Figure 2: The smeared density of states (DOS) for the monkey saddle and cusp
singularities, within the first Born (1BA), full Born (FBA) and self-
consistent Born approximations (SBA) are compared to the DOS of the bare
singularity in panels (a) and (c) respectively. We see that the different
approximations coincide at large energies $\epsilon$. Near the singularity,
the 1BA DOS goes to zero, the FBA DOS diverges, while the SBA alone captures
the finiteness of the smeared singularity. In the panels (b) and (d), we
replot these approximations in $\log-\log$ coordinates to make the power law
dependence explicit. The asymmetric nature of the cusp singularity manifests
as the non-trivial ratio $g(-|\epsilon|)/g(|\epsilon|)=\sqrt{2}$. The
rationale behind the choice of numerical parameters used for generating the
plot is explained in Appendix C.

###  II.6 Validity of the expansion

Before we proceed to analyse the various Born approximations outlined above,
we briefly discuss the validity of the diagrammatic technique. We do this in
the context of the full Born approximation since the first Born and tree level
are contained within this. Furthermore, the precise topology of the
$m^{\text{th}}$ order Born diagram is well known, making it possible to
estimate its size. In contrast, the self consistent Born evades a careful
analysis of such sorts.

The $m^{\text{th}}$ order Born diagram contains precisely one impurity vertex
$n_{\text{imp}}^{\phantom{{\dagger}}}$, $m$ internal Fermion lines and $m+1$
impurity lines (see the first row of Fig 1). The impurity lines simply
contribute respective factors $u_{\mathbf{q}}^{\phantom{{\dagger}}}\sim
u_{0}^{\phantom{{\dagger}}}$ while the Fermion lines each entail a momentum
sum that can be approximated by a momentum integral. Similar to the procedure
shown above for the first Born approximation, each such integral roughly
yields a factor $g(\mu)$ so that the $m^{\text{th}}$ order Born diagram can be
estimated as being proportional to
$n_{\text{imp}}^{\phantom{{\dagger}}}\,u_{0}^{\phantom{{\dagger}}}\,(u_{0}^{\phantom{{\dagger}}}\,g(\mu))^{m}$.
It is easy to check that this has a dimension $\mathcal{E}$ of energy since
$[n_{\text{imp}}^{\phantom{{\dagger}}}]=l^{-2}$,
$[u_{0}^{\phantom{{\dagger}}}]=\mathcal{E}\,l^{2}$ and
$[g]=\mathcal{E}^{-1}\,l^{-2}$ (since $u_{0}^{\phantom{{\dagger}}}$, the
zeroth Fourier coefficient is simply the integral of the scattering potential
over all space and we use the density of states per unit area. See Ref Bruus
and Flensberg, 2004 for more details). As the energy factor
$n_{\text{imp}}^{\phantom{{\dagger}}}\,u_{0}^{\phantom{{\dagger}}}$ is common
to all diagrams, the dimensionless parameter determining the validity of the
expansion is $u_{0}^{\phantom{{\dagger}}}\,g(\mu)$. Thus, we would expect the
full Born approximation to be valid in the regime
$u_{0}^{\phantom{{\dagger}}}\,g(\mu)\ll 1$. Although we are unable to make
such precise estimates in the case of self consistent Born approximation, we
can compute the hierarchy of energies for each of the above approximations,
that is the energy beyond which SBA coincides with FBA, and the energy beyond
which FBA coincides with 1BA. This is done below.

###  II.7 Analysis of the self energy in self consistent Born approximation

The self energy in the SBA satisfies the following relation as seen from the
discussion above

$\Sigma^{\text{SBA}}(\mu)=\frac{n_{\text{imp}}^{\phantom{{\dagger}}}\,u_{0}^{\phantom{{\dagger}}}}{1-\varphi\,(\mu-\Sigma^{\text{SBA}}(\mu))^{-\nu}},$
(24)

with the constant $\varphi$ defined in Eq 22. When we are in the regime where
$|\mu|\gg|\Sigma^{\text{SBA}}(\mu)|$, we can essentially ignore
$\Sigma^{\text{SBA}}(\mu)$ to write
$(\mu-\Sigma^{\text{SBA}}(\mu))^{-\nu}\approx\mu^{-\nu}$. Let us focus on
$\mu>0$. It is straightforward to extend the discussion to $\mu<0$. We now
pull out a factor $D_{+}^{\phantom{{\dagger}}}=r\,D_{-}^{\phantom{{\dagger}}}$
from $\varphi$ to obtain

$\Sigma^{\text{SBA}}(\mu)\approx\frac{n_{\text{imp}}^{\phantom{{\dagger}}}\,u_{0}^{\phantom{{\dagger}}}}{1-u_{0}^{\phantom{{\dagger}}}\,\zeta\,g(\mu)},$
(25)

with the constant $\zeta$ being characteristic of the singularity and
independent of impurity concentration and strength

$\zeta=\frac{2\pi i}{r}\left(\frac{e^{-i\pi\nu}-r}{1-e^{-i2\pi\nu}}\right).$
(26)

This approximate expression for self energy coincides with that of the full
Born approximation so that the condition for SBA to pass into FBA is
$|\mu|\gg|\Sigma^{\text{SBA}}(\mu)|$. The inverse lifetime is now obtained by
taking the imaginary part

$\displaystyle\tau^{-1}(\mu)=$
$\displaystyle\,-2\,\text{Im}\,\Sigma^{\text{SBA}}(\mu),$ (27a)
$\displaystyle=$
$\displaystyle\frac{2\,n_{\text{imp}}^{\phantom{{\dagger}}}\,u_{0}^{2}\,g(\mu)\,\,\text{Im}\,\zeta^{*}}{\left|1-u_{0}^{\phantom{{\dagger}}}\,\zeta\,g(\mu)\right|^{2}}.$
(27b)

Now for the standard singularities, we can verify explicitly that $\zeta$ is a
small number taking a numerical value near $1$ (for example,
$|\zeta|=\sqrt{2}\,\pi$ for the cusp). Therefore, when we have
$u_{0}^{\phantom{{\dagger}}}\,g(\mu)\ll 1$, we can ignore it in the
denominator and we have

$\displaystyle\tau^{-1}(\mu)\approx$
$\displaystyle\,(2\,n_{\text{imp}}^{\phantom{{\dagger}}}\,u_{0}^{2}\,\,\text{Im}\,\zeta^{*})\,g(\mu),$
(28a) $\displaystyle=$
$\displaystyle\,(2\,n_{\text{imp}}^{\phantom{{\dagger}}}\,u_{0}^{2}\,\pi)\,g(\mu).$
(28b)

To get to the last step we use the fact that $\text{Im}\,\zeta^{*}=\pi$ for
any possible value of $\nu$ and $r$. This is the inverse lifetime given by the
first Born approximation as well, as can be checked by comparing to Eq 5.
Therefore, the condition $u_{0}^{\phantom{{\dagger}}}\,g(\mu)\ll 1$ allows us
to transition from FBA to 1BA. Notice that this condition was also shown to
validate the FBA in Sec II.6. To summarize, we have

$\text{SBA}\xrightarrow{|\mu|\gg|\Sigma^{\text{SBA}}(\mu)|}\text{FBA}\xrightarrow{u_{0}^{\phantom{{\dagger}}}\,g(\mu)\ll
1}\text{1BA}.$ (29)

Furthermore, in the 1BA regime, $\tau^{-1}(\mu)\propto g(\mu)$.

Let us examine the case $u_{0}^{\phantom{{\dagger}}}\,g(\mu)\gg 1$ but with
$|\mu|\gg|\Sigma^{\text{SBA}}(\mu)|$ in some detail. While applicability of
the FBA may be questionable here, we rely on the SBA for analysis. Since
$u_{0}^{\phantom{{\dagger}}}\,g(\mu)\gg 1$, we ignore the $1$ in the
denominator of Eq 27b and obtain

$\tau^{-1}(\mu)\approx\frac{2\,n_{\text{imp}}^{\phantom{{\dagger}}}\,\,\text{Im}\,\zeta^{*}}{|\zeta|^{2}\,g(\mu)}=\frac{2\,n_{\text{imp}}^{\phantom{{\dagger}}}\,\pi}{|\zeta|^{2}\,g(\mu)}.$
(30)

Thus, in the regime where $|\mu|\gg|\Sigma^{\text{SBA}}(\mu)|$ and
$u_{0}^{\phantom{{\dagger}}}\,g(\mu)\gg 1$ the lifetime, rather than the
inverse lifetime has the signature of the DOS, i.e $\tau(\mu)\propto g(\mu)$.
Our numerical studies indicate that the first condition is always satisfied
even slightly above $\mu=0$. In fact the real and imaginary parts of the SBA
self energy approach finite constants in the large $\mu$ limit. Furthermore,
the finiteness and smallness of $\text{Im}\,\Sigma^{\text{SBA}}$ in the
respective regimes of interest is consistent with Eq 28b and Eq 30.

Therefore, depending on the strength of the scattering potential with respect
to the density of states, either the resistivity (proportional to $\tau^{-1}$)
or conductivity (proportional to $\tau$) can show quantitative signatures of
the underlying singularity. Notice that this result is _independent_ of the
concentration of the impurities. While we may object that the condition
$u_{0}^{\phantom{{\dagger}}}\,g(\mu)\gg 1$ renders FBA inapplicable, we have
relied here on the SBA rather than FBA to draw this conclusion. We should be
cautious in carrying over the validity estimates of FBA, an essentially
perturbative approximation, over to the SBA, that is a partly non-perturbative
treatment. The precise relation between conductivity/resistivity and the
lifetime is elaborated in the section below.

##  III Conductivity

To compute the conductivity using the Kubo formula for tight-binding models,
we use the current operator defined as

$\hat{\mathbf{j}}=\frac{e}{\mathcal{A}}\sum_{\mathbf{k}}\nabla_{\mathbf{k}}^{\phantom{{\dagger}}}\xi_{\mathbf{k}}^{\phantom{{\dagger}}}\,\hat{c}_{\mathbf{k}}^{{\dagger}}\hat{c}_{\mathbf{k}}.$
(31)

Although this lacks momentum or spatial dependence, we can use this current to
obtain the essential features of the DC conductivity, including the conjoined
effects of the singularity and impurities. To simplify the notation, we have
made only the charge explicit in the above expression and have suppressed
other dimensionfull factors (like mass). These can be reinstated later when
necessary.

An essential component of this calculation is the retarded current correlation
function that takes the form Bruus and Flensberg (2004); Mahan (2000)

$\displaystyle\Pi_{ij}^{R}(iq_{n}^{\phantom{{\dagger}}})=$
$\displaystyle\frac{e^{2}}{\mathcal{A}}\sum_{\mathbf{p}}(\nabla_{\mathbf{p}}^{\phantom{{\dagger}}}\xi_{\mathbf{p}}^{\phantom{{\dagger}}})_{i}^{\phantom{{\dagger}}}(\nabla_{\mathbf{p}}^{\phantom{{\dagger}}}\xi_{\mathbf{p}}^{\phantom{{\dagger}}})_{j}^{\phantom{{\dagger}}}$
$\displaystyle\times\frac{1}{\beta}\sum_{m}\mathcal{G}(\mathbf{p},ip_{m}^{\phantom{{\dagger}}}+iq_{n}^{\phantom{{\dagger}}})\,\mathcal{G}(\mathbf{p},ip_{m}^{\phantom{{\dagger}}}),$
(32)

where $\mathcal{A}$ is the area (or volume) and $i$ and $j$ are the directions
indicating the component of the gradient to be used. The dc conductivity is
then obtained by using the Lehmann representation for the Green’s functions
(in terms of the spectral function), performing the Matsubara sums,
analytically continuing
$iq_{n}^{\phantom{{\dagger}}}\rightarrow\omega+i\delta$, and evaluating the
$\omega\rightarrow 0$ limit to get Mahan (2000)

$\displaystyle\sigma_{ij}=e^{2}\int\frac{d^{2}p}{(2\pi)^{2}}\,\Bigg{[}$
$\displaystyle\frac{\partial\xi_{\mathbf{p}}^{\phantom{{\dagger}}}}{\partial
p_{i}^{\phantom{{\dagger}}}}\,\frac{\partial\xi_{\mathbf{p}}^{\phantom{{\dagger}}}}{\partial
p_{j}^{\phantom{{\dagger}}}}$
$\displaystyle\times\int\frac{d\epsilon}{2\pi}A^{2}_{\mathbf{p}}(\epsilon)\left(-\frac{dn_{\text{F}}^{\phantom{{\dagger}}}(\epsilon)}{d\epsilon}\right)\Bigg{]}.$
(33)

Let us assume that the inverse of the scattering lifetime
$\tau_{\mathbf{p}}^{-1}(\omega,\mu)$ is finite and small. If we further assume
that it is approximately independent of $\mathbf{p}$ and that the dependence
on $\omega$ is “weak” in comparison the strong $\omega^{2}$ of the spectral
function $A_{\mathbf{p}}^{\phantom{{\dagger}}}(\omega)$, we can altogether
drop $\omega$ and $\mathbf{p}$ dependence so as to justify a relabelling as
$\tau^{-1}(\mu)$ (since it still depends on the Fermi level). We briefly
elaborate on the latter point, that is the weak $\omega$-dependence of
$(\tau^{-1}(\omega,\mu))^{2}$ as compared to $\omega^{2}$, both of which
appear in the denominator of the spectral function. We first note that within
each of the Born approximations discussed above, $(\tau^{-1}(\omega,\mu))^{2}$
depends on $\omega$ and $\mu$ in the same way (more precisely, it depends on
$\omega+\mu$). From panels (b), (d), (f) and (h) of Fig 4, we see that as
$\omega$ (or $\mu$ in the figure) changes from $0.1$ to $1$ in the numerical
calculations, $\tau_{\mathbf{p}}^{-1}(\omega,\mu)$ at the most doubles in
value. Thus, a change of $\omega^{2}$ by two orders of magnitude corresponds
to a fourfold change in $(\tau^{-1}(\omega,\mu))^{2}$, which still has a value
much smaller than unity. Using this property, we drop the $\omega$ dependence
of $\tau^{-1}(\omega,\mu)$. The spectral function then effectively becomes a
Cauchy distribution. We can then use the properties of Cauchy distribution and
the delta function to approximate the above expression as Mahan (2000)

$\sigma_{ij}^{\phantom{{\dagger}}}=2e^{2}\,\tau(\mu)\int\frac{d^{2}p}{(2\pi)^{2}}\,\frac{\partial\xi_{\mathbf{p}}^{\phantom{{\dagger}}}}{\partial
p_{i}^{\phantom{{\dagger}}}}\,\frac{\partial\xi_{\mathbf{p}}^{\phantom{{\dagger}}}}{\partial
p_{j}^{\phantom{{\dagger}}}}\,\delta(\xi_{\mathbf{p}}^{\phantom{{\dagger}}}-\mu).$
(34)

At this point we have just recovered the conventional expression for
conductivity that is proportional to the electron lifetime at the Fermi level.
Nevertheless, the main outstanding issue is to ensure that the integral is
finite in the $\mu\rightarrow 0$ limit for HOVHS, since we are dealing with
systems that have infrared divergences in the clean limit. We do this in
Appendix E, where we show that the integral is finite and non-zero in the
$\mu\rightarrow 0$ limit.

Figure 3: Within the self-consistent Born approximation, the impurity smeared
density of states is finite at the singularity, but has a power law tail
resembling the bare singularity. The resemblance to the bare singularity
improves as the concentration of impurities $n_{\text{imp}}^{\,}$ and the
strength of the scattering potential $u_{0}^{\,}$ become smaller. This is
depicted in the figure above where, in panels (a) and (c), we fix
$n_{\text{imp}}^{\,}$ and vary $u_{0}^{\,}$, while in panels (b) and (d), we
fix $u_{0}^{\,}$ and vary $n_{\text{imp}}^{\,}$. We observe that the smeared
DOS and bare DOS become indistinguishable at either low concentrations or weak
scattering strength. The choice of numerical parameters used is explained in
Appendix C. We note that the condition $u_{0}^{\,}\,g(\mu)\ll 1$, mentioned in
Sec II.7, is trivially satisfied in the situations treated above, with
$\epsilon$ instead of $\mu$. Figure 4: Signatures of the power law diverging
DOS can manifest in the inverse life-time $\tau^{-1}(\mu)$ of the electrons at
the Fermi level. The primary requirement for this, derived in Sec II.7, is
that the self energy within the self-consistent Born approximation has to be
much smaller than the chemical potential $\mu$, that measures the energy
distance of the Fermi level from the singularity. Mathematically this reads
$|\Sigma^{\text{SBA}}(\mu)/\mu|\ll 1$. In addition to this, when the product
$u_{0}^{\,}\,g(\mu)$, of the ‘strength’ of the scattering potential and the
density of states at the Fermi level is much smaller than unity (see Fig 3),
the inverse lifetime of the electrons at the Fermi level becomes approximately
proportional to $g(\mu)$. In the pairs of panels (a), (b) and (c), (d), we fix
the impurity concentration and vary $u_{0}^{\,}$ in a system hosting the
monkey saddle, perturbed by impurities. While the condition for power-law
signatures to appear in $\tau^{-1}(\mu)$ is analysed graphically in the left
panels (a) and (c), the inverse lifetime, weighed by appropriate factors is
compared to the DOS in the right panels (b) and (d). A similar scheme with
$u_{0}^{\,}$ fixed and varying $n_{\text{imp}}^{\,}$ is shown in the pairs of
panels (e), (f) and (g), (h). We see that the signatures of the power law DOS
are seen well in $\tau^{-1}(\mu)$, particularly as either $u_{0}^{\,}$ or
$n_{\text{imp}}^{\,}$ become small. The tight binding model where the
singularity is obtained, including the choice of numerical of parameters, is
discussed in Appendix C.

##  IV Results and discussion

In the previous sections, we have outlined the mathematical procedure involved
in analysing higher order singularities in two dimensional bands, in the
presence of impurities. We have elucidated various practical approximations
(the hierarchy of Born approximations) that make an analytic and numerical
analysis tractable. Furthermore, keeping in mind the requirements of numerical
calculations and the physical analysis of weak impurity scattering made in the
preceding sections, we worked out a scheme to choose and/or analyse the
impurity concentration and strength in the context of lattice tight binding
models (see Sec I and Appendices A, C). This broad framework allows us to
treat various lattice models hosting HOVHS and perform pertinent calculations.
In this section, we analyse the results of such a calculation.

The primary quantity of interest when dealing with a band singularity is the
DOS, which provides a measurable diagnosis of the singularity. To this end, we
analyse the impurity induced smearing of DOS due in two HOVHS, the cusp
singularity (with dispersion $k_{x}^{4}-k_{y}^{2}$) and the monkey saddle
singularity (having dispersion
$k_{y}^{3}-3\,k_{y}^{\phantom{{\dagger}}}\,k_{x}^{2}$). While an analysis of
all of the seventeen singularities part of the catastrophe theory
classification (Chandrasekaran  _et al._ , 2020) is possible, we restrict the
presentation to these two for brevity, since the qualitative features are
similar. In fact, the monkey saddle being an odd, three-fold rotationally
symmetric singularity with particle-hole symmetry and symmetric DOS and, the
cusp being an even, two-fold rotationally symmetric singularity having
asymmetric DOS and no particle hole symmetry are reasonable representatives of
the rest of the singularities.

From Fig 2, we see that the smeared DOS under various Born approximations
coincides with the bare DOS at large energies (away from the singularity)
while somewhat closer to the singularity, the full Born and self-consistent
Born coincide and very close to the singularity, only the SBA has a finite
peak. All this sits well with the analysis of Sec II.7. As the $\log-\log$
plot reveals, the power-law tail of the DOS and the asymmetric ratio of
prefactors are present even in the smeared DOS (of the various approximations)
under suitable conditions. This important fact addresses one of our primary
questions as to whether any features of the singularity will survive impurity
averaging. Furthermore, as Fig 3 reveals, the smeared DOS in SBA reflects the
bare DOS better and better as the concentration and scattering strength of the
impurities decreases.

In Sec II.7, we laid out the conditions under which the inverse lifetime of
the electrons at the Fermi level may be proportional to the DOS. This
condition along with the inverse lifetime is analysed graphically in Fig 4. We
observe that so long as the impurity strength and concentration remain small,
the inverse lifetime does become proportional to the bare DOS slightly away
from the singularity. Since the inverse lifetime plays an important role in
the electrical conductivity, we subsequently analysed the conductivity in Sec
III. As mentioned above, the integrals multiplying the lifetime in the
expression for conductivity (in Eq 34) are shown to be finite and non-zero in
the $\mu\rightarrow 0$ limit in Appendix E, allowing us to conclude that the
signatures of the singularity may potentially appear in the conductivity as
well.

In this paper, we made a diagrammatic analysis of the problem of higher order
singularities in unclean systems, using well known approximations. We found
that in so far as impurity scattering is weak, it allows much of the
quantitative signatures of the HOVHS to survive in measurable physical
quantities such as DOS and electrical conductivity. While potential comparison
to real materials and experiments may be possible even at this level of
description, a logical and necessary next step would be to go beyond the
approximations applied in this work. The effect on other response functions
like magnetic and thermal susceptibilities is left for future work. An
investigation of such a nature would entail the discussion of the combined
effects of both disorder and electronic interactions, and this may lead to
novel phases driven by instabilities, apart from non-trivial renormalizations
of the Fermi surface geometry itself.

###### Acknowledgements.

We would like to thank Claudio Chamon, Siddhant Das, Dima Efremov, Mark
Greenaway, Garry Goldstein and Alex Shtyk for useful discussions. The work has
been supported by the EPSRC grant EP/T034351/1.

##  Appendix A Choosing the energy scale $E_{s}^{\,}$

The energy scale $E_{s}^{\,}$ is chosen to give a rough energy window around
the singularity, where the polynomial dispersion of the pristine singularity
(and its power law DOS) describe the original dispersion ‘adequately’. There
are a number of reasonable strategies we can employ to choose this. An obvious
strategy is to look for other critical points around the higher order critical
point. These must necessarily exist since the dispersion is periodic. The
minimum of the absolute difference between the energy of the singularity and
the energies of other critical points, provides a natural scale at which the
polynomial singularity fails to describe the actual dispersion.
Mathematically, if the singularity is located at
$\mathbf{k}_{0}^{\phantom{{\dagger}}}$,

$E_{s}^{\,}=\inf\limits_{\begin{subarray}{c}\mathbf{k}^{\prime},\,\text{s.
t}\\\
\nabla\xi_{\mathbf{k}^{\prime}}^{\,}=0\end{subarray}}\left|\xi_{\mathbf{k}_{0}}^{\phantom{{\dagger}}}-\xi_{\mathbf{k}^{\prime}}^{\phantom{{\dagger}}}\right|.$
(35)

A simplification of this procedure is to restrict the search to only the high
symmetry points of the BZ (These are already constrained by symmetry to be
critical points. See Ref Chandrasekaran  _et al._ , 2020).

##  Appendix B Discussion on the density of states

###  B.1 Scaling and DOS

The power law dependence of the DOS in a higher order singularity can be
obtained by rescaling the momentum integration variables
$k_{x}^{\phantom{{\dagger}}}$ and $k_{y}^{\phantom{{\dagger}}}$ in the
following integral

$g(\epsilon)=\int\limits_{\text{u.c}}\frac{d^{2}k}{(2\pi)^{2}}\,\,\delta(\xi_{\mathbf{k}}-\epsilon),$
(36)

where the integration is over a unit cell centred at the singularity, and the
precise scaling transformation $k_{x}^{\phantom{{\dagger}}}=\epsilon^{p}u$ and
$k_{y}^{\phantom{{\dagger}}}=\epsilon^{q}v$ is characteristic of the
respective singularities. Doing so we obtain the power law divergent DOS
taking the form $g(\epsilon)\sim\epsilon^{p+q-1}$.

Such a scaling transformation of course modifies the domain of integration.
But the leading order divergent term can be obtained by extending the domain
of integration over the entire $k$-plane, and by doing so we pick up only a
finite error (Chandrasekaran  _et al._ , 2020). For the prefactor above and
below the singular energy, we still have to evaluate the integral. (See Refs
Chandrasekaran  _et al._ , 2020 and Yuan and Fu, 2020 for such calculations).
For the pristine forms of the singularities (i.e the canonical forms, such as
$k_{x}^{4}-k_{y}^{2}$ etc), the precise values of the coefficients can been
evaluated. In real systems however, we expect the each of the terms in the
dispersion to have non-trivial coefficients, for example
$\xi_{\mathbf{k}}^{\phantom{{\dagger}}}=\alpha\,k_{x}^{4}-\beta\,k_{y}^{2}$.
Let us assume for simplicity that $\alpha$ and $\beta$ are positive.
Appropriate modifications can be made in the forthcoming procedure when either
or both of them are negative. We first do a slightly different scale
transformation $k_{x}^{\phantom{{\dagger}}}=\alpha^{-p}\epsilon^{p}\,u$ and
$k_{y}^{\phantom{{\dagger}}}=\beta^{-q}\epsilon^{q}\,v$ and the density of
states for this case, denoted by $\tilde{g}(\epsilon)$, differs from the DOS
for the pristine singularity $g(\epsilon)$ by simply an overall constant, i.e

$\tilde{g}(\epsilon)=\alpha^{-p}\beta^{-q}\,g(\epsilon).$ (37)

We now state the exact expressions for the DOS of the singularities treated in
the paper, near the critical energy in Table 2.

Singularity | Dispersion ($\alpha,\beta>0$) | $g(\epsilon)$
---|---|---
Fold | $\alpha\,k_{x}^{3}-\beta k_{y}^{2}$ | $\alpha^{-1/3}\,\beta^{-1/2}\,\frac{2\,\Gamma\left(1/3\right)}{2\,\pi^{3/2}\Gamma\left(5/6\right)}\left(\Theta(-\epsilon)+\frac{\Theta(\epsilon)}{\sqrt{3}}\right)\,|\epsilon|^{-1/6}$
Cusp | $\alpha\,k_{x}^{4}-\beta k_{y}^{2}$ | $\alpha^{-1/4}\,\beta^{-1/2}\,\frac{2\,\Gamma^{2}\left(5/4\right)}{\pi^{5/2}}\left(\Theta(-\epsilon)+\frac{\Theta(\epsilon)}{\sqrt{2}}\right)\,|\epsilon|^{-1/4}$
| Monkey
---
saddle
$\alpha(k_{x}^{3}-3\,k_{x}^{\phantom{{\dagger}}}\,k_{y}^{2})$ | $\alpha^{-2/3}\,\frac{\Gamma\left(7/6\right)}{2\,\pi^{3/2}\,\Gamma\left(2/3\right)}|\epsilon|^{-1/3}$
$X_{9}^{\phantom{{\dagger}}}$ | $\alpha(k_{x}^{4}+k_{y}^{4}-6\,k_{x}^{2}\,k_{y}^{2})$ | $\alpha^{-1/2}\,\frac{\Gamma\left(1/4\right)}{16\,\pi^{3/2}\,\Gamma\left(3/4\right)}|\epsilon|^{-1/2}$
Table 2: The density of states per unit area for some of the commonly
occurring higher order singularities are tabulated above. We have used a
slightly generic form of the ideal dispersions, with arbitrary coefficients
$\alpha$ and $\beta$ (We have treated the case where both $\alpha$ and $\beta$
are positive for the sake of simplicity. The discussion can be extended easily
to other situations). The presence of non-trivial coefficients simply modifies
the power law diverging DOS by an overall constant. The symbol $\Gamma$
specified above is the familiar gamma function that extends the notion of
factorials. Keeping in mind the realistic situations where the monkey saddle
and $X_{9}^{\,}$ occur respectively at 3-fold and 4-fold rotationally
symmetric points, we have used general expressions that respect these
symmetries.

###  B.2 Signature of DOS: $\tau(\mu)$ vs $\tau^{-1}(\mu)$

In Sec II.7, we laid out the conditions to be satisfied so that either
$\tau(\mu)$ or $\tau^{-1}(\mu)$ may show signatures of the power law diverging
DOS. Here we work out the energy window where the signatures may be observed.
For $\tau(\mu)\sim g(\mu)$, we need
$u_{0}^{\phantom{{\dagger}}}\,\tilde{g}(\mu)\gg 1$ so that we have

$\mu_{\text{min}}^{\phantom{{\dagger}}}\ll\mu\ll\text{min}\left\\{\left(\frac{D_{+}^{\phantom{{\dagger}}}\,u_{0}^{\phantom{{\dagger}}}}{\alpha^{p}\beta^{q}}\right)^{1/\nu},E_{s}^{\,}\right\\}.$
(38)

For $\tau^{-1}(\mu)\sim g(\mu)$ we instead need

$\text{max}\left\\{\mu_{\text{min}}^{\phantom{{\dagger}}},\left(\frac{D_{+}^{\phantom{{\dagger}}}\,u_{0}^{\phantom{{\dagger}}}}{\alpha^{p}\beta^{q}}\right)^{1/\nu}\right\\}\ll\mu\ll
E_{s}^{\,},$ (39)

where $E_{s}^{\,}$ is the energy scale described in Appendix A, where the
description in terms of the pure singularity breaks down. Lastly,
$\mu_{\text{min}}^{\phantom{{\dagger}}}$ is the scale at which

$\mu_{\text{min}}^{\phantom{{\dagger}}}\approx\left|\Sigma^{\text{SBA}}(\mu_{\text{min}}^{\phantom{{\dagger}}})\right|.$
(40)

(Recall that $\mu_{\text{min}}^{\phantom{{\dagger}}}\ll\mu$ was needed to
transition out of the SBA regime into the FBA regime).

##  Appendix C Choosing parameters for numerics

In order to perform numerical calculations, we need to carefully choose
numerical values of $n_{\text{imp}}^{\phantom{{\dagger}}}$ and
$u_{0}^{\phantom{{\dagger}}}$. To do so, we first obtain the monkey saddle and
cusp singularities by tuning tight binding models, and then use the parameters
of the models to choose $n_{\text{imp}}^{\phantom{{\dagger}}}$ and
$u_{0}^{\phantom{{\dagger}}}$. (Recall from the discussion in Sec I, that we
require $n_{\text{imp}}^{\phantom{{\dagger}}}\ll
n_{\text{el}}^{\phantom{{\dagger}}}$ and $u_{0}^{\phantom{{\dagger}}}\ll
A_{\text{u.c}}^{\phantom{{\dagger}}}\,E_{s}^{\,}$, where
$A_{\text{u.c}}^{\phantom{{\dagger}}}$ is the area of the real lattice unit
cell).

###  C.1 Monkey saddle

A single monkey saddle can be obtained in the Haldane model (Haldane, 1988)
defined on the hexagonal lattice. The nearest neighbor 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}.$ (41)

Here we have set the lengths to unity. 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:
$H(\mathbf{k})=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

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

This Hamiltonian 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}}.$
(43)

Consider the following, three-fold rotation consistent high symmetry point

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

We perform the following tuning of the staggered chemical potential

$M\rightarrow\frac{t_{1}^{2}-18\,t_{2}^{2}}{2\,\sqrt{3}\,t_{2}^{\phantom{{\dagger}}}}.$
(45)

Under this, the Taylor expanded upper band $(n=2)$ dispersion around the
$K_{-}^{\phantom{{\dagger}}}$ point reads

$\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)$
$\displaystyle+O\left(k^{4}\right).$ (46)

This is clearly the monkey saddle that we sought to obtain. Let us choose
$t_{1}^{\phantom{{\dagger}}}=1$ and $t_{2}^{\phantom{{\dagger}}}=1/2$. By
comparing the energy at $\mathbf{K}_{-}^{\phantom{{\dagger}}}$ with the energy
at the other high symmetry points viz. $\mathbf{K}_{+}^{\phantom{{\dagger}}}$,
$M$ and $\Gamma$ points, we can calculate $E_{s}^{\,}$, which we find to be
approximately $1.68$. The area of the unit cell is $3\sqrt{3}/2$ so that we
need

$u_{0}^{\phantom{{\dagger}}}\ll
A_{\text{u.c}}^{\phantom{{\dagger}}}\,E_{s}^{\phantom{{\dagger}}}\approx
4.36.$ (47)

With two atoms per unit cell, the atomic density is $4/(3\sqrt{3})\approx
0.77$ and we should choose $n_{\text{imp}}^{\phantom{{\dagger}}}\ll 0.77$. For
generating data for panels (a) and (b) in Fig 2, we used
$n_{\text{imp}}^{\phantom{{\dagger}}}=0.1$ and
$u_{0}^{\phantom{{\dagger}}}=1$. The overall coefficient $\alpha$ that
multiplies the singularity was defined earlier in Appendix B.1. In the present
case, it is found to be $\alpha=3\sqrt{3}/4$.

###  C.2 Cusp

One of the simplest models that yields the cusp singularity is a one band
tight binding model defined on a square lattice, with asymmetric $x$ and $y$
nearest and next nearest neighbor hoppings $t_{1x}^{\phantom{{\dagger}}}$,
$t_{1y}^{\phantom{{\dagger}}}$, $t_{2x}^{\phantom{{\dagger}}}$ and
$t_{2y}^{\phantom{{\dagger}}}$. The dispersion takes the form

$\displaystyle\varepsilon(\mathbf{k})=$
$\displaystyle-2\,t_{1x}^{\phantom{{\dagger}}}\cos
k_{x}^{\phantom{{\dagger}}}-2\,t_{1y}^{\phantom{{\dagger}}}\cos
k_{y}^{\phantom{{\dagger}}}$
$\displaystyle-2\,t_{2x}^{\phantom{{\dagger}}}\cos
2k_{x}^{\phantom{{\dagger}}}-2t_{2y}^{\phantom{{\dagger}}}\cos
2k_{y}^{\phantom{{\dagger}}}.$ (48)

By setting $t_{2y}^{\phantom{{\dagger}}}=t_{1y}^{\phantom{{\dagger}}}/4$ and
$t_{2x}^{\phantom{{\dagger}}}=0$, and series expanding around the point
$\mathbf{X}=(0,\pi)$, we obtain the cusp singularity

$\displaystyle\varepsilon(\mathbf{k}+\mathbf{X})\approx\,$
$\displaystyle\frac{1}{8}\left(-15+6\,k_{x}^{2}-2\,t_{1y}^{\phantom{{\dagger}}}\,\left(-6+k_{y}^{4}\right)\right)$
$\displaystyle+O(k^{6}).$ (49)

Let us choose $t_{1y}^{\phantom{{\dagger}}}=3$, so that we get
$\alpha=\beta=3/4$. Computing the energy difference between the $\mathbf{X}$
point and the $M$ and $\Gamma$ points (i.e $(\pi,\pi)$ and $(0,0)$), we find
$E_{s}^{\,}=4$. The density of atoms is unity so that we need to choose
$n_{\text{imp}}^{\phantom{{\dagger}}}\ll 1$ and
$u_{0}^{\phantom{{\dagger}}}\ll 4$. For generating panels (b) and (d) of Fig
2, we set $n_{\text{imp}}^{\phantom{{\dagger}}}=0.1$ and
$u_{0}^{\phantom{{\dagger}}}=0.75$.

##  Appendix D Evaluation of $I(\omega,\mu,z_{0}^{\,})$

Figure 5: To evaluate the integral in Eq 9, we use the contour depicted above.
It has four distinct pieces: two circular contours respectively of radius
$\mathfrak{r}\rightarrow 0$ and $R\rightarrow\infty$ and two contours slightly
above and below the positive $x$ axis. It can be shown that the circular
contours have zero contribution in their respective limits. The contours along
the positive $x$ axis are separated by the branch cut and can be added and
expressed in terms of the original integral of interest. We then apply residue
theorem to evaluate this.

To evaluate the integral in Eq 9, we split it into two parts, one for positive
$\epsilon$ and one for negative $\epsilon$ and convert the negative integral
by a substitution $z=-\epsilon$ to obtain

$\displaystyle I(\omega,\mu,z_{0}^{\phantom{{\dagger}}})=$ $\displaystyle-
D_{+}^{\phantom{{\dagger}}}\int_{0}^{\infty}dz\,z^{-\nu}\,\frac{1}{z-(\omega+\mu+z_{0}^{\phantom{{\dagger}}})}$
$\displaystyle+D_{-}^{\phantom{{\dagger}}}\int_{0}^{\infty}dz\,z^{-\nu}\,\frac{1}{z+(\omega+\mu+z_{0}^{\phantom{{\dagger}}})}.$
(50)

Since $\nu$ is a rational fraction, we need a branch cut to evaluate these
integrals. We choose the $+x$ axis for this purpose and use the contour shown
in Fig 5. Since the integrals respectively have poles at
$z_{0}^{\phantom{{\dagger}}}+\omega+\mu$ and
$-(z_{0}^{\phantom{{\dagger}}}+\omega+\mu)$, enclosed by the contours, we
apply the residue theorem yielding ultimately

$I(\omega,\mu,z_{0}^{\phantom{{\dagger}}})=\frac{2\pi
i(\omega+\mu+z_{0}^{\phantom{{\dagger}}})^{-\nu}}{1-e^{-i2\pi\nu}}\,[e^{-i\pi\nu}D_{-}^{\phantom{{\dagger}}}-D_{+}^{\phantom{{\dagger}}}].$
(51)

Notice that this is justified because
$\text{Im}\,z_{0}^{\phantom{{\dagger}}}\neq 0$ by our assumption and the poles
will lie away from the $x$-axis, safely enclosed by the contour.

##  Appendix E Finiteness of the integral

###  E.1 Cusp

Figure 6: The constant energy contours of the cusp singularity, for a positive
and negative energy are respectively depicted in (a) and (b). The region with
energy below the singularity is shaded blue, while the region with energy
above the singularity is shaded red. The cutoff procedure we will use is a box
cutoff with $\Lambda$ for $k_{x}^{\,}$ and $\Lambda^{2}$ for $k_{y}^{\,}$.

The cusp singularity has a canonical dispersion
$\xi_{\mathbf{k}}^{\phantom{{\dagger}}}=k_{x}^{4}-k_{y}^{2}$. To evaluate the
integral in Eq 34 we will use a box boundary with a
$k_{x}^{\phantom{{\dagger}}}$\- cutoff $\Lambda$ and a corresponding
$k_{y}^{\phantom{{\dagger}}}$ cutoff $\Lambda^{2}$ as shown in Fig 6. First
let us consider the case $\mu>0$. Here we can restrict the integration to the
the regions I and III. The integral then takes the form

$I_{ij}^{\phantom{{\dagger}}}(\mu,\Lambda)=\int\limits_{-\Lambda}^{\Lambda}\frac{dk_{x}^{\phantom{{\dagger}}}}{(2\pi)}\int\limits_{-k_{x}^{2}}^{k_{x}^{2}}\frac{dk_{y}^{\phantom{{\dagger}}}}{(2\pi)}(\partial_{i}^{\phantom{{\dagger}}}\xi)(\partial_{j}^{\phantom{{\dagger}}}\xi)\,\delta\left(k_{x}^{4}-k_{y}^{2}-\mu\right).$
(52)

The delta function can be expressed as

$\displaystyle\delta\bigg{(}k_{x}^{4}$ $\displaystyle-
k_{y}^{2}-\mu\bigg{)}=\Theta\left(k_{x}^{4}-\mu\right)$
$\displaystyle\times\left[\frac{\delta\left(k_{y}^{\phantom{{\dagger}}}-\sqrt{k_{x}^{4}-\mu}\right)}{|-2k_{y}^{\phantom{{\dagger}}}|}+\frac{\delta\left(k_{y}^{\phantom{{\dagger}}}+\sqrt{k_{x}^{4}-\mu}\right)}{|-2k_{y}^{\phantom{{\dagger}}}|}\right]$
(53)

The step function ensures that the roots are real (they always lie in the
integration domain when real). Now consider $i=x$ and $j=x$. We can easily
evaluate the $k_{y}^{\phantom{{\dagger}}}$ integral using the above delta
function identity to obtain

$\displaystyle I_{xx}^{\phantom{{\dagger}}}(\mu,\Lambda)=$
$\displaystyle\int\limits_{-\Lambda}^{\Lambda}\frac{dk_{x}^{\phantom{{\dagger}}}}{(2\pi)^{2}}\frac{16\,k_{x}^{6}}{\sqrt{k_{x}^{4}-\mu}}\,\Theta\left(k_{x}^{4}-\mu\right)$
(54a) $\displaystyle=$
$\displaystyle\frac{8}{\pi^{2}}\int\limits_{\mu^{1/4}}^{\Lambda}dk_{x}^{\phantom{{\dagger}}}\frac{k_{x}^{6}}{\sqrt{k_{x}^{4}-\mu}}.$
(54b)

At this point it is clear why we needed the UV cutoff $\Lambda$. Without it,
the integral in Eq 54b diverges. We now make the substitution
$t=\mu/k_{x}^{4}$. The integral becomes

$I_{xx}^{\phantom{{\dagger}}}=-\frac{2\mu^{5/4}}{\pi^{2}}\int\limits_{1}^{\mu/\Lambda^{4}}dt\,\,t^{-9/4}(1-t)^{-1/2}.$
(55)

This integral can be easily rewritten in the form of an incomplete Beta
function, using the substitution $l=1-t$:

$\displaystyle I_{xx}^{\phantom{{\dagger}}}=$
$\displaystyle\frac{2\mu^{5/4}}{\pi^{2}}\int\limits_{0}^{1-\mu/\Lambda^{4}}dl\,\,l^{1/2-1}(1-l)^{-5/4-1}$
(56a) $\displaystyle=$
$\displaystyle\frac{2\mu^{5/4}}{\pi^{2}}B_{1-\mu/\Lambda^{4}}^{\phantom{{\dagger}}}(1/2,-5/4).$
(56b)

We now series expand this around $\mu=0$:

$\displaystyle
I_{xx}^{\phantom{{\dagger}}}(\mu,\Lambda)\approx\,\Lambda^{5}\bigg{(}$
$\displaystyle 0.16+0.41\,\mu/\Lambda^{4}-0.29\,(\mu/\Lambda^{4})^{5/4}$
$\displaystyle+O\left((\mu/\Lambda^{4})^{2}\right)\bigg{)}.$ (57)

Thus, in the $\mu\rightarrow 0+$ limit, the integral goes to a finite constant
with a small linear correction to leading order. For $i=x$ and $j=y$ (or vice
versa), the sum over the two roots
$k_{y}^{\phantom{{\dagger}}}=-\sqrt{k_{x}^{4}-\mu}$ and
$k_{y}^{\phantom{{\dagger}}}=+\sqrt{k_{x}^{4}-\mu}$, causes the $k_{y}$
integral to Vanish since
$(\partial_{x}^{\phantom{{\dagger}}}\xi)(\partial_{y}^{\phantom{{\dagger}}}\xi)=-8k_{x}^{3}k_{y}^{\phantom{{\dagger}}}$
is odd in $k_{y}^{\phantom{{\dagger}}}$. Therefore we continue to investigate
the case $i=y$, $j=y$. After applying the identity in Eq 62 and simplifying,
the integral in this case reads

$\displaystyle I_{yy}^{\phantom{{\dagger}}}(\mu,\Lambda)=$
$\displaystyle\frac{1}{\pi^{2}}\int\limits_{-\Lambda}^{\Lambda}dk_{x}^{\phantom{{\dagger}}}\,\sqrt{k_{x}^{4}-\mu}$
(58a) $\displaystyle=$
$\displaystyle\frac{2}{\pi^{2}}\int\limits_{\mu^{1/4}}^{\Lambda}dk_{x}^{\phantom{{\dagger}}}\,\sqrt{k_{x}^{4}-\mu}$
(58b)

Once again we use the substitution $x^{4}=\mu/t$ to obtain

$I_{yy}^{\phantom{{\dagger}}}(\mu,\Lambda)=-\frac{\mu^{3/4}}{2\pi^{2}}\int\limits_{1}^{\mu/\Lambda^{4}}dt\,\,t^{-7/4}(1-t)^{1/2}.$
(59)

The substitution $t=1-l$ puts this in the form of an incomplete Beta function

$\displaystyle I_{yy}^{\phantom{{\dagger}}}(\mu,\Lambda)=$
$\displaystyle\frac{\mu^{3/4}}{2\pi^{2}}\int\limits_{0}^{1-\mu/\Lambda^{4}}dl\,\,l^{3/2-1}(1-l)^{-3/4-1}$
(60a) $\displaystyle=$
$\displaystyle\frac{\mu^{3/4}}{2\pi^{2}}B_{1-\mu/\Lambda^{4}}^{\phantom{{\dagger}}}(3/2,-3/4)$
(60b)

We can now expand around $\mu=0$

$\displaystyle I_{yy}(\mu,\Lambda)\approx\,\,\Lambda^{3}\bigg{(}$
$\displaystyle 0.07-0.18\,(\mu/\Lambda^{4})^{3/4}+0.10\,\mu/\Lambda^{4}$
$\displaystyle+O\left(\left(\mu/\Lambda\right)^{7/4}\right)\bigg{)}$ (61a)

We examine $\mu<0$. The contours for this case are shown in Fig 6(b). The
cutoff is implemented through the $k_{y}^{\phantom{{\dagger}}}$ integral that
has limits $\pm\Lambda^{2}$. The delta function identity to use in this case
is

$\displaystyle\delta\bigg{(}k_{x}^{4}$ $\displaystyle-
k_{y}^{2}-\mu\bigg{)}=\Theta\left(\Lambda^{2}-\sqrt{x^{4}+|\mu|}\right)$
$\displaystyle\times\left[\frac{\delta\left(k_{y}^{\phantom{{\dagger}}}-\sqrt{k_{x}^{4}+|\mu|}\right)}{|-2k_{y}^{\phantom{{\dagger}}}|}+\frac{\delta\left(k_{y}^{\phantom{{\dagger}}}+\sqrt{k_{x}^{4}+|\mu|}\right)}{|-2k_{y}^{\phantom{{\dagger}}}|}\right]$
(62)

Although the roots are always real (since $\mu<0$), we need the step function
to ensure that they fall within the integral domain, i.e
$\tilde{k}_{y}^{\pm}=\pm\sqrt{k_{x}^{4}+|\mu|}\in(-\Lambda^{2},\Lambda^{2})$.
This leads to a corresponding restriction on $k_{x}^{\phantom{{\dagger}}}$ as
$|k_{x}^{\phantom{{\dagger}}}|<(\Lambda^{4}-|\mu|)^{1/4}$. Let us treat
$I_{xx}^{\phantom{{\dagger}}}$ first. After some simplification, we get

$I_{xx}^{\phantom{{\dagger}}}(\mu,\Lambda)=\frac{8}{\pi^{2}}\int\limits_{0}^{(\Lambda^{4}-|\mu|)^{1/4}}dk_{x}^{\phantom{{\dagger}}}\frac{k_{x}^{6}}{\sqrt{k_{x}^{4}+|\mu|}}.$
(63)

By a series of substitutions, first $t=|\mu|/(k_{x}^{4}+|\mu|)$ and then
$l=1-t$ we get

$\displaystyle I_{xx}^{\phantom{{\dagger}}}(\mu,\Lambda)=$
$\displaystyle\frac{2|\mu|^{5/4}}{\pi^{2}}\int\limits_{0}^{1-|\mu|/\Lambda^{4}}dl\,\,l^{7/4-1}(1-l)^{-5/4-1}$
(64a) $\displaystyle=$
$\displaystyle\frac{2|\mu|^{5/4}}{\pi^{2}}B_{1-|\mu|/\Lambda^{4}}(7/4,-5/4).$
(64b)

Series expanding this, we obtain to leading order in the $\mu\rightarrow 0$
limit

$\displaystyle
I_{xx}^{\phantom{{\dagger}}}(\mu,\Lambda)\approx\Lambda^{5}\bigg{(}$
$\displaystyle
0.16-0.61\,|\mu|/\Lambda^{4}+0.41\,\left(|\mu|/\Lambda^{4}\right)^{5/4}$
$\displaystyle+O\left(\left(|\mu|/\Lambda^{4}\right)^{2}\right)\bigg{)}.$ (65)

Lastly, we compute $I_{yy}^{\phantom{{\dagger}}}$. Here again we first make
the substitution $t=|\mu|/(k_{x}^{4}+|\mu|)$ followed by the substitution
$l=1-t$ to get

$\displaystyle I_{yy}^{\phantom{{\dagger}}}(\mu,\Lambda)=$
$\displaystyle\frac{|\mu|^{3/4}}{2\pi^{2}}\int\limits_{0}^{1-|\mu|/\Lambda^{4}}dl\,\,l^{1/4-1}(1-l)^{1/4-1}$
(66a) $\displaystyle=$
$\displaystyle\frac{|\mu|^{3/4}}{2\pi^{2}}B_{1-|\mu|/\Lambda^{4}}(-1/4,-1/4).$
(66b)

Expanding this to leading order in $|\mu|$ we get

$\displaystyle I_{yy}(\mu,\Lambda)\approx\,\,\Lambda^{3}\bigg{(}$
$\displaystyle 0.07+0.25\,(\mu/\Lambda^{4})^{3/4}-0.15\,\mu/\Lambda^{4}$
$\displaystyle+O\left(\left(\mu/\Lambda\right)^{7/4}\right)\bigg{)}$ (67)

We see that integral is not symmetric between $\mu<0$ and $\mu>0$ although
both cases have the same $\mu\rightarrow 0$ limit. Also
$I_{xx}^{\phantom{{\dagger}}}(\mu,\Lambda)$ and
$I_{yy}^{\phantom{{\dagger}}}(\mu,\Lambda)$ are not equal. In fact we have in
the $\mu\rightarrow 0$ limit

$\frac{I_{xx}^{\phantom{{\dagger}}}(\mu,\Lambda)}{I_{yy}^{\phantom{{\dagger}}}(\mu,\Lambda)}\,\xrightarrow{\mu\rightarrow
0}\,\frac{12\Lambda^{2}}{5}.$ (68)

###  E.2 Monkey saddle

Figure 7: In evaluating the integral in Eq 9, the monkey saddle requires a
careful choice of UV cutoff since the box cutoff makes analytic treatment
difficult. As before, we color the regions with energy below and above the
singularity respectively blue and red. For any energy $\mu$, there are three
disconnected contours that are related by $2\pi/3$ rotations. This serves as
the motivation for using the hexagonal cutoff scheme for the integral.

We begin by noting the following

$\displaystyle\xi_{\mathbf{k}}^{\phantom{{\dagger}}}=:\xi(k_{x}^{\phantom{{\dagger}}},k_{y}^{\phantom{{\dagger}}})=k_{x}^{3}-3k_{x}^{\phantom{{\dagger}}}k_{y}^{2},$
(69) $\displaystyle\implies\,$
$\displaystyle\frac{\partial\xi_{\mathbf{k}}^{\phantom{{\dagger}}}}{\partial
k_{x}^{\phantom{{\dagger}}}}=3k_{x}^{2}-3k_{y}^{2},$ (70)
$\displaystyle\frac{\partial\xi_{\mathbf{k}}^{\phantom{{\dagger}}}}{\partial
k_{y}^{\phantom{{\dagger}}}}=-6k_{x}^{\phantom{{\dagger}}}k_{y}^{\phantom{{\dagger}}}.$
(71)

The constant energy contours at an energy $\mu$ are

$k_{y}^{\phantom{{\dagger}}}=\pm\sqrt{\frac{k_{x}^{3}-\mu}{3k_{x}^{\phantom{{\dagger}}}}}.$
(72)

For $\mu>0$, there are three disconnected pieces that are related by a
$2\pi/3$ rotation. These reside respectively in regions I, III and V in Fig
6(a). Thus for $\mu>0$, the integral can be restricted to these regions. In
region I, a cutoff $\Lambda$ on $k_{x}^{\phantom{{\dagger}}}$ serves as a
cutoff on $k_{y}$ as well, so that the region of integration becomes bounded
(a triangular region to be precise). However in III and V, a cutoff on
$k_{x}^{\phantom{{\dagger}}}$ does not automatically restrict the
$k_{y}^{\phantom{{\dagger}}}$ integral. This is because the
$k_{x}^{\phantom{{\dagger}}}\rightarrow 0$ limit corresponds to
$k_{y}^{\phantom{{\dagger}}}\rightarrow\pm\infty$. Now we might attempt to
additionally put a hard cutoff on $k_{y}^{\phantom{{\dagger}}}$, say
$\Lambda^{\prime}$. But this will also make the range of
$k_{x}^{\phantom{{\dagger}}}$ integration restricted and dependant on $\mu$
and $\Lambda^{\prime}$ (since every $k_{x}^{\phantom{{\dagger}}}$ will not
give a valid $k_{y}^{\phantom{{\dagger}}}$ on the $\mu$-contour). More
precisely, the lower and upper limits of the $k_{x}^{\phantom{{\dagger}}}$
integral will be the solutions of

$k_{x}^{3}-3k_{x}\Lambda^{\prime 2}=\mu.$ (73)

This is obviously not easy to work with. We therefore use a different cutoff
procedure, with a hexagonal region as depicted in Fig 7. The reason for this
is that, regions III and V can be ‘transformed’ into region I by simple
rotations of respectively $-2\pi/3$ and $2\pi/3$. Such rotations will not
change the dispersion, but only the derivatives
$\partial_{i}^{\phantom{{\dagger}}}\xi$ and
$\partial_{j}^{\phantom{{\dagger}}}\xi$, that change by linear combinations of
$\partial_{x}^{\phantom{{\dagger}}}\xi$ and
$\partial_{y}^{\phantom{{\dagger}}}\xi$ (By the chain rule since we have a
linear transformation of the coordinates).

Let us call the rotated variables $u$ and $v$. Let us
$\rho_{\theta}^{\phantom{{\dagger}}}$ denote the rotation matrix for angle
$\theta$. As mentioned above, we have
$\xi(u,v)=\xi(\rho_{\theta}^{\phantom{{\dagger}}}(k_{x}^{\phantom{{\dagger}}},k_{y}^{\phantom{{\dagger}}}))=\xi(k_{x}^{\phantom{{\dagger}}},k_{y}^{\phantom{{\dagger}}})$
but the derivatives become

$\displaystyle\frac{\partial\xi}{\partial k_{x}^{\phantom{{\dagger}}}}$
$\displaystyle\xrightarrow{\rho_{\pm
2\pi/3}^{\phantom{{\dagger}}}}-\frac{3}{2}(u^{2}-v^{2}\pm 2\sqrt{3}\,uv)$ (74)
$\displaystyle\frac{\partial\xi}{\partial k_{y}^{\phantom{{\dagger}}}}$
$\displaystyle\xrightarrow{\rho_{\pm
2\pi/3}^{\phantom{{\dagger}}}}\mp\frac{3\sqrt{3}}{2}(u^{2}-v^{2}\mp
2\,uv/\sqrt{3})$ (75)

Now the integrals over the three regions can be combined into a single
integral over region I using this scheme. We now compute
$(\partial_{i}^{\phantom{{\dagger}}}\xi)(\partial_{j}^{\phantom{{\dagger}}}\xi)$
summed over the three regions

$\displaystyle(\partial_{x}^{\phantom{{\dagger}}}\xi)^{2}|_{\textup{I}}^{\phantom{{\dagger}}}+(\partial_{x}^{\phantom{{\dagger}}}\xi)^{2}|_{\textup{II}}^{\phantom{{\dagger}}}+(\partial_{x}^{\phantom{{\dagger}}}\xi)^{2}|_{\textup{III}}^{\phantom{{\dagger}}}\rightarrow\frac{27}{2}(u^{2}+v^{2})^{2}$
(76)
$\displaystyle(\partial_{x}^{\phantom{{\dagger}}}\xi)(\partial_{y}^{\phantom{{\dagger}}}\xi)|_{\textup{I}}^{\phantom{{\dagger}}}+(\partial_{x}^{\phantom{{\dagger}}}\xi)(\partial_{y}^{\phantom{{\dagger}}}\xi)|_{\textup{II}}^{\phantom{{\dagger}}}+(\partial_{x}^{\phantom{{\dagger}}}\xi)(\partial_{y}^{\phantom{{\dagger}}}\xi)|_{\textup{III}}^{\phantom{{\dagger}}}\rightarrow
0$ (77)
$\displaystyle(\partial_{y}^{\phantom{{\dagger}}}\xi)^{2}|_{\textup{I}}^{\phantom{{\dagger}}}+(\partial_{y}^{\phantom{{\dagger}}}\xi)^{2}|_{\textup{II}}^{\phantom{{\dagger}}}+(\partial_{y}^{\phantom{{\dagger}}}\xi)^{2}|_{\textup{III}}^{\phantom{{\dagger}}}\rightarrow\frac{27}{2}(u^{2}+v^{2})^{2}$
(78)

We have to thus evaluate the following integral

$I(\mu,\Lambda)=\frac{27}{2}\int\limits_{0}^{\Lambda}\frac{du}{(2\pi)}\int\limits_{-\sqrt{3}x}^{\sqrt{3}x}\frac{dv}{(2\pi)}(u^{2}+v^{2})^{2}\delta(u^{3}-3uv^{2}-\mu)$
(79)

We proceed as in the previous section, first by using a delta function
identity to calculate the $v$ integral then use appropriate substitution to
recast it in terms of the incomplete Beta functions. We obtain

$\displaystyle I(\mu,\Lambda)=$
$\displaystyle\frac{2}{\sqrt{3}\pi^{2}}\,\frac{\mu}{\Lambda^{3}}\,B_{1-\mu/\Lambda^{3}}^{\phantom{{\dagger}}}\left(\frac{1}{2},-1\right)$
$\displaystyle-\frac{1}{\sqrt{3}\pi^{2}}\,\frac{\mu}{\Lambda^{3}}\,B_{1-\mu/\Lambda^{3}}\left(\frac{1}{2},0\right)$
$\displaystyle+\frac{1}{4\sqrt{3}\pi^{2}}\,\frac{\mu}{\Lambda^{3}}\sqrt{1-\frac{\mu}{\Lambda^{3}}}.$
(80)

The series expansion around $\mu=0$ reads

$\displaystyle I(\mu,\Lambda)\approx$
$\displaystyle\,0.12-0.04\,(\mu/\Lambda^{3})-0.02\,(\mu/\Lambda^{3})^{2}$
$\displaystyle+O(\,(\mu/\Lambda^{3})^{3}).$ (81)

It is easy to see by a simple change of variables that this derivation and
final result hold for $\mu<0$ as well.

## References

  * Lifshitz (1960) I. Lifshitz, Sov. Phys. JETP 11, 1130 (1960).
  * Van Hove (1953) L. Van Hove, Phys. Rev. 89, 1189 (1953).
  * Abrikosov (2017) A. A. Abrikosov,  _Fundamentals of the Theory of Metals_ (Courier Dover Publications, 2017).
  * Chandrasekaran  _et al._ (2020) A. Chandrasekaran, A. Shtyk, J. J. Betouras, and C. Chamon, Phys. Rev. Research 2, 013355 (2020).
  * Efremov  _et al._ (2019) D. V. Efremov, A. Shtyk, A. W. Rost, C. Chamon, A. P. Mackenzie, and J. J. Betouras, Phys. Rev. Lett. 123, 207202 (2019).
  * Shtyk  _et al._ (2017) A. Shtyk, G. Goldstein, and C. Chamon, Phys. Rev. B 95, 035137 (2017).
  * Yuan and Fu (2020) N. F. Q. Yuan and L. Fu, Phys. Rev. B 101, 125120 (2020).
  * Yuan  _et al._ (2019) N. F. Yuan, H. Isobe, and L. Fu, Nature communications 10, 1 (2019).
  * Kerelsky  _et al._ (2019) A. Kerelsky, L. J. McGilly, D. M. Kennes, L. Xian, M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, _et al._ , Nature 572, 95 (2019).
  * Aoki  _et al._ (2016) D. Aoki, G. Seyfarth, A. Pourret, A. Gourgout, A. McCollam, J. A. N. Bruin, Y. Krupko, and I. Sheikin, Phys. Rev. Lett. 116, 037202 (2016).
  * Barber  _et al._ (2019) 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, and I. I. Mazin, Phys. Rev. B 100, 245139 (2019).
  * Benhabib  _et al._ (2015) 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, and A. Forget, Phys. Rev. Lett. 114, 147001 (2015).
  * Khan and Johnson (2014) S. N. Khan and D. D. Johnson, Phys. Rev. Lett. 112, 156401 (2014).
  * Coldea  _et al._ (2019) 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 4, 1 (2019).
  * Okamoto  _et al._ (2010) Y. Okamoto, A. Nishio, and Z. Hiroi, Phys. Rev. B 81, 121102 (2010).
  * Sherkunov and Betouras (2018) Y. Sherkunov and J. J. Betouras, Phys. Rev. B 98, 205151 (2018).
  * Slizovskiy  _et al._ (2015) S. Slizovskiy, A. V. Chubukov, and J. J. Betouras, Phys. Rev. Lett. 114, 066403 (2015).
  * Stewart (1982) I. Stewart, Reports on Progress in Physics 45, 185 (1982).
  * Yelland  _et al._ (2011) E. Yelland, J. Barraclough, W. Wang, K. Kamenev, and A. Huxley, Nature physics 7, 890 (2011).
  * Classen  _et al._ (2020) L. Classen, A. V. Chubukov, C. Honerkamp, and M. M. Scherer, Phys. Rev. B 102, 125141 (2020).
  * Lin and Nandkishore (2020) Y.-P. Lin and R. M. Nandkishore, Phys. Rev. B 102, 245122 (2020).
  * Bi and Fu (2021) Z. Bi and L. Fu, Nature communications 12, 1 (2021).
  * Guerci  _et al._ (2021) D. Guerci, P. Simon, and C. Mora, “Higher-order Van hove singularity in magic-angle twisted trilayer graphene,” (2021), arXiv:2106.14911 [cond-mat.mes-hall] .
  * Oriekhov  _et al._ (2021) D. O. Oriekhov, V. P. Gusynin, and V. M. Loktev, Phys. Rev. B 103, 195104 (2021).
  * Rosenzweig  _et al._ (2020) P. Rosenzweig, H. Karakachian, D. Marchenko, K. Küster, and U. Starke, Phys. Rev. Lett. 125, 176403 (2020).
  * Zhou  _et al._ (2021) H. Zhou, Y. Saito, L. Cohen, W. Huynh, C. L. Patterson, F. Yang, T. Taniguchi, K. Watanabe, and A. F. Young, “Isospin magnetism and spin-triplet superconductivity in bernal bilayer graphene,” (2021), arXiv:2110.11317 [cond-mat.mes-hall] .
  * Castrigiano and Hayes (2019) D. P. Castrigiano and S. A. Hayes,  _Catastrophe theory_ (CRC Press, 2019).
  * Poston and Stewart (2014) T. Poston and I. Stewart,  _Catastrophe theory and its applications_ (Courier Corporation, 2014).
  * Bruce  _et al._ (1992) J. W. Bruce, J. W. Bruce, and P. Giblin,  _Curves and Singularities: a geometrical introduction to singularity theory_ (Cambridge university press, 1992).
  * Keldysh  _et al._ (1965) L. V. Keldysh _et al._ , Sov. Phys. JETP 20, 1018 (1965).
  * Kamenev and Andreev (1999) A. Kamenev and A. Andreev, Phys. Rev. B 60, 2218 (1999).
  * Efetov (1999) K. Efetov,  _Supersymmetry in disorder and chaos_ (Cambridge university press, 1999).
  * Mézard  _et al._ (1987) M. Mézard, G. Parisi, and M. A. Virasoro,  _Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications_ , Vol. 9 (World Scientific Publishing Company, 1987).
  * Bruus and Flensberg (2004) H. Bruus and K. Flensberg,  _Many-body quantum theory in condensed matter physics: an introduction_ (OUP Oxford, 2004).
  * Mahan (2000) G. D. Mahan,  _Many-Particle Physics_ (Springer Science & Business Media, 2000, 2000).
  * Haldane (1988) F. D. M. Haldane, Phys. Rev. Lett. 61, 2015 (1988).