\section{lattice models}
\emph{Winding numbers and geometry of the loops for a lattice model.-}
In order to reveal the topological properties described by the NL winding number, we next consider a lattice model described by an anticommuting set of Dirac matrices $\mathbf{\Gamma}=(\sigma_x s_0,\sigma_z s_0,\sigma_y s_x,\sigma_y s_y,\sigma_y s_z)$, as
\begin{eqnarray}
H&=&\mathbf{h}(k)\cdot\mathbf{\Gamma},\label{total_H}
\end{eqnarray}
with
\begin{eqnarray}
\mathbf{h}(k)&=&(h_{x0},h_{z0},h_{yx},h_{yy},h_{yz}),\nonumber\\
h_{x0}&=&\mu-t_{\parallel}(\cos{k_x}+\cos{k_y}),\nonumber\\
h_{z0}&=&-t_{\perp}\cos{k_z},\nonumber
\end{eqnarray}
which form two NLs in $k_x-k_y$ plane for $k_z=\pi/2$ and $-\pi/2$ when $0<|\mu/t_{\parallel}|<2$.
The position and shape of the NLs are only associated with the ratio of
$\mu$ and $t_{\parallel}$, hence we can choose $t_{\parallel}=t_{\perp}=1$ for
the sake of simplicity.
In the following we only consider the case with positive $\mu$, with the center of the NLs
given by $k_x=k_y=0$. For negative $\mu$, the center of the NLs is at $k_x=k_y=\pi$,
and a similar discussion applies.
The two NLs of $k_z=\pm\pi/2$ give two independent winding numbers, $\nu_{NL}^{\pm}$,
respectively, and we define the total winding number of the system as
\begin{eqnarray}
\nu_{sum}=\nu_{NL}^+ +\nu_{NL}^-.
\end{eqnarray}
Without loss of generality, here we choose $h_{yx}=0$ to preserve a chiral symmetry
with the operator $S=\sigma_y s_x$. The other two gap terms of $\mathbf{h}_y$ are functions of $k$,
and $\mathbf{h}_y=0$ gives 1D lines in the BZ.
We consider the following form of $\mathbf{h}_y$ that breaks time reversal symmetry:
\begin{eqnarray}
h_{yy}&=&-t_y\sin{k_y}\\
h_{yz}&=&\mu_s-t_x\sin{k_x}-t_z\sin{k_z},
\end{eqnarray}
and the system falls into the AIII class.
$|\mathbf{h}_y|=0$ gives some 1D lines in the $k_x-k_z$ plane with
$k_y=0$ or $\pi$, which may or may not be enclosed by the NLs.
We would also like to point out that although the gap terms need be small
for the system to preserve a NL like structure,
{\it the topological properties are not related to the exact value of these terms,
but only to the ratios between them}.
For the sake of simplicity, we choose $t_y=t_x=1$ and positive $\mu_s$ hereafter.

We first consider a simple case with $t_z=0$. In this case, the 1D solution of
$\mathbf{h}_y=0$ gives four lines perpendicular to $k_x-k_y$ plane,
two with $k_y=0$ and two with $k_y=\pi$. For positive $\mu$, the pair of lines with
$k_y=\pi$ are always outside the loop.
By tuning $\mu_s$ and $\mu$, the NLs may enclose
$2$, $1$ or $0$ lines with $k_y=0$, as shown by Fig.\ref{fig2}(a)-(c).
However, the windings of these lines in the $k_x-k_y$ plane have opposite values
(as in Fig.\ref{fig1}),
and the NL
enclosing either $0$ or $2$ lines will result in $\nu_{NL}=0$.
On the other hand,
as the system preserves a reflection symmetry along $z$ direction,
each line will be enclosed by either two or none of the NLs,
hence the total winding number, $\nu_{sum}$, in this case is always even.

