Chiral Majorana edge states in HgTe quantum wells

Firstly, we show that BIA-SOC is sufficient in order to lift the (Kramer) spin-degeneracy [3] by freezing out half the degrees of freedom (unless k = 0), which otherwise would make the spectrum doubly degenerate. Under periodic boundary conditions, k_x and k_y are good quantum numbers and the energy dispersion for the standard BHZ Hamiltonian supplemented with the BIA term can be calculated to

$$\begin{equation} E(k_x,k_y)=\epsilon(k^2)\pm\sqrt{(A|k|\pm h)^2+ M(k^2)^2}. \end{equation} \tag{ 4 }$$

Clearly, the spectrum is non-degenerate for |k| ≠ 0. By additionally introducing the Zeeman terms, the time-reversal symmetry is broken and thus the Kramer degeneracy at |k| = 0 is lifted.

The combined effect of BIA and Zeeman terms on the bulk energy dispersion is visualized in figure 2 (left). In the following, we choose the Fermi energy E_F to be zero so that the Fermi level is situated within the gap between the E1-like bands opened by the Zeeman splitting. The doping parameter is tuned in such a way that C ≈ −M. As a result, only a single pair of Fermi points exists and thus the system has effectively spinless dynamics.

Zoom In Zoom Out Reset image size

When, in addition, the proximity of the s-wave superconductor is taken into account, electrons with energy around the Fermi level and momenta k and −k couple, opening a gap around the Fermi level as shown in figure 2 (right). The special conditions under which this gap closes will be discussed in the following. The Fermi points at which this happens are the key to the emergence of the chiral Majorana edge states [32].

Before calculating the Fermi points, we will briefly review the symmetry properties of the two E1-like bands closest to the Fermi level. Their bulk energy spectrum is presented in figure 3 in the presence of different SOC terms: BIA or Rashba-SOC or the cooperation of BIA and Rashba-SOC. Interestingly, although the spectrum is rotationally symmetric when either Rashba-SOC or BIA are present, this rotational symmetry is broken when Rashba-SOC and BIA are simultaneously present.

Zoom In Zoom Out Reset image size

This can be understood as follows [33]: while the effective magnetic fields due to Rashba-SOC are oriented clockwise and perpendicular to the wave vector k, the effective magnetic fields due to the BIA spin splitting have tetrahedral symmetry, as shown in figure 3. When both BIA and Rashba-SOC are present, the effective magnetic fields add vectorially, decreasing or enhancing each other. The energy of the lowest-energy band is very sensitive to the effective k-dependent in-plane magnetic field. Therefore, it exhibits tetrahedral symmetry when BIA- and Rashba-SOC terms interfere with each other.

In the following, we will identify the Fermi points of the spectrum. To simplify analytical calculations, we project the full Hamiltonian (1) onto an effective Hamiltonian. When M defines the dominant energy scale and the doping is chosen as C ≈ −M, four of the eight bands become irrelevant. Using quasi-degenerate perturbation theory [33] ('Löwdin partitioning'), the eight-band Hamiltonian (1) can be reduced to an effective four-band Hamiltonian acting on the E1-like electron and hole states only (see appendix)

$$\begin{eqnarray}&&\fl \skew3\tilde{H}=\mu(\hat{k}^2)\tau_z+B_E\sigma_z+\Delta_E\tau_x+R_0\hat{k}_x\sigma_y\tau_z-R_0\hat{k}_y\sigma_x\tau_z +\tilde{h}\hat{k}_x\left(\cos{2\theta}\sigma_x-\sin{2\theta}\sigma_y\right)\tau_z\nonumber \\ &&-\tilde{h}\hat{k}_y\left(\cos{2\theta}\sigma_y+\sin{2\theta}\sigma_x\right)\tau_z, \end{eqnarray} \tag{ 5 }$$

where the new parameters $\mu ({\hat k}^2)=\epsilon ({\hat k}^2)+M({\hat k}^2) +(A^2{\hat k}^2+h^2)/(2M)$ and $\tilde {h}=Ah/M$ have been introduced.

