# Analytic density of states of two-dimensional Chern insulator

Vera Uzunova  Institute of Theoretical Physics, Faculty of Physics, University
of Warsaw, ul. Pasteura 5, 02-093 Warszawa, Poland  Institute of Physics of
the National Academy of Sciences of Ukraine, 46 Nauky Ave., 03039 Kyiv,
Ukraine  Krzysztof Byczuk  Institute of Theoretical Physics, Faculty of
Physics, University of Warsaw, ul. Pasteura 5, 02-093 Warszawa, Poland

(August 28, 2024)

###### Abstract

We present analytic expressions for the density of states and its consistent
derivation for the two-dimensional Qi-Wu-Zhang (QWZ) Hamiltonian, a generic
model for the Chern topological insulators of class A. This density of states
is expressed in terms of elliptical integrals. We discuss and plot special
cases of the dispersion relations and the corresponding densities of states.
Spectral moments are also presented. The exact formulae ought to be useful in
determining physical properties of the non-interacting Chern insulators and
within the dynamical mean-field theory for interacting fermions with the QWZ
Hamiltonian in the non-interacting limit.

##  I Introduction

A topological insulator (TI) is a common name for the novel class of systems
with non-trivial topological properties Kane ; Zhang . Historically, the first
example of TI was a two-dimensional electron gas in a strong magnetic field
where the integer quantum Hall effect was observed Laughlin . After predicting
and later discovering of other examples of TIs nature-mat  the subject becomes
in a main stream of condensed matter physics cond-mat , of cold atoms in
optical lattices cold , of photonics photonics , and even of electric
engineering PT .

One of a possible path to investigate TIs is to study tight binding models
defined on particular lattices. Among various interesting examples are either
the Su-Schrieffer-Heeger model SSH  and the Rice-Mele model Rice  in one-
dimension or the Haldane model Haldane  and the Qi-Wu-Zhang (QWZ) model in two
dimensions QWZ . The latter one is defined on a square lattice and the
corresponding two dimensional Brillouin zone whereas the former one is
formulated on a hexagonal lattice.

In particular, the QWZ model is a well-known system in studying physics of
fermions such as bulk and edge properties, different topological states,
thermodynamics, transport, and many others Short . This model is also used as
a non-interacting part of the many-body Hamiltonian where the two-particle
interaction is included. Its further generalization to arbitrary dimensions
and even to the limit of infinite dimension proved that topological insulators
are possible in interacting systems as well Potthoff .

In spite of such a broad interest in the QWZ model its density of states (DOS)
is not yet determined analytically. Although the DOS by itself is not
sufficient to provide topological classification of a system, it is a basic
and very important quantity necessary to investigate thermodynamics,
thermodynamic phases, response of the system on different probes, and many
other quantities. To fill in this gap, in this article we derive analytical
formulae of the DOS in terms of elliptic integrals. We discuss details of the
derivation and basic properties of the DOS when the relevant model parameter
is varied.

The DOS, denoted here by $\rho(\Omega)$, counts the number of states in a
vicinity of a particular value of energy $\Omega$, i.e.
$dN=\rho(\Omega)d\Omega$. It can be obtained from the Green’s function.
Analytic derivation of the DOS, even for two-dimensional (2D) systems, is
typically a challange and thus there are very few known analytical results.
One of the first example was obtained in 1953 by Hobson and Nierenberg Hobson
. It is an analytical expression for the DOS of graphene with the nearest-
neighbor hopping and represented by the complete elliptical integrals. The
consequent derivation of this result is presented in Ovchynnikov . The DOS for
some others 2D latices are also obtained analytically Kogan . In particular,
it can be done for square, triangular, honeycomb, Kagome, and Lieb lattices.
In this paper the DOS is analytically derived for the QWZ Hamiltonian modeling
a Chern insulator on the square lattice.

Our presentation is organized as follows: In section II we define the QWZ
model and discuss the dispersion relations, in Section III we introduce the
DOS and present its analytic derivation for the QWZ Hamiltonian, Section IV is
devoted to the discussion of the DOS in different parameter regimes, in
Section V we show some additional features in the DOS, the spectral moments
are discussed in Section VI, and we close our presentation with Section VII,
where we offer our conclusions and outlooks. In Appendix A we provide
mathematical definitions of the elliptic integrals and in the Appendix B we
give selected details on calculating the spectral moments.

##  II QWZ model Hamiltonian in two dimensions

A generic form of the two-band Hamiltonian for a 2D noninteracting system in
the momentum space can be written as

$\hat{H}=\sum_{\bf k}\hat{H}_{\bf k}=\sum_{\bf k}\mathbf{h}({\bf
k})\cdot\hat{\boldsymbol{\sigma}},$ (1)

where $\mathbf{k}=(k_{x},k_{y})$ ($-\pi/a_{L}<$ $k_{x},$ $k_{y}$
$\leq\pi/a_{L}$) is a 2D wave vector in the first Brillouin zone corresponding
to the 2D square lattice with the lattice constant $a_{L}$, ${\bf h}({\bf k})$
is a vector with three components being given functions of $\bf k$, and
$\hat{\boldsymbol{\sigma}}$ is the vector with components represented by the
three Pauli matrices $\hat{\sigma}_{x}$, $\hat{\sigma}_{y}$, and
$\hat{\sigma}_{z}$.

The Hamiltonian (1) describes a two level system, corresponding to either the
two orbitals or the spin $1/2$ degrees of freedom. The vector ${\bf h}({\bf
k})$ is interpreted as a Zeeman like magnetic field with some (perhaps non-
trivial) dependence on the wave vector $\bf k$. This model breaks the time
reversal symmetry and belongs to the class A in the ten fold way
classification scheme ten-way . The Hamiltonian (1) can be easily
diagonalized, giving a two band energy spectrum
$\epsilon_{\pm}(\mathbf{k})=\pm h({\bf k})$, where $h({\bf k})=|{\bf h}({\bf
k})|$ is the length of the vector $\mathbf{h}(\mathbf{k})$.

In the following we consider a particular parametrization where the length of
${\bf h}({\bf k})$ is given by $h({\bf
k})^{2}=m^{2}+2t^{2}+2t^{2}\cos(k_{x}a_{L})\cos(k_{y}a_{L})+2mt[\cos(k_{x}a_{L})+\cos(k_{y}a_{L})]$.
It corresponds to the following representation of the vector $\mathbf{h}({\bf
k})$:

$\displaystyle h_{x}({\bf k})=t\sin(k_{x}a_{L}),$ $\displaystyle h_{y}({\bf
k})=t\sin(k_{y}a_{L}),$ $\displaystyle h_{z}({\bf
k})=m+t\cos(k_{x}a_{L})+t\cos(k_{y}a_{L}),$ (2)