\begin{figure}
\includegraphics[width=0.8\linewidth]{fig2.eps}
\caption{The NLs and the 1D lines given by $|\mathbf{h}_y|=0$,
represented by blue and red lines respectively.
For the latter, we only show the lines with $k_y=0$,
as the ones for $k_y=\pi$ lie outside the NLs.
(a) $\mu=0.2$, $\mu_s=0.8$ and $t_z=0$; (b) $\mu=1$, $\mu_s=0.8$ and $t_z=0$; (c) $\mu=1.8$, $\mu_s=0.8$ and $t_z=0$; (d) $\mu=1.6$, $\mu_s=0.5$ and $t_z=0.2$; (e) $\mu=1.6$, $\mu_s=0.5$ and $t_z=0.4$; (f) $\mu=1.6$, $\mu_s=0.5$ and $t_z=0.6$. }\label{fig2}
\end{figure}

In the presence of a nonzero $t_z$,
the reflection symmetry is broken, and
the lines of $|\mathbf{h}_y|=0$ will change shape with $t_z$ and eventually form a closed ring,
as shown in Fig.\ref{fig2}(d)-(f).
In this case, an enclosed line will cross one of the NLs at some point,  resulting
in a topological phase transition.
After this transition, the system has an odd winding number, $\nu_{sum}=1$,
as only one of the NLs encloses a singularity.

\emph{edgestates and phase diagram.-}
The topological properties of a 3D topological insulator can be represented by the number of
Dirac cones in the edge states under OBC.
Next, we apply the method in Ref.\cite{edge_state} to study the edge states in our model.
The existence of edge states and their eigenenergies,
under OBC,  are associated with the bulk topology
of the system, which can be seen by the trajectory of $\mathbf{h}(k)$ in the 5-component
vector space formed by the Dirac matrices $\mathbf{\Gamma}$.
Here we choose a surface plane perpendicular to the $x$ direction by fixing $k_y$ and $k_z$,
and study the corresponding edge states as an example.
Edge states in the other two directions can also be studied in this way,
and give similar results are obtained.
The Hamiltonian terms associated with $k_x$ give an elliptical
 trajectory of $\mathbf{h}$
in the 1-5 plane in $\mathbf{\Gamma}$ space,
\begin{eqnarray}
\mathbf{h}(k_x)=(-\cos{k_x},0,0,0,-\sin{k_x}).
\end{eqnarray}
The remaining Hamiltonian terms,
\begin{eqnarray}
\mathbf{h}^0=(\mu-\cos{k_y},-\cos{k_z},0,-\sin{k_y},\mu_s-t_z\sin{k_z}),
\end{eqnarray}
can be viewed as the vector from the origin of the vector space to the center of
the ellipse $\mathbf{h}(k_x)$. The parallel and perpendicular components of $\mathbf{h}^0$ to the 1-5 plane are given by
\begin{eqnarray}
\mathbf{h}^0_{\parallel}&=&(\mu-\cos{k_y},0,0,0,\mu_s-t_z\sin{k_z})\,,\\
\mathbf{h}^0_{\perp}&=&(0,-\cos{k_z},0,-\sin{k_y},0).
\end{eqnarray}
The existence of edge states depends on whether the ellipse
$\mathbf{h}(k_x)$ encloses the point $\mathbf{h}^0_{\parallel}$,
and this condition reads
\begin{eqnarray}
|\mathbf{h}^0_{\parallel}|=\sqrt{(\mu-\cos{k_y})^2+(\mu_s-t_z\sin{k_z})^2}<1\,.\label{condition}
\end{eqnarray}
Provided that Eq.(\ref{condition}) holds, the edge state energies
are given by
\begin{eqnarray}
E_{\pm}=|\mathbf{h}^0_{\perp}|=\pm\sqrt{\cos^2{k_z}+\sin^2{k_y}}.
\end{eqnarray}