In the absence of BIA-SOC, the reduced form (5) is analogous to the effective Hamiltonians of different semiconductor systems, which have already been shown to host chiral Majorana edge states [9–13]. Note that the Rashba-SOC energy U_R ≡ 2m*R²₀ [18] is kept quite small here intentionally in order to study the combined effect of BIA and Rashba-SOC: for the values in table 1, it amounts to $U_{\mathrm { R}}=R_{0}^{2}/(\frac {A^2}{2M}-D-B) \approx 0.3\,\mathrm {meV}$. On the other hand, the BIA-SOC leads to an effective Dresselhaus-SOC in the reduced model (5), which amounts to $U_{\mathrm {D}}=(Ah/M)^{2}/(\frac {A^2}{2M}-D-B) \approx 0.7\,\mathrm {meV}$. In [34], electric fields up to $\mathcal {E}{\sim }100\,\mathrm {mV}\,\mathrm {nm}^{-1}$ were estimated to be possible, which could enhance the Rashba-SOC energy in principle by a factor of 30, without affecting the effective Dresselhaus energy.

The effective Hamiltonian (5) yields (for θ = π/4 and Δ_E ≠ 0) zero energy level crossings at

$$\begin{equation} (R_0-\tilde{h})^2k_x^2+(R_0+\tilde{h})^2k_y^2=-\left(|\Delta_E|\pm\sqrt{B_E^2-\mu(k^2)^2}\right)^2. \end{equation} \tag{ 6 }$$

This implies that for non-vanishing SOC, Fermi points exist (for B²_E − μ(k²)² > 0) at

$$\begin{equation} k_x=0, \quad k_y=0 \ \textrm{and} \ |\Delta_E|=\sqrt{B_E^2-\mu_0^2}, \end{equation} \tag{ 7 }$$

$$\begin{equation} k_x=0,\quad k_y\neq 0, R_0=-\tilde{h} \ \textrm{and} \ \textrm{(6),} \end{equation} \tag{ 8 }$$

$$\begin{equation} k_x\neq 0,\quad k_y=0, R_0=\tilde{h} \ \textrm{and} \ \textrm{(6),} \end{equation} \tag{ 9 }$$

where μ₀ = μ(k² = 0) = C + M + h²/(2M). While the Fermi point (7) exists also in other semiconductor systems that have been proposed to host chiral Majorana edge states [9–13], the Fermi points (8) and (9) are characteristic of a system with two different SOC terms.

In this section, we identify the topological phases of the non-trivial momentum space topology of Hamiltonian (1). The fundamental Hamiltonian (1) belongs to the Cartan–Altland–Zirnbauer symmetry class D, which in two spatial dimensions is characterized by a topological invariant $\in \mathbb {Z}$ [35]. The first Chern number which is commonly used as the topological invariant can be obtained for a non-interacting Hamiltonian with non-degenerate spectrum as the integral of the Berry curvature over the first Brillouin zone

$$\begin{equation} \mathcal{C}_1=\frac{1}{2\pi}\int_{\mathrm{BZ}} \mathcal{F}\,\rm\mathrm{d}^2k, \end{equation} \tag{ 10 }$$

where the Berry curvature $\mathcal {F}=\partial _{k_x}\mathcal {A}_{k_y}-\partial _{k_y}\mathcal {A}_{k_x}$ is obtained from the the Bloch functions via $\mathcal {A}_{k_\alpha }=-\mathrm {i}\sum _j\langle u _j(k)|\partial _{k_\alpha }|u_j(k)\rangle$, α = x,y, the index j running over the occupied states, as reviewed in [36].

We evaluate the general formula of the first Chern number (10) for the parameterized Hamiltonian (1) both numerically and analytically in the parameter regime outlined previously⁶.

Numerically, the first Chern number is computed from a discretized version [37] of (10) for a regularized version⁷ of Hamiltonian (1) using MATLAB at several illustrative values of the parameters on each side of the topological transitions. The topological transitions have been obtained by identifying the zero energy level crossings of the effective four-band Hamiltonian, see (7)–(9). This provides a good approximation for the topological transitions of the full model (1) in the considered parameter regime.

The analytical calculation [38] of the first Chern number is performed in the limit where an effective two-band description is possible [39] (see appendix).

We find that the first Chern number characterizing the topological phases of our physical system assumes the values 0,1,−1 in the considered parameter regime, depending on the specific values of the parameters in the following way:

$$\begin{equation} \mathcal{C}_1=\Theta{\left(B_E^2-\mu_0^2-\Delta_E^2 \right)} \mathrm{sign}\,{\left(R_0^2-\tilde{h}^2\right)}\mathrm{sign}\,{\left(B_E\right)}, \end{equation} \tag{ 11 }$$

where Θ(x) denotes the Heaviside step function.

The Chern numbers 1 and −1 indicate non-trivial topology, whereas the Chern number 0 indicates trivial topology which is adiabatically connected to a vacuum.

While the general form of (11) is analogous to known topological invariants of similar systems, such as spin–orbit coupled quantum wires [10, 12, 13] or 2D spin–orbit coupled semiconductors [9–11], this constitutes, to our knowledge, the first such calculation for a HgTe-QW in proximity to a superconductor. Importantly, this is also the first time that the BIA term h is included. When both Rashba-SOC and BIA are present, the sign of the topological invariant $\mathcal {C}_1$ can be reversed due to the interplay between both SOCs. This is similar to the combined effect of Rashba and Dresselhaus SOCs that was described for non-centrosymmetric superconductors [40].

In the absence of BIA, the relevant parameters for a topological transition to occur are the proximity coupling Δ_E, the Zeeman energy B_E and the sum of the mass gap and doping, μ₀ = C + M. When BIA is taken into account, the topological transition is shifted toward higher doping, as described by the term $\mu _0=C+M+\frac {h^2}{2M}$ (in the limit that M is the dominant energy scale and C ≈ −M). The topological transition is presented as a function of Zeeman energy B_E and doping C in figure 4.

Zoom In Zoom Out Reset image size

Under certain simplifying assumptions, the edge dispersion of the boundary problem can be calculated to first order in k using perturbation theory. We assume that it is sufficient to calculate the edge dispersion in the reduced model (5) and require that quadratic terms ∝k²_x,k²_y may be neglected (assumption 1). Consistently, we also match only the wave functions and not the derivatives at the boundary.

We first prove that the boundary conditions are fulfilled for E = 0, k_x = 0 and calculate the wave function of the corresponding bound state.

Imposing the condition E = k_x = 0 in (5) yields four possible values for λ = ik_y:

$$\begin{equation} \pm\lambda_\pm=\frac{\pm\left( |\Delta_E| \pm\sqrt{B_E^2-\mu_0^2}\right)}{\sqrt{R_0^2+\tilde{h}^2+2R_0\tilde{h}\sin{2\theta}}} . \end{equation} \tag{ 13 }$$

Depending on the parameter values, the signs of λ can be determined as follows under the assumption that |B_E| > |μ₀|:

The respective eigenvectors are

$$\begin{equation} |\lambda_+\rangle = \frac{1}{2}\left( \begin{array}{@{}c@{}} \mathrm{e}^{\mathrm{i}\chi}\sqrt{1-\mu_0/B_E}\\ -\mathrm{sign}\,{(B_E)}\sqrt{1+\mu_0/B_E}\\ \mathrm{e}^{\mathrm{i}\chi}\mathrm{sign}\,{(B_E)}\sqrt{1+\mu_0/B_E}\\ \sqrt{1-\mu_0/B_E} \end{array} \right) \end{equation} \tag{ 14 }$$

with $\mathrm {e}^{\mathrm {i}\chi }=(\mathrm {i}R_0+\tilde {h}\mathrm {e}^{2\theta \mathrm {i}})/\sqrt {R_0^2+\tilde {h}^2+2R_0\tilde {h}\sin {2\theta }}$ and where |λ₋〉, | − λ₊〉 and | − λ₋〉 can be obtained from |λ₊〉 by a clockwise translation of the minus sign in front of the second entry.

If B²_E > μ²_l + Δ²_E on the left side of the boundary (topologically non-trivial) and B²_E < μ²_r + Δ²_E on the right side of the boundary (topologically trivial), the matching condition requires that the eigenvectors on the left- and right-hand sides of the boundary are linearly dependent

$$\begin{equation} \det{\left(|\lambda^{\mathrm{l}}_+\rangle , |-\lambda^{\mathrm{l}}_-\rangle,|-\lambda^{\mathrm{r}}_+\rangle,|-\lambda^{\mathrm{r}}_-\rangle\right)} =0 \end{equation} \tag{ 15 }$$

yielding the solution