where $t$ is the hopping amplitude. In the momentum space this vector
$\mathbf{h}(\mathbf{k})$ has a Skyrmion configuration for $0<|m|/t<2$ Short .
In other words, the system is a topological insulator with the finite Chern
number $\pm 1$ and conducting surface states at the half-filling. Hamiltonian
(1) can be interpreted as a tight-binding model of a magnetic semiconductor
with the Rashba-type spin-orbit interaction and a uniform magnetization along
the $z$-axis Short . In the following we take $t=1$, which sets the energy
unit. We also use $a_{L}=1$ for the length unit.

For $m=0$, $0.5$, $1$, $2$, and $2.5$ the corresponding eigenvalues (energy
bands) of the Hamiltonian (1) are plotted in Figs. 1-5, respectively. We see
that at $m=0$ (Fig. 1) and $m\pm 2$ (Fig. 4) the band gap is closed at
$\Gamma$, $X$ ($Y$) and $M$ special points in the square Brillouin zone,
respectively, and characteristic Dirac cone(s) is (are) formed. For $m\neq 0$,
$\pm 2$ the bands are separated by the gap as seen in Figs. 2, 3, and 5.

Figure 1:  Dispersion relations $\epsilon_{\pm}(\mathbf{k})$ of the QWZ model
at $m=0$. The gap is closed at $X=(\pm\pi,0)$ and $Y=(0,\pm\pi)$ points.

Figure 2: Dispersion relations $\epsilon_{\pm}(\mathbf{k})$ of the QWZ model
at: (a) $m=0.5$, and (b) $m=-0.5$.

Figure 3:  Dispersion relations $\epsilon_{\pm}(\mathbf{k})$ of the QWZ model
at: (a) $m=1$, and (b) $m=-1$.

Figure 4:  Dispersion relations $\epsilon_{\pm}(\mathbf{k})$ of the QWZ model
at: (a) $m=2$, and (b) $m=-2$. For $m=2$ the gap is closed at
$M=(\pm\pi,\pm\pi)$ and $(\pm\pi,\mp\pi)$ points and for $m=-2$ it is closed
at $\Gamma=(0,0)$ point.

Figure 5:  Dispersion relations $\epsilon_{\pm}(\mathbf{k})$ of the QWZ model
at: (a) $m=2.5$, and (b) $m=-2.5$.

##  III Analytic derivation of DOS

###  III.1 DOS and its symmetry

The band resolved DOS per lattice site for the lower
$\epsilon_{-}(\mathbf{k})=-h({\bf k})$ (upper
$\epsilon_{+}(\mathbf{k})=+h({\bf k})$) band is defined as usual

$\rho_{\pm}(\Omega)=\frac{1}{N_{L}}\sum_{\bf
k}\delta(\Omega-\epsilon_{\pm}(\mathbf{k})),$ (3)

where $\delta(x)$ is the Dirac distribution function and $N_{L}$ is the number
of lattice sites. Due to the symmetry of the Dirac function one can easily
show that

$\rho_{+}(\Omega)=\rho_{-}(-\Omega).$ (4)

The total DOS per lattice site is defined as a sum of them, i.e.,

$\rho(\Omega)=\rho_{+}(\Omega)+\rho_{-}(\Omega).$ (5)

Since the lower and upper bands do not overlap footnote  the band resolved DOS
can be directly extracted form the total DOS, i.e.,

$\rho_{\pm}(\Omega)=\theta(\pm\Omega)\rho(\Omega),$ (6)

where $\theta(x)$ is a step function. Therefore, in the following we focus on
deriving a formula for the total DOS.

###  III.2 Retarded Green’s function and total DOS

Equivalently, the total DOS per lattice site is determined as a trace in band
indices ($\pm$) of the imaginary part of the retarded Green’s function

$\rho(\Omega)=-\frac{1}{\pi}\Im[\operatorname{Tr}\mathbf{G}(\Omega)].$ (7)

The diagonal matrix elements of the Green’s function $\mathbf{G}(\Omega)$ are

$G_{\pm}(\Omega)=\frac{1}{N_{L}}\sum_{\mathbf{k}}\frac{1}{\Omega-\epsilon_{\pm}(\mathbf{k})+\imath
0^{+}},$ (8)

and for a while the small imaginary part $\imath 0^{+}$ will not be written
explicitly. Summation over k can be replaced by the continuous integral in the
first Brillouin zone
$\sum_{\mathbf{k}}\rightarrow\frac{L^{2}}{(2\pi)^{2}}\int_{BZ}dk_{x}dk_{y}$,
where $L$ is the length of the system and $L^{2}=N_{L}a_{L}^{2}$, with
$a_{L}=1$. Then the trace of the Green’s function has the form

$\operatorname{Tr}\mathbf{G}(\Omega)=\frac{\Omega}{2\pi^{2}}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}\frac{dk_{x}dk_{y}}{\Omega^{2}-h({\bf
k})^{2}}.$ (9)

Integrations with respect to $k_{x}$ and $k_{y}$ are symmetric. We first
preform the integration of the function $(\Omega^{2}-h(k_{x},k_{y})^{2})^{-1}$
with respect to $k_{x}$. It is convenient to denote

$\displaystyle a=\Omega^{2}-(m^{2}+2+2m\cos k_{y}),$ $\displaystyle
b=-2(m+\cos k_{y}),$ (10)

and then to use the identity (2.558.4) in Gradshteyn

$\displaystyle\int\frac{dk_{x}}{a+b\cos
k_{x}}=\frac{2\pi}{\sqrt{a^{2}-b^{2}}}\arctan\left[\frac{a-b}{\sqrt{a^{2}-b^{2}}}\tan\left(\frac{k_{x}}{2}\right)\right],$
(11)

where $a^{2}-b^{2}>0$. For the case $a^{2}-b^{2}\leq 0$ the result of the
integration (11) can be formally rewritten as

$\int_{-\pi}^{\pi}\frac{dk_{x}}{\Omega^{2}-h^{2}}=\frac{-2\pi\imath}{\sqrt{b^{2}-a^{2}}}\operatorname{csgn}\left[\frac{\imath(b-a)}{\sqrt{b^{2}-a^{2}}}\right],$
(12)

where it is taken into account that $\arctan(\pm\tan\pi/2)=\pm\pi/2$. The
complex signum function, abbreviated as $\operatorname{csgn}$, is equal to the
sign (function $\operatorname{sgn}$) of the real part of the argument, and
sign of the imaginary part if the real part is zero. We analytically continue
this result by taking $\Omega\rightarrow\Omega+\imath 0^{+}$, where we add
infinitesimally small imaginary part as it is required in the retarded Green’s
function. Then $a=(\Omega+\imath 0^{+})^{2}+...$ gets small imaginary part
$2\imath\Omega 0^{+}$ defining the value of $\operatorname{csgn}$. Taking the
limit of infinitisimal $0^{+}$ we arrive at the following expression for the
trace of the Green’s function

