a TI thin film of infinite dimensions along the $x$ and $y$ directions, and small finite thickness along the $z$ direction, can be written as $$H = (\vec{k}\times\vec{\sigma})\cdot\hat{z}\tau_z + M_y\sigma_y + \lambda \tau_x. \label{tfHam0}$$

We use units where $e=\hbar=v_f=1$.

We first highlight the influence of the inter-surface coupling term $\lambda$ on the energy spectrum. Consider the limit where $\lambda \rightarrow 0$,$M_y \neq 0$. In this case, the upper and lower surfaces may be considered separately, and the energy spectrum consists of two Dirac cones. The states localized near the upper surface have $\langle \tau_z \rangle = +1$, while the state localized near the lower surface have the opposite sign of $\langle \tau_z \rangle$. The $M_y\sigma_y$ term, however, has the same sign for both the upper and lower surfaces. The Dirac points for the Dirac cones for the upper surface states and the lower surface states are hence displaced in opposite directions in $k$ space.

![ Panel (a) shows the dispersion relations for the two values of $\lambda$ indicated in the legend at $k_y = 0$ for $\vec{M} = 0.5\hat{y}$. The three horizontal dotted lines correspond to the values of energies at which the EECs in panels (b) to (d) at 