$$\begin{equation} \fl \psi(y)=\frac{1}{\mathcal{N}}\left(\Theta{(-y)} \mathrm{e}^{-\lambda^{\mathrm{l}}_-y}|-\lambda^{\mathrm{l}}_-\rangle +\Theta{(y)}\left(c_1 \mathrm{e}^{- \lambda^{\mathrm{r}}_+y}|-\lambda^{\mathrm{r}}_+\rangle + c_2\mathrm{e}^{- \lambda^{\mathrm{r}}_-y}|-\lambda^{\mathrm{r}}_-\rangle\right)\right). \end{equation} \tag{ 16 }$$

Here $\mathcal {N}$ is a normalization factor such that $\int \rm d y|\psi (y)|^2=1$. The coefficients c₁ and c₂ are determined by the matching condition.

This proves that there is a bound state for k_x = E = 0. The antisymmetry of the Hamiltonian (1) under particle–hole transformation Ξ^†HΞ = −H, where Ξ = iKσ_yτ_y, has the important consequence that ψ(y) for k_x = E = 0 constitutes a Majorana excitation, as can be easily confirmed by explicit calculation.

Now we consider the case of non-zero k_x. When k_x is small, the k_x-dependent BIA and Rashba-SOC terms in equation (5) can be treated as a perturbation V . The expectation value of this perturbation with respect to the k_x = E = 0 state gives the energy in terms of k_x to first order in perturbation theory

$$\begin{equation} E=\int {\rm d} y \psi^*(y)V\psi(y)=\frac{|\Delta_E|}{B_E}\frac{R_0^2-\frac{A^2h^2}{M^2}}{\sqrt{R_0^2+\frac{A^2h^2}{M^2}+2R_0\frac{Ah}{M}\sin{2\theta}}}k_x. \end{equation} \tag{ 17 }$$

The resulting edge energy dispersions that were calculated numerically from (1) and using (17) are presented and compared to the bulk energy dispersions in figure 6.

Zoom In Zoom Out Reset image size

The boundary between topologically trivial and non-trivial regions results in a single chiral MBS of approximately linear energy dispersion that smoothly connects to the occupied and non-occupied bulk bands. The existence of this edge state is independent of which boundary condition (soft or hard wall) or which inversion asymmetry term (BIA, Rashba-SOC or both) is considered. As expected, its existence depends on the topological bulk properties on both sides of the boundary—the edge state disappears when the Chern numbers on both sides of the boundary are identical.

The formula (17) describes the numerically obtained energy dispersions of figure 6 well, regardless of whether hard-wall or soft-wall boundary conditions have been employed. A small discrepancy between the numerical calculation on the basis of the Hamiltonian (1) and (17) is mostly due to the limited applicability of assumption 1. Yet, (17) describes the salient features of the edge dispersion well.

As predicted by (17), the edge state velocity depends most strongly on the parameters θ,A,h,R₀,B_E,Δ_E. We have checked that the parameters D,B,C,Δ_H,B_H have comparably little effect on the energy dispersion.

Interestingly, we find that the group velocity $v=\frac {1}{\hbar }\frac {\partial E}{\partial k_x}$ of the chiral edge state is tunable by the interplay between the Rashba-SOC and BIA spin–orbit interactions. The edge state is right-moving when the Rashba-SOC term dominates and left-moving when the BIA term dominates (see figure 7). This is consistent with the change of the first Chern number in (11) between $\mathcal {C}_1=\pm 1$ by means of the index theorem [32]. When applied to the hard-wall boundary condition, the index theorem equates the topological charge $\mathcal {C}_1$ of the bulk with the sign ν of the slope of the chiral edge mode. It is therefore easy to see that as the Rashba-SOC parameter R₀ changes from $R_0^2{<}\tilde {h}^2$ to $R_0^2{>}\tilde {h}^2$, the first Chern number changes from −1 to +1 and the edge state from left to right moving. This could be experimentally verified by changing the Rashba-SOC parameter, while the intrinsic BIA-SOC remains fixed. Another interesting observation is that the group velocity of the chiral Majorana edge state depends on the orientation of the edge with respect to the crystal axis when both Rashba-SOC and BIA-SOC terms are present. The absolute values of the velocity are smallest along the [110] (θ = π/4) and largest along the $[\bar {1}10]$ (θ = 3π/4) direction. In HgTe-QWs, the edge structures are, for experimental reasons, prepared in parallel to the cleavage planes; therefore, 95% of the current carrying edges are orientated either in the [110] or $[\bar {1}10]$ direction [41]. Due to the simultaneous presence of Rashba-SOC and BIA-SOC terms, these two experimentally accessible edge orientations are expected to have distinct chiral Majorana edge state properties.