$\operatorname{Tr}\mathbf{G}(\Omega)=\frac{2\Omega}{\pi}\int_{0}^{\pi}dk_{y}\begin{cases}\frac{\operatorname{sgn}{a}}{\sqrt{a^{2}-b^{2}}}&\text{for}\;\;a^{2}>b^{2},\\\
\imath\frac{\operatorname{sgn}{\Omega}}{\sqrt{b^{2}-a^{2}}}&\text{for}\;\;a^{2}\leq
b^{2},\end{cases}$ (13)

where we used the fact that the integrand is an even function of $k_{y}$. Here
the values $a$ and $b$ are the functions of the variables $k_{y}$ and $\Omega$
and the parameter $m$ as determined by Eqs. (III.2).

The DOS is proportional to the imaginary part of $\operatorname{Tr}G(\Omega)$,
therefore it is nonzero only in the region of $\Omega$ where

$a^{2}-b^{2}\leq 0.$ (14)

Outside the region (14) $\operatorname{Tr}G$ is real and the DOS is equal to
zero, i.e., $\rho(\Omega)=0$.

The next step in evaluation of Eq. (13) is to perform integration over
$k_{y}$. It is convenient to replace $y=\cos k_{y}$ and
$dk_{y}=-dy(1-y^{2})^{-1/2}$. The boundaries of the integration are $-1\leq
y\leq 1$, where changing the order of the integration boundaries results an
additional minus sign. The integration over $y$ is performed differently
depending on the value of the parameter $m$. Therefore, in what follows we are
considering three cases with $|m|=1$, $|m|>1$, and $|m|<1$, separately.

####  III.2.1 Case $|m|=1$

Calculations are simpler in the special case of $|m|=1$. The function in the
denominator of the integral (13) can be written explicitly as
$a^{2}-b^{2}=(\Omega^{2}-1)(\Omega^{2}-5-4my)$. It is linear in $y$ and
changes the sign only once at the point $y=y_{0}\operatorname{sgn}{m}$, where
$y_{0}=(\Omega^{2}-5)/4$. Solving the condition (14) with respect to $\Omega$
and using the fact that $|y|\leq 1$ we obtain that DOS is nonzero only for
$1\leq|\Omega|\leq 3$.

The condition (14) considered with respect to $y$ determines the boundaries of
integration in Eq. (13). For $m=1$ it is satisfied for $y_{0}\leq y\leq
1$.Then the imaginary part of Eq. (13) in terms of variable $y$ gives the DOS
in the implicit form

$\displaystyle\rho(\Omega)=\frac{|\Omega|}{\pi^{2}\sqrt{\Omega^{2}-1}}\int_{y_{0}}^{1}\frac{dy}{\sqrt{(1-y^{2})(y-y_{0})}},$
(15)

where the definition (7) is used. The similar expression will appears for
$m=-1$: the boundaries of integration are $-1\leq y\leq-y_{0}$ and denominator
in the integrand is $\sqrt{(1-y^{2})(-y-y_{0})}$. The corresponding DOS can be
transformed into the result Eq. (15) by replacing the integration variable
$y\rightarrow-y$. So Eq. (15) is valid for both cases $m=\pm 1$.

The integral over $y$ can be calculated in terms of the complete elliptical
integral of the first kind $K(x)$ using the identity (3.131.5) from Gradshteyn

$\displaystyle\int_{u_{2}}^{u_{3}}\frac{dy}{\sqrt{(u_{3}-y)(y-u_{2})(y-u_{1})}}$
$\displaystyle=\frac{2}{\sqrt{u_{3}-u_{1}}}K\left[\sqrt{\frac{(u_{3}-u_{2})}{(u_{3}-u_{1})}}\right],$
(16)

where $u_{3}>u_{2}>u_{1}$. In Eq. (15) these parameters are $u_{1}=-1$,
$u_{2}=y_{0}$, $u_{3}=1$ and the result of integration is equal to
$\sqrt{2}K\left[\sqrt{(1-y_{0})/2}\right]$. Then the DOS for $|m|=1$ is

$\rho_{|m|=1}(\Omega)=\begin{cases}\frac{\sqrt{2}}{\pi^{2}}\frac{|\Omega|}{\sqrt{\Omega^{2}-1}}K\left[\sqrt{\frac{9-\Omega^{2}}{8}}\right],\;\;\;\;\text{if}\;\;1\leq\Omega^{2}\leq
9,\\\
0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\text{otherwise}.\end{cases}$
(17)

For convenience, in the Appendix A we provide the definitions of the elliptic
integrals.

####  III.2.2 Case $|m|>1$

To calculate the DOS for all other values of the parameter $m$ we need first
to analyze the function in denominator of Eq. (13). For $|m|\neq 1$ we can
write $a^{2}-b^{2}=4(m^{2}-1)(y-y_{1})(y-y_{2})$. Here the left and right
zeros of the function are denoted as $y_{1}=\min(\tilde{y}_{1},\tilde{y}_{2})$
and $y_{2}=\max(\tilde{y}_{1},\tilde{y}_{2})$, respectively, where

$\displaystyle\tilde{y}_{1}=\frac{\Omega^{2}-1-(m+1)^{2}}{2(m+1)},$
$\displaystyle\tilde{y}_{2}=\frac{\Omega^{2}-1-(m-1)^{2}}{2(m-1)}.$ (18)

So we have to set $y_{1}=\tilde{y}_{1}$, $y_{2}=\tilde{y}_{2}$ if the the
following condition is met

$\displaystyle\frac{\Omega^{2}+m^{2}-2}{m^{2}-1}>0,$ (19)

and to set $y_{1}=\tilde{y}_{2}$, $y_{2}=\tilde{y}_{1}$, otherwise.

Lets consider the case $|m|>1$, then factor $m^{2}-1$ is positive. Rewriting
Eq. (13) in terms of variable $y$ and substituting it in Eq. (7), we obtain
the nonzero DOS in the implicit form

$\displaystyle\rho(\Omega)=\frac{|\Omega|}{\pi^{2}\sqrt{m^{2}-1}}\int\frac{dy}{\sqrt{(1-y^{2})(y_{1}-y)(y-y_{2})}}.$
(20)

Here the boundaries of integration are determined by Eq. (14), which reads
$y_{1}\geq y\geq y_{2}$, implying $|y|\leq 1$. There are two cases, in which
these conditions can be satisfied by the integration variable $y$:

$\displaystyle|y_{1}|\leq 1,|y_{2}|>1,$ $\displaystyle|y_{1}|>1,|y_{2}|\leq
1.$ (21)

First one leads to the integration region $y_{1}\leq y\leq 1$, and the second
one to $-1\leq y\leq y_{1}$. To perform integration we use the identity
(3.147.5) in Gradshteyn , namely

