\documentclass[aps,prl,twocolumn]{revtex4}
%\documentclass[aps,prb,showpacs]{revtex4}
\usepackage{bm,color,amsmath,amssymb,mathrsfs,latexsym,graphicx,psfrag}
%\addtolength{\textwidth}{3pt}
%\addtolength{\textheight}{0pt}
\renewcommand{\theequation}{S\arabic{equation}}
% boldsymbol (requires amsmath)
\newcommand{\bs}[1]{\boldsymbol{#1}}

% A command for inner product and bras and kets
\newcommand{\braket}[2]{\left\langle #1 | #2 \right\rangle}
\newcommand{\bra}[1]{\left\langle#1\right|}
\newcommand{\ket}[1]{\left|#1\right\rangle}
\newcommand{\bigket}[1]{\bigl|#1\bigr\rangle}
\newcommand{\textket}[1]{|#1\rangle}

% Various bracketing commands
\newcommand{\of}[1]{\!\left(#1\right)}
\newcommand{\sqof}[1]{\left[#1\right]}
\newcommand{\cuof}[1]{\left\{#1\right\}}

% commutator and anticommutator
\newcommand{\comm}[2]{\left[#1,#2\right]}
\newcommand{\anticomm}[2]{\left\{#1,#2\right\}}
% sum on nearest neighbor bonds
\newcommand{\bond}{\left\langle i, j \right\rangle}
%\newcommand{\bondsum}{\sum_{\left\langle i, j \right\rangle}}
\newcommand{\nbond}{\left\langle\left\langle i, j \right\rangle\right\rangle}
% 1/2
\newcommand{\half}{$\frac{1}{2}$ }
% simplifies using the up and down arrows to denote spin
\newcommand{\up}{\uparrow}
\newcommand{\dw}{\downarrow}
% Theta function
\newcommand{\tfunc}{\vartheta_1}
% notation for vacuum, an empty set inside a ket
\newcommand{\vac}{\left|\,0\,\right\rangle}
% Absolute value
\newcommand{\abs}[1]{\left|#1\right|}

% Roman functions for real and imaginary parts
\newcommand{\re}{\mathrm{Re}}
\newcommand{\im}{\mathrm{Im}}

% Sets of up-spin and down-spin locations
\newcommand{\bket}{\left\{z_1 \cdots z_{\num}\right\}}
\newcommand{\wket}{\left\{w_1 \cdots w_{\num}\right\}}

%Expectation values
\newcommand{\expect}[1]{\left\langle#1\right\rangle}

\def\ie{{\it i.e.},\ }
\def\eg{{\it e.g.}\ }
\def\ea{{\it et al.}}

\newcommand{\Hamil}{{\cal H}}
\newcommand{\bk}{{\mathbf k}}
\newcommand{\bq}{{\mathbf q}}
\newcommand{\br}{{\mathbf{r}}}
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}

\def\be{\begin{equation}}
\def\ee{\end{equation}}
\def\bea{\begin{eqnarray}}
\def\eea{\end{eqnarray}}
\def\Cal{\cal}

\begin{document}

\section{Derivation of $h_1(\bq)$ using the little group at $D_1$}
\label{sec:littlegroup}
The full symmetry group of the thin film in the absence of applied fields is $D_{4h}\otimes{T}$. The little group at a Dirac point $D_i$ is the subgroup of all operations that leave $D_i$ invariant. The little group therefore consists of a mirror plane that passes $\bar\Gamma{}D_i$, $C_{2T}$ and their combinations. Taking $D_1$ as example, the little group is generated by $M_{1\bar10}$ and $C_{2T}$. In a general spin-$1/2$ system we have: $M^2_{1\bar10}=C_2^2=T^2=-1$ and $\{M_{1\bar10},C_2\}=[M^2_{1\bar10},T]=[C_2,T]=0$, where $C_2$ is the 180-rotation about $[001]$-direction. Therefore the two generators satisfy (i) $M_{1\bar10}^2=-C_{2T}^2=-1$ (ii) $\{M_{1\bar10},C_{2T}\}=0$. There is only one 2D irreducible representation up to a basis rotation: $M_{1\bar10}=i\sigma_y$ and $C_{2T}=K\sigma_x$. Physically, $M_{1\bar10}$ relates the Hamiltonian $h_1(q_1,q_2)$ to $h_1(q_1,-q_2)$ and $C_{2T}$ commutes with $h_1(\bq)$; or mathematically, $M_{1\bar10}h_1(q_1,q_2)M^{-1}_{1\bar10}=h_1(q_1,-q_2)$ and $[C_{2T},h_1(\bq)]=0$. The irreducible representation of the little group along with the symmetry constraints determine the form of $h_1(\bq)$ shown in Eq.(1).

In general, the $k\cdot{p}$ model is given by
\bea
h_1(\bq)=d_0(\bq)I_{2\times2}+d_x(\bq)\sigma_x+d_y(\bq)\sigma_y+d_z(\bq)\sigma_z,
\eea
which must satisfy the symmetry constraints:
\bea
M_{1\bar10}h_1(q_1,q_2)M^{-1}_{1\bar10}&=&h_1(q_1,-q_2),\\
\nonumber
[C_{2T},h_1(\bq)]=0.
\eea
These symmetry constraints give that (1) $d_{0,y}$ is even under $q_2\rightarrow-q_2$, (2) $d_x$ is odd under $q_2\rightarrow-q_2$ and (3) $d_z=0$ to arbitrary order. We expand them to the second order in $|q|$:
\bea
d_0(q_1,q_2)&=&v_0q_1+\frac{q_1^2}{2m_1}+\frac{q_2^2}{2m_2},\\
\nonumber
d_x(q_1,q_2)&=&v_2q_2+\frac{q_1q_2}{2m_3},\\
\nonumber
d_y(q_1,q_2)&=&v_1q_1+\frac{q_1^2}{2m_4}+\frac{q_2^2}{2m_5}.
\eea
These terms make the dispersion deviate from perfectly linear and may be understood as the `warping' terms; they also make corrections to the wave functions at each $\bq$. It should be noted that $d_z=0$ holds up to arbitrary orders and this means there is no out-of-plain pseudo-spin component at any $\bq$. While including higher order terms explains the shape-changing of the equal energy contours from perfect ellipsoids, the Lifshitz transition cannot be described in the framework of any two-band theory. To do so, the model must be extended a four-band one, in order to account for the hybridization between nearest cones, as discussed in Ref.[21].

\section{Relating the four Dirac cones by $C_4$ symmetry}
\label{sec:fourDirac}
In the main text, we mention that by 90-degree rotations the effective theories for the four cones can be related. This is an intuitive statement yet to be made precise. In fact, $k\cdot{p}$ theories are always written with respect to a chosen basis, which is our case is furnished by (the periodic part of) the two Bloch states that are degenerate at the Dirac point. Due to the degeneracy, there is a gauge degree of freedom in the choice. Here the choice is made by fixing the little group representation at $D_1$: $M_{1\bar10}=i\sigma_y$ and $C_{2T}=K\sigma_x$. If we denote the two basis states by $|u_{1\uparrow}\rangle$ and $|u_{1\downarrow}\rangle$, we then fix the bases at $D_{2,1',2'}$ to be $\{|u_{2\uparrow}\rangle,|u_{2\downarrow}\rangle\}=\{\tilde C^2_4|u_{1\uparrow}\rangle,\tilde C^2_4|u_{1\downarrow}\rangle\}$, $\{|u_{1'\uparrow}\rangle,|u_{1'\downarrow}\rangle\}=\{\tilde C_4|u_{1\uparrow}\rangle,\tilde C_4|u_{1\downarrow}\rangle\}$ and $\{|u_{2'\uparrow}\rangle,|u_{2'\downarrow}\rangle\}=\{\tilde C^3_4|u_{1\uparrow}\rangle,\tilde C^3_4|u_{1\downarrow}\rangle\}$, respectively. Mark that here $\tilde C_4$ is the matrix representing the 90-degree rotation in both orbital space (including spin). Defining the Bloch wave function at $\mathbf{D}_i+\bq$ as $|\psi_{i\uparrow/\downarrow}(\bq)\rangle=e^{i(\mathbf{D_i}+\bq)\cdot\br}|u_{i\uparrow/\downarrow}\rangle$, it is easy to check that $|\psi_{2}(\bq)\rangle=\hat{C}^2_4|\psi_1(-\bq)\rangle$, $\hat{C}_4|\psi_{1'}(\bq)\rangle=\hat{C}_4|\psi_1(-q_2,q_1)\rangle$ and $|\psi_{2'}(\bq)\rangle=\hat{C}_4^3|\psi_1(q_2,-q_1)\rangle$. Here $\hat{C}_4$ is the single particle operator acting in the Hilbert space, which is the combination of the orbital rotation $\tilde{C}_4$ plus rotation $(x,y)\rightarrow(-y,x)$, where $(x,y)$ is a lattice point and the rotation center is also placed at a lattice point. The full single Hamiltonian, projected to the states at the vicinities of the four Dirac points, is given by
\bea
\hat{H}=\sum_{\bq,i=1,2,1',2',\alpha,\beta=\uparrow,\downarrow}(h_i(\bq))^{\alpha\beta}|\psi_{i\alpha}(\bq)\rangle\langle\psi_{i\beta}(\bq)|.
\eea
$C_4$ symmetry implies $[\hat{C}_4,\hat{H}]=0$, which immediately leads to $h_2(q_1,q_2)=h_1(-q_1,-q_2)$, $h_{1'}(q_1,q_2)=h_1(-q_2,q_1)$ and $h_{2'}(q_1,q_2)=h_1(q_2,-q_1)$, confirming the intuitive relations appearing in the main text.

\section{Calculation of the Chern number of the top/bottom surface}
\label{sec:Chern}
In the text we refer to the Chern number contributed by one massive Dirac cone, which is not mathematically well-defined. In fact, the integrated Berry's curvature of a gapped Dirac cone is non-quantized in any finite $\bk$-space, hence possesses no well-defined Chern number. The Chern number of a whole 2D surface (top surface for example) is, however, a well-defined quantity (if periodic boundary is taken for the other two directions), which may be calculated. Suppose we are interested in the Chern number, $C$, at some Zeeman field $\Delta_Z=\Delta_0>0$. Then since time-reversal reverses the Chern number, we know for $\Delta_Z=-\Delta_0$, the Chern number must be $-C$. Consider a 3D space spanned by $q_{1,2}$ and $\Delta_Z$, then from Gauss's law, the Chern number change from $\Delta_Z=-\Delta_0$ to $\Delta_0$ equals the total monopole charge between these two planes in the 3D parameter space. The monopole, or gap closing point, is always at $(q_1,q_2,\Delta_Z)=0$, around which the Hamiltonian is that of 3D Weyl fermions: $h(q_1,q_2,q_3)=\sum_{i,j=1,2,3}A_{ij}\sigma_iq_j$, where $q_3\equiv{\Delta_Z}$. The charge of such a monopole is $\textrm{sign}\det(A)$, and since there are in total four such monopoles between $\Delta_Z=\pm\Delta_0$, we have the difference in Chern number $C-(-C)=4\textrm{sign}(\det{}A)$, or $C=2\textrm{sign}(\det{}A)$. All Chern numbers obtained in the text are derived using this method.

\section{Diagonalizing the Hamiltonian in Eq.(4)}
\label{sec:dispersion}
A Hamiltonian that describes isolated top and surface states around $D_i$ is
\bea
\tilde{H}_i=\left(\begin{matrix}
H^t_i & 0\\
0 & H^b_i
\end{matrix}\right),
\eea
and hybridization is equivalent to adding an off-diagonal block term, resulting in, to the lowest order in $|q|$,
\bea
\tilde{H}_i=\left(\begin{matrix}
H^t_i & \Delta_HI_{2\times2}\\
\Delta_HI_{2\times2} & H^b_i
\end{matrix}\right).
\eea
Diagonalizing $\tilde{H}_1$ directly, we obtain four bands:\begin{widetext}
\bea
E_1(\bq)&=&v_0q_1+\sqrt{{\Delta^t_1}^2+{\Delta^t_1}^2+2\Delta_H^2+2q_1^2v_1^2+2q_2^2v_2^2+\sqrt{(\Delta^t_1-\Delta^b_1)^2+4[(\Delta^t_1+\Delta^b_1)^2+4v_1^2q_1^2+4v_2^2q_2^2]}},\\
\nonumber
E_2(\bq)&=&v_0q_1+\sqrt{{\Delta^t_1}^2+{\Delta^t_1}^2+2\Delta_H^2+2q_1^2v_1^2+2q_2^2v_2^2-\sqrt{(\Delta^t_1-\Delta^b_1)^2+4[(\Delta^t_1+\Delta^b_1)^2+4v_1^2q_1^2+4v_2^2q_2^2]}},\\
\nonumber
E_3(\bq)&=&v_0q_1-\sqrt{{\Delta^t_1}^2+{\Delta^t_1}^2+2\Delta_H^2+2q_1^2v_1^2+2q_2^2v_2^2-\sqrt{(\Delta^t_1-\Delta^b_1)^2+4[(\Delta^t_1+\Delta^b_1)^2+4v_1^2q_1^2+4v_2^2q_2^2]}},\\
\nonumber
E_4(\bq)&=&v_0q_1-\sqrt{{\Delta^t_1}^2+{\Delta^t_1}^2+2\Delta_H^2+2q_1^2v_1^2+2q_2^2v_2^2+\sqrt{(\Delta^t_1-\Delta^b_1)^2+4[(\Delta^t_1+\Delta^b_1)^2+4v_1^2q_1^2+4v_2^2q_2^2]}}.
\eea\end{widetext}
Straightforward algebraic work shows that the \emph{only} solution for $E_2(\bq)=E_3(\bq)$, i.e., a gap-closing point, exists at $q_1=q_2=0$ when $|\Delta_H|=\sqrt{\Delta_1^t\Delta_1^b}$.

Parallel discussion for $D_{2,1',2'}$ proceeds and we conclude that a topological phase transition happens when 
\bea\label{eq:transition}
|\Delta_H|=\sqrt{\Delta^t_i\Delta_i^b},
\eea
whereas the Chern number contributed by the cone at $D_i$ changes from $\pm1$, depending on the sign of $\Delta_i^{t,b}$, to zero. Mark that on the right hand side of Eq.(\ref{eq:transition}), if $\Delta^t_i\Delta^b_i<0$, the transition cannot happen at any $\Delta_H$.

\section{Derivation of Table I}

Here we provide details for deriving Table I. 

First we note the following two simple facts: (i) the strain tensor is a rank-2 tensor so its general transformation under $O(3)$ takes the form:
\bea\label{eq:tensor}
\epsilon_{ab}\rightarrow{R}_{aa'}R_{bb'}\epsilon_{a'b'},
\eea
where $R$ is the three-by-three rotation matrix,
and (ii) it is invariant under time-reversal. All point group operators including $C_{2,4}$ and $M_{110,1\bar10}$ are elements of $O(3)$ and due to (ii), $C_{2T}$ transforms the tensor in the same way as $C_2$ does. The rotation matrix for these operations are given by
\bea\label{eq:Rmatrices}
R_{C_2}&=&\left(\begin{matrix} % or pmatrix or bmatrix or Bmatrix or ...
      -1 & 0 & 0\\
      0 & -1 & 0\\
      0 & 0 & 1\\
   \end{matrix}\right),
R_{M_{110}}=\left(\begin{matrix} % or pmatrix or bmatrix or Bmatrix or ...
      -1 & 0 & 0\\
      0 & 1 & 0\\
      0 & 0 & 1\\
   \end{matrix}\right),\\\nonumber
R_{M_{1\bar10}}&=&\left(\begin{matrix} % or pmatrix or bmatrix or Bmatrix or ...
      1 & 0 & 0\\
      0 & -1 & 0\\
      0 & 0 & 1\\
   \end{matrix}\right),
R_{C_4}=\left(\begin{matrix} % or pmatrix or bmatrix or Bmatrix or ...
      0 & -1 & 0\\
      1 & 0 & 0\\
      0 & 0 & 1\\
   \end{matrix}\right).
\eea
Substituting Eq.(\ref{eq:Rmatrices}) back into Eq.(\ref{eq:tensor}), we have the upper part of Table I.

We then can use the upper part of Table I to determine the terms that represent spin in the effective theories for the four Dirac cones. The strain terms for the states around $D_{i=1,2,1',2'}$ are generally written as
\bea\nonumber
\delta{H}_{ab,i}=m^0_{ab,i}\sigma_0+m^x_{ab,i}\sigma_x+m^y_{ab,i}\sigma_y+m^z_{ab,i}\sigma_z+O(|q|).\\
\eea
The little group of $D_i$ gives constraints on $m^\alpha_{ab,i}$. Hereafter we ignore $m^0_{ab,i}$ because this does not change the wavefunction, thus preserving the topology, of the surface states. $\epsilon_{11,22,33}$ are unchanged under $C_{2T}$, $M_{110}$ or $M_{1\bar10}$, from which we know $m^x_{11/22/33,i}=m^z_{11/22/33,i}=0$. Since $C_4$ relates the four Dirac points and sends $\epsilon_{11/22}$ to $\epsilon_{22/11}$, we know $m^y_{11/22,1}=m^y_{22/11,1'}=m^y_{11/22,2}=m^y_{22/11,2'}$ and $m^y_{22/11,1}=m^y_{11/22,1'}=m^y_{22/11,2}=m^y_{11/22,2'}$. So far we have derived the first three columns of the lower part of the Table.

$\epsilon_{12}$ changes sign under $M_{110,1\bar10}$ and invariant under $C_{2T}$, from which we have $m^y_{12,i}=m^z_{12,i}=0$. Again from $C_4$, we obtain $m^x_{12,1}=-m^x_{12,1'}=m^x_{12,2}=-m^x_{12,2'}$. This is the fourth column.

$\epsilon_{13/23}$ are invariant under $M_{1\bar10/110}$ but changes sign under $M_{110/1\bar10}$ and $C_{2T}$. The little group at $D_{1,2}$ is $\{M_{1\bar10},C_{2T}\}$, so we have $m^{x,z}_{13,1}=m^{x,z}_{13,2}=0$ and $m^{x,y}_{23,1}=m^{x,y}_{23,2}=0$. Using that $C_2$ maps $D_{1/2}$ to $D_{2/1}$, we have $m^y_{13,1}=-m^y_{13,2}$ and $m^z_{13,1}=-m^z_{13,2}$. Then we notice that since $C_4$ maps $D_{1/1'}$ to $D_{2/2'}$, and $\epsilon_{13}$ ($\epsilon_{23}$) to $\epsilon_{23}$ ($-\epsilon_{13}$), there are $m^{x,y,z}_{23,1'/2'}=m^{x,y,z}_{13,1/2}$. This gives the last two columns.

The above relations give in total six free masses: $(m^y_{11,1},m^y_{22,1},m^y_{33,1},m^x_{12,1},m^y_{13,1},m^z_{23})$, defined as $(\lambda_{11},\lambda_{22},\lambda_{33},\lambda_{12},\lambda_{13},\lambda_{23})$.

\end{document}