Zoom In Zoom Out Reset image size

If the mass gap M is the dominant energy scale of the problem and C ≈ −M, leading to only four relevant E1-like energy bands, this allows us to reduce the eight-band Hamiltonian (1) to an effective four-band Hamiltonian.

To this end, we divide the Hamiltonian into three parts:

$$\begin{equation}\fl H=H_0+H_1+H_2, \end{equation} \tag{ A.1 }$$

$$\begin{equation} \fl H_0=\left(C+Ms_z\right)\tau_z, \end{equation} \tag{ A.2 }$$

$$\begin{equation}\fl H_1=R_0\frac{\left(s_z+1\right)}{2}(\hat{k}_x\sigma_y-\hat{k}_y\sigma_x)\tau_z +\left(\Delta_++\Delta_-s_z\right)\tau_x-\left(B+D\right)\hat{k}^2+\left(B_++B_-s_z\right)\sigma_z,\end{equation} \tag{ A.3 }$$

$$\begin{equation}\fl H_2=A\hat{k}_xs_x\sigma_z\tau_z-A\hat{k}_ys_y\tau_z +h\sigma_ys_y\tau_z , \end{equation} \tag{ A.4 }$$

where for clarity, in contrast to the more general Hamiltonian (1), we set θ = 0.

In the following, H₀ (containing only s_z) is considered as the unperturbed Hamiltonian and H₁ (containing only s_z) and H₂ (containing s_y) as the block-diagonal and block-offdiagonal perturbations, respectively. The perturbations are small in the sense that the mass gap M is the dominant energy scale, i.e.

$$\begin{equation} |\Delta_\pm|, \quad|B_\pm|, \quad|h|, \quad|A||k|, \quad|R_0||k|, \quad|B||k|^2, \quad|D||k|^2 \ll |M|. \end{equation} \tag{ A.5 }$$

In order to approximately diagonalize the Hamiltonian H

$$\begin{equation} \skew3\tilde{H}=\mathrm{e}^{-S}H\mathrm{e}^S \end{equation} \tag{ A.6 }$$

we apply quasi-degenerate perturbation theory ('Löwdin partitioning') [33].

The important matrix for the transformation, S, consists of different orders of the perturbation parameter, $S=S^{(1)}+\mathcal {O}(2)$. It can be shown that

$$\begin{equation} \skew3\tilde{H}=H_0+H_1+\textstyle\frac{1}{2}[H_2,S^{(1)}]+\mathcal{O}(3) \end{equation} \tag{ A.7 }$$

is block-diagonal up to the second order if S⁽¹⁾ is determined in such a way that

$$\begin{equation} [H_0,S^{(1)}]=-H_2. \end{equation} \tag{ A.8 }$$

For the Hamiltonian (A.1) it is straightforward to show⁸ that

$$\begin{equation}\fl S^{(1)}=+\mathrm{i}\frac{h}{2M}s_x\sigma_y -\mathrm{i}\frac{A\hat{k}_x}{2M}s_y\sigma_z-\mathrm{i}\frac{A\hat{k}_y}{2M}s_x, \end{equation} \tag{ A.9 }$$

$$\begin{eqnarray} &&\fl \tilde{H}=\mu(s_z,\hat{k}^2)\tau_z+B(s_z)\sigma_z+\Delta(s_z)\tau_x +\tilde{h}\hat{k}_x\sigma_x\tau_z-\tilde{h}\hat{k}_ys_z\sigma_y\tau_z +R(s_z)\hat{k}_x\sigma_y\tau_z-R(s_z)\hat{k}_y\sigma_x\tau_z,\nonumber\\ \end{eqnarray} \tag{ A.10 }$$

where the new parameters

$$\begin{equation} \mu(s_z,\hat{k}^2)=\epsilon(\hat{k}^2) +\left(M(\hat{k}^2)+\frac{A^2\hat{k}^2}{2M}+\frac{h^2}{2M} \right)s_z, \end{equation} \tag{ A.11 }$$