$\displaystyle\int_{u_{2}}^{u_{3}}\frac{dy}{\sqrt{(u_{4}-y)(u_{3}-y)(y-u_{2})(y-u_{1})}}$
$\displaystyle=\frac{2}{\sqrt{(u_{4}-u_{2})(u_{3}-u_{1})}}K\left[\sqrt{\frac{(u_{3}-u_{2})(u_{4}-u_{1})}{(u_{4}-u_{2})(u_{3}-u_{1})}}\right],$
(22)

where $u_{4}>u_{3}>u_{2}>u_{1}$. The first line of Eq. (III.2.2) corresponds
to $u_{1}=-1$, $u_{2}=y_{1}$, $u_{3}=1$ and $u_{4}=y_{2}$ and the second line
corresponds to $u_{1}=y_{1}$, $u_{2}=-1$, $u_{3}=y_{2}$ and $u_{4}=1$. The
results of the integration in both cases are identical

$\displaystyle\rho(\Omega)=\frac{\sqrt{2}}{\pi^{2}}\frac{|\Omega|}{\sqrt{m^{2}-1}}\frac{1}{\sqrt{y_{2}-y_{1}}}K\left[\sqrt{\frac{(1-y_{1})(1+y_{2})}{2(y_{2}-y_{1})}}\right].$
(23)

The last step of this evaluation is to express the obtained formula in terms
of $m$ and $\Omega$. As follows from Eq. (19) for all $\Omega^{2}>2-m^{2}$ we
can set $y_{1}=\tilde{y}_{1}$, $y_{2}=\tilde{y}_{2}$, where $\tilde{y}_{1}$,
$\tilde{y}_{2}$ are given by Eq. (III.2.2). Then the solution of the
inequalities (III.2.2) is given by the condition
$(|m|-2)^{2}<\Omega^{2}<(|m|+2)^{2}$. After making sure that
$2-m^{2}<(|m|-2)^{2}$ we substitute the same $y_{1}$, $y_{2}$ into Eq. (23) to
obtain

$\displaystyle\rho_{\text{I}}(\Omega)=\frac{\sqrt{2}}{\pi^{2}}\frac{|\Omega|}{\sqrt{\Omega^{2}+m^{2}-2}}\cdot$
$\displaystyle\cdot
K\left[\sqrt{\frac{[\Omega^{2}-(m-2)^{2}][(m+2)^{2}-\Omega^{2}]}{8(\Omega^{2}+m^{2}-2)}}\right],$
(24)

where we denote this result by $\rho_{\text{I}}(\Omega)$.

The final expression for the DOS for the case $|m|>1$ can be presented as

$\rho_{|m|>1}(\Omega)=\begin{cases}\rho_{\text{I}}(\Omega),\;\;\;\;\text{if}\;\;(|m|-2)^{2}\leq\Omega^{2}\leq(|m|+2)^{2},\\\
0,\;\;\;\;\;\;\;\;\;\;\;\text{otherwise}.\end{cases}$ (25)

Note, that this expression reduces to the case $m^{2}=1$. In this case $|m|=1$
Eq. (25) coincides exactly with Eq. (17).

####  III.2.3 Case $|m|<1$

We repeat all the reasoning of the previous subsection for the case $|m|<1$.
In this case the factor $|m|-1$ is negative, and the nonzero DOS obtained from
Eqs. (13) and (7) in terms of the variable $y$ is of the form

$\displaystyle\rho(\Omega)=\frac{|\Omega|}{\pi^{2}\sqrt{1-m^{2}}}\int\frac{dy}{\sqrt{(1-y^{2})(y-y_{1})(y-y_{2})}},$
(26)

where $|y|\leq 1$ by the definition. The boundaries of the integration are
determined by Eq. (14), that gives $y\leq y_{1}$ or $y\geq y_{2}$. There are
three different possibilities for $y$ to satisfy these conditions,

$\displaystyle|y_{1}|\leq 1,|y_{2}|>1,$ $\displaystyle|y_{1}|>1,|y_{2}|\leq
1,$ $\displaystyle|y_{1}|<1,|y_{2}|<1.$ (27)

We are considering each of them.

To perform integration in Eq. ( 26) we use the identities (3.147.3) and
(3.147.7) in Gradshteyn  respectively, that for the case of complete
elliptical integrals have the same right hand side, namely

$\displaystyle\int_{u_{1}}^{u_{2}}\frac{dy}{\sqrt{(u_{4}-y)(u_{3}-y)(u_{2}-y)(y-u_{1})}}=$
(28a)
$\displaystyle=\int_{u_{3}}^{u_{4}}\frac{dy}{\sqrt{(u_{4}-y)(y-u_{3})(y-u_{2})(y-u_{1})}}=$
(28b)
$\displaystyle=\frac{2}{\sqrt{(u_{4}-u_{2})(u_{3}-u_{1})}}K\left[\sqrt{\frac{(u_{4}-u_{3})(u_{2}-u_{1})}{(u_{4}-u_{2})(u_{3}-u_{1})}}\right],$

where $u_{4}>u_{3}>u_{2}>u_{1}$.

The first line of Eq. (III.2.3) gives the boundaries of integration as $-1\leq
y\leq y_{1}$ in the integral (26). The integration can be performed by using
Eq. (28a) with $u_{1}=-1$, $u_{2}=y_{1}$, $u_{3}=1$ and $u_{4}=y_{2}$. The
second line of Eq. (III.2.3) gives the boundaries of integration as $y_{2}\leq
y\leq 1$ in the integral (26). The integration can be performed by using Eq.
(28b) with $u_{1}=y_{1}$, $u_{2}=-1$, $u_{3}=y_{2}$ and $u_{4}=1$. The results
in these two cases are identical

$\displaystyle\rho(\Omega)=\frac{\sqrt{2}}{\pi^{2}}\frac{|\Omega|}{\sqrt{1-m^{2}}}\frac{1}{\sqrt{y_{2}-y_{1}}}K\left[\sqrt{\frac{(1+y_{1})(y_{2}-1)}{2(y_{2}-y_{1})}}\right].$
(29)

The third line of Eq. (III.2.3) splits integral (26) into the sum of two
integrals with boundaries $-1\leq y\leq y_{1}$ and $y_{2}\leq y\leq 1$. The
first one is performed by using Eq. (28a) and the second one is performed by
using Eq. (28b), where we put $u_{1}=-1$, $u_{2}=y_{1}$, $u_{3}=y_{2}$ and
$u_{4}=1$ for both integrals. The results of integrations are identical and
summation is reduces to multiplication by $2$. The final expression reads

$\displaystyle\rho(\Omega)=\frac{2\sqrt{2}}{\pi^{2}}\frac{|\Omega|}{\sqrt{1-m^{2}}}\frac{1}{\sqrt{(1-y_{1})(1+y_{2})}}\cdot$
$\displaystyle\cdot
K\left[\sqrt{\frac{(1+y_{1})(1-y_{2})}{(1-y_{1})(1+y_{2})}}\right].$ (30)