Candidate Dirac cones at  $(k_y,k_z)=(0,\pm\frac{\pi}{2})$
must satisfy the inequality (\ref{condition}).
When $t_z=0$, Eq.(\ref{condition}) becomes
\begin{eqnarray}
(\mu-1)^2<1-\mu_s^2\,,
\end{eqnarray}
for both $k_z=\pm\frac{\pi}{2}$. In other words, it gives either no Dirac cone, or a pair of
Dirac cones at  $(k_y,k_z)=(0,\pm\frac{\pi}{2})$.
For nonzero $t_z$, there may exist $0$, $1$ or $2$ Dirac cones depending on the parameters, as the condition (\ref{condition}) becomes
\begin{eqnarray}
(\mu_s\mp t_z)^2<2\mu-\mu^2\,,
\end{eqnarray}
for $k_z=\pm\frac{\pi}{2}$, respectively. In Fig.\ref{fig2} we display phase diagrams showing
the number of Dirac cone edge states. These results are also
in consistence with the winding numbers of the NLs, as shown in the figure.

\begin{figure}
\includegraphics[width=0.45\linewidth]{fig3_1.eps}
\includegraphics[width=0.45\linewidth]{fig3_2.eps}
\caption{Phase diagram with $t_z=0$ (left) and $\mu=1$ (right). The number of Dirac cones in edge states equals to the total winding number of the model, $\nu_{sum}=\nu_{NL}^+ +\nu_{NL}^-$.
In the left panel, the yellow region with $\nu_{NL}^+=\nu_{NL}^-=0$ includes three different situations: i) two lines are enclosed by the NLs when $\mu_s<1$ and $\mu<1$; ii) all the lines are out of the NLs when $\mu_s<1$ and $\mu>1$; and iii) not any line exsits when $\mu_s>1$.}\label{fig3}
\end{figure}

In order to visualize the edge states, next we choose open boundary condition in $x$ direction and rewrite the Hamiltonian as
a tight-binding between planes
\begin{eqnarray}
H&=&\sum_{n}\hat{c}^{\dagger}_{n}U_d\hat{c}_{n}+\sum_{n}\frac{1}{2}\hat{c}^{\dagger}_{n+1}U_{od}\hat{c}_{n}+h.c.,
\end{eqnarray}
with $\hat{c}_n$ is a vector of annihilation operators $\hat{c}_{n,\sigma,s}$ on
plane $n$,
$\hat{c}_n=(\hat{c}_{n,+,+},\hat{c}_{n,+,-},\hat{c}_{n,-,+},\hat{c}_{n,-,-})_n$, and
\begin{eqnarray}
U_d&=&(\mu-\cos{k_y})\sigma_x s_0-\cos{k_z}\sigma_z s_0\nonumber\\
&&-\sin{k_y} \sigma_y s_y+(\mu_s-t_z\sin{k_z})\sigma_y s_z\,,\\
U_{od}&=&-\sigma_x s_0+i\sigma_y s_z\,.
\end{eqnarray}
We numerically diagonalize this Hamiltonian and show the four closest
doubly degenete energy bands above and below zero energy in Fig.\ref{fig4}, with
OBC along $x$ from (a) to (c), and periodic boundary condition from (d) to (f) as comparison.
Panels (a) and (b) are for the topologically nontrivial phases with $\nu_{NL}^+=\nu_{NL}^-=1$
and $\nu_{NL}^+=1,~\nu_{NL}^-=0$,
and the spectra show edge states with two or one Dirac cone, respectively.
Panel (c) is for the topologically trivial phase with $\nu_{NL}^+=\nu_{NL}^-=0$,
where there is no edge state connecting the conduction and valence bands.

\begin{figure}
\includegraphics[width=1\linewidth]{fig4.eps}
\caption{The spectra of the two doublet bands nearest to $E=0$, (a)-(c) for open boundary condition along $x$, and (d)-(f) for period boundary condition. The three columns are for three topologically different phases, while the left two are topologically nontrivial. The parameters are $\mu=1$ and (a), (d) $\mu_s=t_z=0$; (b), (e) $\mu_s=t_z=1$; and (c), (f) $\mu_s=2$, $t_z=0$.}\label{fig4}
\end{figure}