$$\begin{equation} B(s_z)=B_+ +B_-s_z, \end{equation} \tag{ A.12 }$$

$$\begin{equation} \Delta(s_z)=\Delta_++\Delta_-s_z, \end{equation} \tag{ A.13 }$$

$$\begin{equation} \tilde{h}=\frac{Ah}{M}, \end{equation} \tag{ A.14 }$$

$$\begin{equation} R(s_z)=R_0\frac{\left(s_z+1\right)}{2} \end{equation} \tag{ A.15 }$$

have been introduced. The advantage of the transformed Hamiltonian (A.10) is that s_z is (approximately) a good quantum number. For further calculations, we neglect the s_z = −1 solutions in the main text (5), and retrieve the angular dependence on θ by the simple replacement rules

$$\begin{equation} h\sigma_y\rightarrow h\left( \cos{2\theta}\sigma_y+\sin{2\theta}\sigma_x\right), \end{equation} \tag{ A.16 }$$

$$\begin{equation} h\sigma_x\rightarrow h\left( \cos{2\theta}\sigma_x-\sin{2\theta}\sigma_y\right) \end{equation} \tag{ A.17 }$$

which amount to a rotation in the xy-plane.

In order to verify (11) of the main text in an analytical way, we project the four-band Hamiltonian

$$\begin{equation} \tilde{H}=\left(\mu_0-\frac{k^2}{2m}\right)\tau_z+B_E\sigma_z+\Delta_E\tau_x+(R_0-\tilde{h})k_x\sigma_y\tau_z-(R_0+\tilde{h})k_y\sigma_x\tau_z \end{equation} \tag{ A.18 }$$

or, equivalently,

$$\begin{eqnarray}&&\fl \mathcal{H}=\int \mathrm{d}x\mathrm{d}y \left[\Psi^\dagger\left( \mu_0-\frac{k^2}{2m}+B_E\sigma_z +(R_0-\tilde{h})k_x\sigma_y-(R_0+\tilde{h})k_y\sigma_x\right)\Psi +(\Delta_E \psi_{\downarrow}\psi_\uparrow+\mathrm{h.c}.)\right]\nonumber\\ \end{eqnarray} \tag{ A.19 }$$

onto a two-band Hamiltonian. Since the precise value of the angle θ does not affect the topological properties of the bulk, we have chosen $\theta =\frac {\pi }{4}$ here for convenience. The projection procedure is similar to [39] and holds in the limit $B_E\gg m(R^2+\tilde {h}^2),|\Delta _E|$ where $m\!=\!\frac {A^2}{2M}\!-\!D\!-\!B$.

This projection relies on the approximations ψ_↓ ≈ ψ₋ and $\psi _\uparrow \approx \frac {(R_0+\tilde {h})k_y+\mathrm {i}(R_0-\tilde {h})k_x}{2B_E}\psi _-$ and yields

$$\begin{equation} \tilde{H}=\left(\mu_{\mathrm{eff}}-\frac{k^2}{2m}\right)s_z+\Delta^y_{\mathrm{eff}}k_ys_x+\Delta^x_{\mathrm{eff}}k_xs_y, \end{equation} \tag{ A.20 }$$

where μ_eff = μ − B_E, $\Delta ^x_{\mathrm {eff}}=\frac {|\Delta _E|}{2|B_E|}(R_0-\tilde {h})$ and $\Delta ^y_{\mathrm {eff}}=\frac {|\Delta _E|}{2|B_E|}(R_0+\tilde {h})$.

After regularization, this can be written as

$$\begin{equation} \tilde{H}=\sum_{i=1}^3 h_is_i, \end{equation} \tag{ A.21 }$$

$$\begin{equation} h_1=-\frac{1}{a}\Delta^y_{\mathrm{eff}}\sin{ak_y},\end{equation} \tag{ A.22 }$$

$$\begin{equation} h_2=\frac{1}{a}\Delta^x_{\mathrm{eff}}\sin{ak_x},\end{equation} \tag{ A.23 }$$

$$\begin{equation} h_3=\mu_{\mathrm{eff}}-\frac{1}{2a^2m}(\cos{ak_x}+\cos{ak_y}-2). \end{equation} \tag{ A.24 }$$