As the last step we are expressing the formulas (29) and (III.2.3) in terms of
$m$ and $\Omega$. Eq. (29) corresponds to the first two lines of Eq.
(III.2.3). As follows from Eq. (19) for all $\Omega^{2}>2-m^{2}$ we can set
$y_{1}=\tilde{y}_{2}$, $y_{2}=\tilde{y}_{1}$, where $\tilde{y}_{1}$,
$\tilde{y}_{2}$ are given by Eq. (III.2.2). Note that the first two lines of
Eq. (III.2.3) are equivalent to Eq. (III.2.2) and are symmetrical with respect
to permutation of $y_{1}$ and $y_{2}$. So we can repeat the same reasoning as
for the case $|m|>1$. The solution of this condition with respect to
$\Omega^{2}$ gives the same region $(|m|-2)^{2}<\Omega^{2}<(|m|+2)^{2}$ and
substitution of $y_{1}$ and $y_{2}$ into Eq. (29) results in the same DOS
$\rho_{\text{I}}(\Omega)$, given by Eq. (25).

To evaluate Eq. (III.2.3) corresponding to the third line of Eq. (III.2.3) we
notice that this line is symmetric with respect to permutation of $y_{1}$ and
$y_{2}$. Rewriting it in terms of $\Omega^{2}$ we obtain the region
$m^{2}>\Omega^{2}>(|m|-2)^{2}$ in which Eq. (III.2.3) is valid. This region is
split into two parts by Eq. (19). In particular, for
$2-m^{2}>\Omega^{2}>(|m|-2)^{2}$ we set $y_{1}=\tilde{y}_{2}$,
$y_{2}=\tilde{y}_{1}$, and substitute them to Eq. (III.2.3) to obtain

$\displaystyle\rho_{\text{II}}(\Omega)=\frac{8}{\pi^{2}}\frac{|\Omega|}{|\Omega^{2}+m^{2}|}\cdot$
$\displaystyle\cdot
K\left[\sqrt{\frac{[\Omega^{2}-(m-2)^{2}][\Omega^{2}-(m+2)^{2}]}{(\Omega^{2}-m^{2})^{2}}}\right],$
(31)

where we denote this result by $\rho_{\text{II}}(\Omega)$. In the region
$m^{2}>\Omega^{2}>2-m^{2}$ we set $y_{1}=\tilde{y}_{1}$,
$y_{2}=\tilde{y}_{2}$, and substitution to Eq. (III.2.3) gives

$\displaystyle\rho_{\text{III}}(\Omega)=\frac{8}{\pi^{2}}\frac{|\Omega|}{\sqrt{[\Omega^{2}-(m-2)^{2}][\Omega^{2}-(m+2)^{2}]}}\cdot$
$\displaystyle\cdot
K\left[\sqrt{\frac{(\Omega^{2}-m^{2})^{2}}{[\Omega^{2}-(m-2)^{2}][\Omega^{2}-(m+2)^{2}]}}\right],$
(32)

where we denote this result by $\rho_{\text{III}}(\Omega)$. The final
expression for the DOS for the case $|m|<1$ can be presented as

$\rho_{|m|<1}(\Omega)=\begin{cases}\rho_{\text{I}}(\Omega),\;\;\;\;\text{if}\;\;(|m|-2)^{2}<\Omega^{2}\leq(|m|+2)^{2},\\\
\rho_{\text{II}}(\Omega),\;\;\;\text{if}\;\;2-m^{2}\leq\Omega^{2}\leq(|m|-2)^{2},\\\
\rho_{\text{III}}(\Omega),\;\;\text{if}\;\;m^{2}\leq\Omega^{2}<2-m^{2},\\\
0,\;\;\;\;\;\;\;\;\;\;\;\text{otherwise}.\end{cases}$ (33)

Note, that the DOS for all values of parameter $m$ does not depend on the sign
of $m$ and, therefore is symmetric with respect to $m$ and $-m$ in contrast to
the dispersion relations.

##  IV Summary of analytical results and plots of total DOS

In this Section we present plots of the total DOS for different $m$,
corresponding to the dispersion relations shown in Section II, and we provide
consistent analyzes of them. The plots are given in Figs. 6-10, that appear in
the same order as the dispersion relations in Figs. 1-5. In these DOS plots
the ranges of the vertical axis are different but all DOS are normalized to
two, corresponding to two bands. We see that the shapes of the DOS for the QWZ
model are much richer and with more additional features as compared, for
example, to the hexagonal lattice with nearest neighbor hopping. Apart of the
symmetry $m$ and $-m$, mentioned above, the total DOS is symmetric with
respect to $\Omega$ and $-\Omega$. It is clearly visible in analytic formulae
of the DOS, Eqs. (17, 25, 33), where $\Omega$ is only present as either
$\Omega^{2}$ or $|\Omega|$.

Figure 6: Total density of states of the QWZ model at $m=0$.  Figure 7: Total
density of states of the QWZ model at $m=\pm 0.5$.

The typical shape of the DOS for $0<|m|<1$ is shown in Fig. 7. We can see that
each nonvanishing part of the plot has three sections, in each of the section
the function is described by $\rho_{\text{I}}$, $\rho_{\text{II}}$, or
$\rho_{\text{III}}$, cf Eq. (33). Specifically, in the plot in Fig. 7 for
$|m|=0.5$, these sections are separated by the points
$|\Omega|=\sqrt{1.75}\approx 1.32$ and $|\Omega|=1.5$. At the energies
$\pm\sqrt{2-m^{2}}$, separating two nearby sections, the DOS has infinite
peaks.

Another characteristic feature of the system with $0<|m|<1$ is the opening of
the band gap of the width $2|m|$. It can be seen in the corresponding
dispersion relations, cf. Fig. 2, and in the plots of DOS, where
$\rho(\Omega)=0$ in the range $-|m|<\Omega<|m|$. At the half-filling such
system is a topological insulator QWZ ; Short .

The gap is closed for $m=0$ as seen in Fig. 6, when the DOS has a pseudogap at
$\Omega=0$ (DOS vanishes at a single point). It corresponds to formation of
the Dirac cones at $X=(\pm\pi,0)$ and $Y=(0,\pm\pi)$ points in the Brillouin
zone, that is easy to see in the plots of the dispersion relations, cf. Fig.
1. At the half-filling such system is a semi-metal.

Figure 8: Total density of states of the QWZ model at $m=\pm 1$.