\emph{A real material example.-}Finally, we apply our method to 3D topological insulators of the
Bi$_2$Te$_3$ family\cite{real}.
These materials possess time reversal symmetry and are characterized by a $Z_2$ number.
Nevertheless, their topological nature is determined by the physics near the
time-reversal-invariant point $\Gamma(0,0,0)$ in the Brillouin zone,
around which the low-energy effective Hamiltonian also satisfies a chiral symmetry.
This Hamiltonian is given by
\begin{eqnarray}
H(k)=\epsilon_0(k)I_{4\times4}+\left(
\begin{array}{cccc}
M(k) & A_1k_z & 0 & A_2k_-\\
A_1k_z &-M(k)  &A_2k_-  & 0\\
0 &A_2k_+  &M(k)  & -A_1k_z\\
A_2k_+  &0  &-A_1k_z  & -M(k)
\end{array}\right),
\end{eqnarray}
with $k_{\pm}=k_x\pm ik_y$, $\epsilon_0(k)=C+D_1k_z^2+D_2|k_{+}|^2$ and
$M(k)=M-B_1k_z^2-B_2|k_+|^2$.
Using Dirac matrices, this Hamiltonian can be written as
\begin{eqnarray}
H(k)=M(k)\sigma_z s_0+A_1k_z\sigma_x s_z+A_2k_x\sigma_x s_x+A_2k_y\sigma_x s_y,
\label{bi2te3}
\end{eqnarray}
where we left out the identity matrix as it only changes the shape of the energy bands,
not the topology of the system.
The chiral operator is given by the absent fifth Dirac matrix,
$S=\sigma_y s_0$.

$P*T$ symmetry is here implemented by $P=\sigma_z$. However, we note that
particle-hole symmetry also exists in this case, which reads
$CH^*(-k)C^{-1}=-H(k)$, with $C=\sigma_y s_y$,
and satisfies $CC^*=1$.
Thus the model Eq.(\ref{bi2te3}) falls into the DIII class,
which is also characterized by a $Z$ invariant in 3D\cite{classes}.
Similar to our previous discussion, we write  $H=H_1+H_2$, with
\begin{eqnarray}
H_{1}(k)&=&M(k)\sigma_z s_0+A_1k_z\sigma_x s_z\,,\\
H_{2}(k)&=&A_2k_x\sigma_x s_x+A_2k_y\sigma_x s_y\,.
\end{eqnarray}
Then $H_{1}$ has a NL in $k_x-k_y$ plane, and the nodes of $H_{2}$
produce a single line enclosed by the loop.
This gives a NL winding number $\nu_{NL}=1$,
which corresponds to the topological properties of the $\Gamma$ point.

\emph{Summary.-}
In summary,
we have studied Hamiltonians with $P*T$ symmetry where two terms define a NL which is gapped
out by the other terms.
In the presence of chiral symmetry, these gap terms can drive the system into a chiral TI,
which can be described by a winding number defined along the NL.
This winding number is associated with the geometric relation between the NL
and the zeroes of the gap terms.
We investigated a  lattice model in detail, which has two NLs in the BZ,
each of them with  a winding number of $1$ or $0$
due to the gap terms. This winding number corresponds to the emergence of a Dirac
cone for the surface states under OBC.
Finally, we applied our method to the 3D topological insulators of the Bi$_2$Te$_3$ family,
and showed the connection between their topological nature and the NL winding number.


\emph{Acknowledgments.-}
Financial support from FCT through grant UID/CTM/04540/2013 is acknowledged.
S. C. is supported by NSFC under Grants No. 11425419, No. 11374354 and No. 11174360.