For this two-band model, the first Chern number is defined as [1]

$$\begin{equation} \mathcal{C}_1=\frac{1}{2\pi}\int_{\rm BZ} \mathcal{F}\mathrm{d}^2k =\frac{1}{4\pi}\int_{\rm BZ} {\hat{\bf h}}\cdot \left(\frac{\partial\hat{{\bf h}}}{\partial k_x}\times \frac{\partial\hat{{\bf h}}}{\partial k_y}\right)\mathrm{d}^2k, \end{equation} \tag{ A.25 }$$

where $\hat {{\bf{ h}}}$ is defined as $\hat {{\bf{ h}}}={\bf{ h}}/|{\bf{ h}}|$.

The explicit integration is cumbersome and can be greatly simplified using the equivalent formula [38] which consists of a sum over the 'Dirac points' $D_3=\{(0,0),(0,\frac {\pi }{a}),(\frac {\pi }{a},0),(\frac {\pi }{a},\frac {\pi }{a})\}$ of the two-Pauli-matrix Hamiltonian h₁s₁ + h₂s₂,

$$\begin{eqnarray} \fl \mathcal{C}_1&=&\frac{1}{2}\sum_{k\in D_3} \mathrm{sign}\,{(\partial_{k_x}{\bf h}\times \partial_{k_y}{\bf h})_3\mathrm{sign}\,(h_3)}\nonumber\\ \fl&=&-\frac{1}{2}\sum_{k\in D_3}\mathrm{sign}\,{\left((R_0^2-\tilde{h}^2)\cos{ak_x}\cos{ak_y}\right)} \mathrm{sign}\,\left(\mu_{\mathrm{eff}}-\frac{1}{2a^2m}(\cos{ak_x}+\cos{ak_y}-2)\right)\nonumber\\&&\fl \end{eqnarray} \tag{ A.26 }$$

$$\begin{eqnarray}\fl \hphantom{\mathcal{C}_{1}} & =&-\frac{1}{2}\mathrm{sign}\,{(R_0^2-\tilde{h}^2)}\left[\mathrm{sign}\,{\mu_{\mathrm{eff}}}-2\mathrm{sign}\,{\left(\mu_{\mathrm{eff}}+\frac{1}{a^2m}\right)}+\mathrm{sign}\,{\left(\mu_{\mathrm{eff}}+\frac{2}{a^2m}\right)}\right]\nonumber\\ \fl \hphantom{\mathcal{C}_{1}} &=&\mathrm{sign}\,{(R_{0}^2-(Ah/M)^2)}\mathrm{sign}\,{\left(\frac{1}{a^2m}+\mu_{\mathrm{eff}}\right)}\Theta{(-\mu_{\mathrm{eff}})}\Theta{\left(-\frac{2}{a^2m}-\mu_{\mathrm{eff}}\right).}\nonumber\\&&\fl \end{eqnarray} \tag{ A.27 }$$

If we assume that the regularization is such that $\frac {1}{a^2m}\gg |\mu _{\mathrm {eff}}|$, then the Chern number simplifies to

$$\begin{equation} \mathcal{C}_1=\mathrm{sign}\,{(R_{0}^2-(Ah/M)^2)}\Theta{(B_E-\mu_0)} \end{equation} \tag{ A.28 }$$

which is for B_E > 0 consistent with formula (11) in the main text.

In the case B_E < 0, the projection relies on the approximations ψ_↑ ≈ ψ₋ and $\psi _\downarrow \approx \frac {-(R_0+\tilde {h})k_y-\mathrm {i}(R_0-\tilde {h})k_x}{2B_E}\psi _-$ and yields

$$\begin{equation} \mathcal{C}_1=-\mathrm{sign}\,{(R_{0}^2-(Ah/M)^2)}\Theta{(-B_E-\mu_0)} \end{equation} \tag{ A.29 }$$

which is for B_E < 0 consistent with formula (11) in the main text.

In summary, we have verified

$$\begin{equation} \mathcal{C}_1=\mathrm{sign}\,(B_E)\mathrm{sign}\,{(R_{0}^2-(Ah/M)^2)}\Theta{(B_E^2-\mu_0^2)} \end{equation} \tag{ A.30 }$$

and therefore equation (11) in the main text.