The special case $|m|=1$, given by Eq. (17), is shown in Fig. 8 and the
corresponding dispersion relations are in Fig. 3. Flat parts in the dispersion
relations, i.e. lines along which ${\bf\nabla}_{\bf k}\epsilon_{\pm}({\bf
k})=0$, give rise to the appearance of sharp peaks in the DOS. Despite of the
fact that these flat parts are different for $m=1$ and $m=-1$, the shape of
the DOS is the same. In this $|m|=1$ case the system is a topological
insulator at the half-filling QWZ ; Short .

Figure 9: Total density of states of the QWZ model at $m=\pm 2$. Figure 10:
Total density of states of the QWZ model at $m=\pm 2.5$.

The typical shape of the DOS for $|m|>1$ (without $|m|=2$) is shown in Fig.
10. The nonzero part of the DOS function is described by the single elliptic
integral, cf. Eq. (25). It has two symmetrical infinite peaks at energies
$\pm|m|$. The system has the band gap of the width $2||m|-2|$ in the
dispersion relation, shown in Fig. 5. In the plot of the DOS the gap
corresponds to $\rho(\Omega)=0$ in the range $-||m|-2|<\Omega<||m|-2|$. At the
half-filling the system with $|m|<2$ is a topological insulator and the system
with $|m|>2$ is a trivial insulator QWZ ; Short .

For $|m|=2$ the band gap closes and a pseudogap appears as it is shown in the
plot of the DOS in Fig. 9. Formation of the Dirac cones can be seen at the
corresponding dispersion relations in Fig. 4. For $m=2$ the gap is closed at
$M=(\pm\pi,\pm\pi)$ and $(\pm\pi,\mp\pi)$ points in the Brillouin zone and for
$m=-2$ it is closed at $\Gamma=(0,0)$ point. At the half filling such system
is a semi-metal.

##  V Subtle features seen in total DOS

In this Section, we discus additional subtle features and general trends that
can be observed for the DOS of the considered QWZ model.

###  V.1 Additional finite peaks

As discussed earlier, the total DOS is symmetric with respect to $m$ and $-m$
and as well as it is symmetric with respect to $\Omega$ and $-\Omega$. It has
two infinite peaks located at $\pm\Omega_{\infty}$, where the elliptical
integral $K(1)=\infty$, i.e.,

$\Omega_{\infty}=\begin{cases}\sqrt{2-m^{2}},\;\;\;\text{if}\;\;|m|<1,\\\
{|m|},\;\;\;\;\;\;\;\;\;\;\;\;\text{if}\;\;|m|\geq 1.\end{cases}$ (34)

We find, that for the values of $|m|$ slightly larger than $1$, two additional
finite peaks appear at the edges of a band gap given by $\pm\Omega_{L}$, for
example, it is shown in Fig. 11 for the case of $|m|=1.15$. At larger $|m|$
these peaks disappear.

Figure 11: Total density of states of the QWZ model at $m=\pm 1.15$.

###  V.2 Widths of the band gap and the bands

Let $\Delta$ be the width of the band gap. Its value for arbitrary $m$ (in
units with $t=1$) is given by the simple expression

$\Delta=2\begin{cases}{|m|},\;\;\;\;\;\;\;\;\;\;\;\text{if}\;\;|m|<1,\\\
{||m|-2|},\;\;\;\text{if}\;\;|m|\geq 1.\end{cases}$ (35)

Edges of the band gap are located at energies $\pm\Omega_{L}$, where
$\Omega_{L}=\Delta/2$.

On the other hand, the upper and lower bounds of the energy spectrum (the
dispersion relations $\epsilon_{\pm}({\bf k})$) are at $\pm\Omega_{R}$, where
$\Omega_{R}=|m|+2$. The width of each of the two bands $\epsilon_{\pm}({\bf
k})$ is defined as $W=\Omega_{R}-\Omega_{L}$, since $\Omega_{R}>\Omega_{L}$.
Its value is explicitly given by

$W=2\begin{cases}1,\;\;\;\;\;\;\;\;\;\;\;\text{if}\;\;|m|<1,\\\
|m|,\;\;\;\;\;\;\;\text{if}\;\;1\geq|m|\geq 2,\\\
2,\;\;\;\;\;\;\;\;\;\;\;\text{if}\;\;|m|>2.\end{cases}$ (36)

Figure 12: The value of the DOS at the edges of a band gap, i.e. at
$\Omega=\Omega_{L}$, vs. the parameter $m$. The Inset shows a dependence of
the gap width vs. $m$. The topological and trivial insulators are indicated in
the figures.

The dependence of the gap width $\Delta$ as a function of $m$ is shown in the
Inset of Fig. 12. It is seen that when $\Delta<2$ the gap of the same width is
opened for the six different values of the parameter $m$, at half-filling four
of these $m$ will correspond to topological insulators and other two $m$ will
correspond to the trivial insulators.

To discuss further we chose the case with the width $\Delta=1$, which is
possible for $|m|=0.5$, $1.5$, and $2.5$. The DOS for $|m|=1.5$ and $2.5$ are
shown in Fig. 13 whereas the DOS for $|m|=0.5$ is in Fig. 7. Interestingly, it
can be seen that the DOS at the edges of the band gap, i.e.,
$\rho(\Omega=\pm\Omega_{L})$, for the topological insulator ($|m|=0.5$ and
$1.5$) is larger than the one for the trivial case $|m|=2.5$. The DOS at the
edges of the band gap can be obtained analytically using the properties of the
elliptical integrals $K(0)=\pi/2$, namely

$\rho(\Omega_{L})=\frac{1}{2\pi}\begin{cases}{2|m|}/{\sqrt{1-m^{2}}},\;\;\;\;\;\;\;\;\;\text{if}\;\;|m|<1,\\\
{||m|-2|}/{({|m|-1})},\;\;\;\text{if}\;\;|m|\geq 1.\end{cases}$ (37)

The dependence of $\rho(\Omega_{L})$ as a function of $m$ is shown in Fig. 12.
For $|m|>2$ it is a finite smooth function. When $|m|<2$ the function is
singular and diverges at $m=\pm 1$.

Figure 13: Comparison of the two DOS of the QWZ model, which are characterized
by the same gap. Total density of states at $m=\pm 1.5$ (black dashed line) at
half-filling corresponds to topological insulator. Total density of states at
$m=\pm 2.5$ (green solid line) at half-filling corresponds to trivial
insulator. Enhancement of the DOS for topological case is seen.

##  VI Spectral moments of total DOS

In this Section we present results on the spectral moments of the total DOS.
The general formula for the spectral moment of the order $n$ reads

$M_{n}=\int{\Omega^{n}\rho(\Omega)d\Omega}.$ (38)

Integrals of the form (38) can be calculated analytically on the base of our
analytic results from section III. Some of the calculations are shown in
Appendix B to demonstrate details of the integration technique.

The moment of the zero order represents normalization integral $M_{0}=2$. All
moments of the odd orders are zero because the DOS is an even function of
$\Omega$.

The moments of the even orders $n$ are found as the polynomials in the
parameter $m$ of the same orders:

$\displaystyle M_{2}=2(m^{2}+2),$ $\displaystyle M_{4}=2(m^{4}+8m^{2}+5),$
$\displaystyle M_{6}=2(m^{6}+18m^{4}+51m^{2}+14),$ $\displaystyle
M_{8}=2(m^{8}+32m^{6}+210m^{4}+284m^{2}+42).$ (39)

This sequence can be continued if needed. Note, that all moments have a factor
of $2$, which corresponds to two symmetrical bands.

##  VII Conclusions and outlooks

In this paper we derived the analytic formulae of the DOS of the QWZ
Hamiltonian, a generic model for topological insulators in two dimensions. The
results are expressed in terms of the complete elliptic integrals. Analytic
expressions for the DOS are rare in general. Our results extend the class of
models where the analytic DOS is known.

We discussed in details the plots of the DOS and compared them with the
dispersion relations for the same value of the parameter $m$. Some additional
finite peaks in the DOS were identified. We provided explicit formulae for the
gap width and for the with of the bands in the QWZ model. We also found that
for the same gap width the topological system has larger DOS at the gap edge
as compared with the trivial case. Apparently, in the topological case a more
spectral weight is redistributed close to the band gap. Finally, we obtained
expressions of the spectral moments of the QWZ model. They are polynomials of
the parameter $m$, controlling the topology, and are of the same orders as the
orders of the corresponding moments.

The analytic DOS should be useful in determining thermodynamics or response of
the QWZ Chern topological insulator. It should also simplify the dynamical
mean-field theory study of the QWZ model when a Hubbard type of the
interaction is added to the Hamiltonian.

We hope that by using similar methods other DOS can be obtained in analytic
form, for example for the Haldane model on a hexagonal lattice.

##  VIII Acknowledgment

We thank for the financial support of the Excellence Initiative - Research
University (IDUB) via the grant under the program New Ideas - Ukraine. K.B.
also acknowledges the support of the Deutsche Forschungsgemeinschaft under the
Transregional Collaborative Research Center TRR360.

##  Appendix A Full elliptic integrals of the first and second kind

Expressions for the elliptic integrals, which we are using, are the following

$\displaystyle
K(q)=\int_{0}^{\pi/2}\frac{d\beta}{\sqrt{1-q^{2}\sin^{2}\beta}},$ (40)
$\displaystyle E(q)=\int_{0}^{\pi/2}{d\beta}{\sqrt{1-q^{2}\sin^{2}\beta}},$
(41)

where the argument is denoted by $q$ and called by mathematicians the modulus
of the elliptic integral and by $q^{\prime}$ we express the complementary
modulus, i.e., $q^{\prime 2}=1-q^{2}$.

##  Appendix B Calculation of spectral moments in details

Here we show how to calculate some of the integrals (38). Though the spectral
moment of any order can be easily calculated numerically, the analytical
derivation is of special value and methodological interest.

We start with the simple case of $|m|=1$. We express
$q^{2}={9-\Omega^{2}}/{8}$ and $q^{\prime 2}=(\Omega^{2}-1)/8$. For all values
of $\Omega^{2}$ from the region where $\rho(\Omega)$ is nonzero the modulus is
$0\leq q\leq 1$. In these terms the integral, normalizing the DOS, takes the
form

$\frac{M_{0}}{2}=\frac{4}{\pi^{2}}\int_{0}^{1}{K(q)dq^{\prime}}.$ (42)

Using the identity (6.141.2) in Gradshteyn ,

$\int_{0}^{1}{K(q^{\prime})dq}=\frac{\pi^{2}}{4}.$ (43)

gives the correct value of normalization.

The second moment takes the form

$\frac{M_{2}}{2}=\frac{4}{\pi^{2}}\int_{0}^{1}(9-8q^{2}){K(q)dq^{\prime}}.$
(44)

It can be rewritten as $M_{2}/2=M_{0}/2+32J/\pi^{2}$, where

$J=-\int_{0}^{1}q^{\prime}K(q)qdq.$ (45)

Here we used $q^{\prime}dq^{\prime}=-qdq$. The integral $J$ can be determined
by parts using the following substitution $dv=K(q)qdq$ and $u=q^{\prime}$. It
allows us to use the indefinite integral (5.112.3) in Gradshteyn , namely

$\displaystyle v=\int K(q)qdq=E(q)-q^{\prime 2}K(q),$ (46)

where $E(q)$ is the elliptic integral of the second kind. Integration by parts
gives $J=\int_{0}^{1}E(q)dq^{\prime}-J$, where it is used that
$E(0)=K(0)=\pi/2$ and $E(1)=1$. Using the identity (6.148.2) in Gradshteyn

$\int_{0}^{1}{E(q^{\prime})dq}=\frac{\pi^{2}}{8}$ (47)

we finally obtain $J={\pi^{2}}/{16}$ and $M_{2}=6$.

The moment of the order $n$ is expressed as a combination of integrals of the
product of an elliptic function and an even power of the modulus,

${M_{n}}=\sum_{i=0}^{n/2}C_{i}\int_{0}^{1}q^{2i}{K(q)}dq^{\prime}.$ (48)

Here $C_{i}$ are real numbers, compare with Eq. (44). So the same procedure
with integration by parts can be used to obtain moment of arbitrary order.

For $|m|>1$ the DOS have the infinite peaks at $\pm\Omega_{\infty}$. Integral
over $\Omega$ in Eq. (38) is split into two: from $\Omega_{L}$ to
$\Omega_{\infty}$ and from $\Omega_{\infty}$ to $\Omega_{R}$. In these regions
the complementary modulus $q^{\prime}$ take values from $1$ to $0$, and from
$0$ to $1$ respectively. Denoting $x=\Omega^{2}$ we can present the moment of
the order $n$ in the form

$\displaystyle{M_{n}}=\frac{\sqrt{2}}{\pi^{2}}\int_{0}^{1}dq^{\prime}{K(q)}\left[\frac{x_{2}^{n/2}}{\sqrt{x_{2}+m^{2}-2}}\left(\frac{dx_{2}}{dq^{\prime}}\right)-\right.$
$\displaystyle\left.\frac{x_{1}^{n/2}}{\sqrt{x_{1}+m^{2}-2}}\left(\frac{dx_{1}}{dq^{\prime}}\right)\right].$
(49)

The substitution $x=x_{1,2}(q)$ is not a single-valued function:
$x_{1,2}=m^{2}+4q^{\prime 2}\mp 4q^{\prime}\sqrt{m^{2}-q^{2}}$, where the sign
$"+"$ corresponds to the first region and $"-"$ to the second one. The
expression in the square brackets is transformed algebraically to a polynomial
function $\sum_{i=0}^{n/2}C_{i}(m^{2})q^{2i}$, so for $n=0$ it is equal to
$4\sqrt{2}$ and for $n=2$ it equal to $4\sqrt{2}(m^{2}+8q^{\prime 2})$. Then
moments are of the form (48), where coefficients $C_{i}=C_{i}(m^{2})$ are the
functions of $m^{2}$. The result of calculations leads to the expressions
(VI).

For the case $|m|<1$ two regions of integration need to be considered
separately: $|m|\leq\Omega\leq 2-|m|$ and $2-|m|\leq\Omega\leq 2+|m|$. In each
of these regions the substitution ${x}={x}({q})$ is not a single-valued, so
each integral splits into two ones, similar to the case with $|m|>1$.

In the first region we obtain the following: For
$|m|\leq\Omega\leq\sqrt{2-m^{2}}$ the substitution is
$x_{1}=m^{2}-4q(q-\sqrt{1-m^{2}q^{\prime 2}})/q^{\prime 2}$. The limits of
integration $x_{1}(q^{\prime}=1)=m^{2}$ and $\lim_{q^{\prime}\rightarrow
0}x_{1}=2-m^{2}$. For $\sqrt{2-m^{2}}\leq\Omega\leq 2-|m|$ the substitution is
$x_{2}=m^{2}+4(1-\sqrt{q^{2}+m^{2}q^{\prime 2}})/q^{\prime 2}$ with the limits
of integration $x_{2}(q^{\prime}=1)=(2-|m|)^{2}$ and
$\lim_{q^{\prime}\rightarrow 0}x_{2}=2-m^{2}$.

In the second region we have $x_{3,4}=m^{2}+4{q^{\prime}}^{2}\mp
4q^{\prime}\sqrt{m^{2}-q^{2}}$ valid for $2-|m|\leq\Omega\leq\Omega_{M}$ and
$\Omega_{M}\leq\Omega\leq 2+|m|$, respectively. Here $\Omega_{M}$ is defined
by the condition $dq/d\Omega=0$ and
$q^{2}={(\Omega^{2}-(m-2)^{2})((m+2)^{2}-\Omega^{2})}/{(\Omega^{2}+m^{2}-2)}/8$.
It can be proven that $q(\Omega_{M})=|m|$.

Thus the formula for the $M_{n}$ takes the form

$\displaystyle{M_{n}}=\frac{8}{\pi^{2}}\int_{0}^{1}dq^{\prime}{K(q)}\left[\frac{x_{2}^{n/2}}{{x_{2}+m^{2}}}\left(\frac{dx_{2}}{dq^{\prime}}\right)-\right.$
$\displaystyle\left.\frac{x_{1}^{n/2}}{\sqrt{(x_{1}-(m-2)^{2})(x_{1}-(m+2)^{2})}}\left(\frac{dx_{1}}{dq^{\prime}}\right)\right]+$
$\displaystyle\frac{\sqrt{2}}{\pi^{2}}\int_{0}^{|m|}dq{K(q)}\left[\frac{x_{3}^{n/2}}{\sqrt{x_{3}+m^{2}-2}}\left(\frac{dx_{3}}{dq}\right)-\right.$
$\displaystyle\left.\frac{x_{4}^{n/2}}{\sqrt{x_{4}+m^{2}-2}}\left(\frac{dx_{4}}{dq}\right)\right].$
(50)

It should be noted, that these integrals are difficult to calculate
analytically. For example, for $M_{0}$ the substitution $x=x(q)$ in Eq. (B)
results in

$\displaystyle({M_{0}}-1)\frac{\pi^{2}}{8}=\int_{0}^{|m|}dq{K(q)}\frac{q}{\sqrt{m^{2}-q^{2}}}+$
$\displaystyle\int_{0}^{1}{dq}\frac{K(q)}{q^{\prime
2}}\left[\frac{q}{\sqrt{m^{2}+q^{2}(1-m^{2})}}-\frac{q}{\sqrt{(1-m)^{2}+m^{2}q^{2})}}\right],$
(51)

where we used the identity (6.144) in Gradshteyn

$\displaystyle\int_{0}^{1}K(q)\frac{1}{1+q}dq=\frac{\pi^{2}}{8}.$ (52)

Knowing the fact that $M_{0}=2$, the expression (B) turns into an interesting
relation between integrals containing $K(q)$.

In order to determine these integrals, the integrands can be expanded into an
infinite series in $q^{2}$ followed by term-by-term integration. But the
easiest way to prove Eq. (VI) for the case $|m|<1$ is the numerical
integration.

## References

  * (1) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * (2) X-L. Qi and S-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * (3) R. B. Laughlin, Phys. Rev. B 23, 5632 (1981).
  * (4) J. Wang and S-C. Zhang, Nat. Materials 16, 1062 (2017).
  * (5) B.Q. Lv, T. Qian, and H. Ding, Rev. Mod. Phys. 93, 025002 (2021).
  * (6) N.R. Cooper, J. Dalibard, and I.B. Spielman, Rev. Mod. Phys. 91, 015005 (2019).
  * (7) L. Lu, J. D. Joannopoulos and M. Soljacic, Nat. Photonics 8, 821 (2014).
  * (8) A.P. Ramirez and B. Skinner, Physics Today 73 (9), 30 (2020).
  * (9) W. P. Su, J. R. Schrieffer, and A. J. Heeger, Phys. Rev. Lett. 42, 1698 (1979).
  * (10) M. J. Rice and E. J. Mele, Phys. Rev. Lett. 49, 1455 (1982).
  * (11) F.D.M. Haldane, Phys. Rev. Lett. 61, 2015 (1988).
  * (12) X-L. Qi, Y-S. Wu, and A-C. Zang, Phys. Rev. B 74, 085308 (2006).
  * (13) J. K. Asboth, L. Oroszlany, and A. Palyi, A Short Course on Topological Insulators, Lecture Notes in Physics (LNP, vol. 919) (Springer, Cham, 2016).
  * (14) D. Krüger and M. Potthoff, Phys. Rev. Lett. 126, 196401 (2021).
  * (15) J.P. Hobson and W. A. Nierenberg, Phys. Rev. 89, 662 (1953).
  * (16) V.O. Ananyev and M. I. Ovchynnikov, Cond. Mat. Phys. 20, 4, 43705 (2017).
  * (17) E. Kogan and G. Gumbs, Graphene 10, 1 (2021).
  * (18) Ch-K. Chiu, J.C.Y. Teo, A.P. Schnyder, and Sh. Ryu, Rev. Mod. Phys. 88, 035005 (2016).
  * (19) As will be seen latter, when the gap is closed at $m=0$ and $|m|=2$ the bands touch each other at separated point(s) and the DOS is zero at the corresponding energy $\Omega=0$. Hence, the band resolved DOS is determined uniquely from the total DOS.
  * (20) I.S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series and Products, (Academic Press is an imprint of Elsevier, USA, 2007).