Supersymmetry dictated topology in periodic gauge fields and realization in strained and twisted 2D materials

In this work we consider SUSY Hamiltonian of the form

$$\begin{aligned} \mathcal{H}_{\text{susy}}=\begin{pmatrix}\mathbb{H}_{+}&0\\ 0&\mathbb{H}_{-}\end{pmatrix}=\begin{pmatrix}MM^{\dagger}&0\\ 0&M^{\dagger}M\end{pmatrix}\end{aligned},$$

(1)

where $M$ denotes a generic matrix. Its connection to the SUSY quantum mechanics Junker (2019) is revealed by rewriting

$$\mathcal{H}_{\text{susy}}=\{Q,\,Q^{\dagger}\},~{}$$

(2)

where

$\displaystyle Q=\begin{pmatrix}0&M\\ 0&0\end{pmatrix}~{}~{}\text{and}~{}~{}Q^{\dagger}=\begin{pmatrix}0&0\\ M^{\dagger}&0\end{pmatrix}$

(3)

are the complex supercharges (or SUSY generators). One can readily verify that

$$\{Q,\,Q\}=\{Q^{\dagger},\,Q^{\dagger}\}=[\mathcal{H}_{\text{susy}},\,Q]=[\mathcal{H}_{\text{susy}},\,Q^{\dagger}]=0.$$

(4)

These relations and Eq. (2) define the superalgebra and a SUSY quantum system Junker (2019). Such a SUSY system contains two parts ($\mathbb{H}_{\pm}$), and particles characterized by $\mathbb{H}_{+}\psi_{+}=\mathcal{E}_{+}\psi_{+}$ and $\mathbb{H}_{-}\psi_{-}=\mathcal{E}_{-}\psi_{-}$ are considered as SPs of each other. The eigenvectors of $\mathcal{H}_{\text{susy}}$ can be expressed in terms of $\psi_{\pm}$ as $\Psi_{+}=(\psi_{+},\,0)^{T}$ and $\Psi_{-}=(0,\,\psi_{-})^{T}$.

A SUSY system has several characteristic properties. First, as $\mathcal{H}_{\text{susy}}=\{Q,\,Q^{\dagger}\}$, its eigenenergies are non-negative. Second, the relation $[\mathcal{H}_{\text{susy}},\,Q]=[\mathcal{H}_{\text{susy}},\,Q^{\dagger}]=0$ dictates that any positive eigenvalue, i.e. $\mathcal{E}>0$, must have a pair of eigenstates from the two SPs respectively, satisfying $\mathcal{H}_{\text{susy}}\Psi_{\pm}=\mathcal{E}\Psi_{\pm}$. They can transform into each other via the supercharges $\Psi_{+}=Q\Psi_{-}/\sqrt{\mathcal{E}}$ and $\Psi_{-}=Q^{\dagger}\Psi_{+}/\sqrt{\mathcal{E}}$, or equivalently,

$$\psi_{+}=\frac{1}{\sqrt{\mathcal{E}}}M\psi_{-},~{}~{}\psi_{-}=\frac{1}{\sqrt{\mathcal{E}}}M^{\dagger}\psi_{+}.$$

(5)

The zero (energy) modes, $\Psi_{+}^{0}=(\psi_{+}^{0},\,0)^{T}$ and $\Psi_{-}^{0}=(0,\,\psi_{-}^{0})^{T}$, of the SUSY system are the ones annihilated by the supercharges, $M^{\dagger}\psi_{+}^{0}=0$ and $M\psi_{-}^{0}=0$. Thus the transformation connecting the SPs in Eq. (5) fails for the zero modes. This indicates that the zero-energy modes may remain unpaired or even be absent com (b). The difference in the number of zero modes between the SPs is the Witten index, which can be employed to characterize the SUSY systems.

A nontrivial scenario where SUSY emerges corresponds to 2D electrons in a gauge field. We consider $M$ and $M^{\dagger}$ given in terms of $(\pi_{x},\,\pi_{y})=\boldsymbol{p}+e\boldsymbol{A}$, where $\boldsymbol{p}$ is the momentum operator and $\boldsymbol{A}$ is the gauge potential:

$$\begin{aligned} M=\alpha_{N}(\pi_{x}-i\pi_{y})^{N}=\alpha_{N}\pi_{-}^{N}\\ M^{\dagger}=\alpha_{N}(\pi_{x}+i\pi_{y})^{N}=\alpha_{N}\pi_{+}^{N}\end{aligned}.$$

(6)

Here $N$ denotes a positive integer, $\alpha_{N}$ is a constant such that $MM^{\dagger}$ and $M^{\dagger}M$ have units of energy, and we have also defined $\pi_{\pm}=\pi_{x}\pm i\pi_{y}$. The SUSY Hamiltonian is then

$$\begin{aligned} \mathcal{H}_{\text{susy}}=\begin{pmatrix}\mathbb{H}_{+}&0\\ 0&\mathbb{H}_{-}\end{pmatrix}=\begin{pmatrix}\alpha_{N}^{2}\pi_{-}^{N}\pi_{+}^{N}&0\\ 0&\alpha_{N}^{2}\pi_{+}^{N}\pi_{-}^{N}\end{pmatrix}\end{aligned}.~{}$$

(7)

The out-of-plane gauge field $\boldsymbol{B}=B\hat{\boldsymbol{z}}$ satisfies $[\pi_{x},\,\pi_{y}]=-ie\hbar B$. We will adhere to the above conventions for assigning $\mathbb{H}_{+}$ and $\mathbb{H}_{-}$ throughout the manuscript. In the following we will assume that the field is Abelian, the discussions will be extended to cover non-Abelian gauge fields later.

In the case of $N=1$, the model describes spin–$\frac{1}{2}$ 2D electrons in an out-of-plane magnetic field,

$$M=\frac{1}{\sqrt{2m}}\pi_{-},~{}~{}M^{\dagger}=\frac{1}{\sqrt{2m}}\pi_{+}~{}~{}(\text{when}~{}N=1),$$

(8)

and

$$\begin{aligned} \mathbb{H}_{\pm}^{(N=1)}=\frac{1}{2m}\pi_{\mp}\pi_{\pm}=\frac{1}{2m}(\boldsymbol{p}+e\boldsymbol{A})^{2}\pm\boldsymbol{\mu}_{B}\cdot\boldsymbol{B}\end{aligned},$$

(9)

where $\boldsymbol{\mu}_{B}=\frac{e\hbar}{2m}\hat{\boldsymbol{z}}$ is the Bohr magneton. Spin-up and spin-down electrons are SPs of each other in this case.

Information of the zero modes, i.e., their degeneracy and wave functions, can be obtained by solving $\pi_{\pm}\psi_{\pm}^{0}=0$, which depend on details of $\boldsymbol{B}$. Here we distinguish two scenarios:
(i) The situation of $\boldsymbol{B}$ with nonzero flux is known: the zero modes are spin-polarized with its degeneracy determined by the flux [see e.g., Fig. 3(a)], and notably, these features persist even when $\boldsymbol{B}$ is spatially nonuniform Aharonov and Casher (1979).
(ii) The scenario of periodic $\boldsymbol{B}$ with zero flux exhibits unique characteristics that have not been uncovered, and we will focus on them in this work. In this particular case, the two SPs each contains one unique zero mode, $\psi_{\pm}^{0}=c_{\pm}e^{\pm\phi}$, where $c_{\pm}$ are normalization constants, and the Coulomb gauge has been used with $\phi$ satisfying $e\boldsymbol{A}/\hbar=-\nabla\times(\phi\hat{\boldsymbol{z}})$ and $\partial_{\boldsymbol{r}}^{2}\phi=eB/\hbar$. Note that, although here we start with the case of $N=1$, the existence of such zero modes remains valid when $N>1$. Similar zero modes have also been found in graphene with a periodic gauge field Snyman (2009); Phong and Mele (2023) (also see Sec. V). Therefore, the spectra of two SPs are exactly the same when $\boldsymbol{B}$ is zero-flux periodic. This isospectral behavior is confirmed by solving Eq. (9) in a zero-flux $\boldsymbol{B}$ field of periodicity $L$ [Fig. 1(a)]: $\boldsymbol{B}(\boldsymbol{r})=\hat{z}B_{0}\sum_{i=1}^{3}\cos\left(\boldsymbol{b}_{i}\cdot\boldsymbol{r}\right)$, $\boldsymbol{b}_{1}=(\frac{\sqrt{3}}{2},\,-\frac{1}{2})\frac{4\pi}{\sqrt{3}L}$ and $\boldsymbol{b}_{i}=C_{3}^{i-1}\boldsymbol{b}_{1}$, where $C_{3}$ denotes the three-fold rotation operation. As expected, $\mathcal{E}_{+}$ and $\mathcal{E}_{-}$ each has one zero mode located at the $\gamma$ point [Fig. 1(b)].

We will consider TMDs as an example below, but the discussions are also applicable for many other 2D semiconductors, e.g., graphene with a staggered potential or hBN. In monolayer TMDs, the low-energy carriers are governed by massive Dirac equations $h^{\text{1L}}_{\tau}=v\boldsymbol{\sigma}_{\tau}\cdot\boldsymbol{p}+\Delta\sigma_{z}$ in two inequivalent valleys labeled as $\tau K$, where $\tau=\pm$ is the valley index, $v$ is the Fermi velocity, and $\boldsymbol{\sigma}_{\tau}=(\tau\sigma_{x},\,\sigma_{y})$ Xiao et al. (2012). By introducing nonuniform strain patterns, giant pseudo-magnetic fields opposite in the two valleys can be generated Vozmediano et al. (2010); Amorim et al. (2016); Naumis et al. (2017). The Dirac Hamiltonian in the presence of strain pseudo-magnetic field becomes

$$H^{\text{D}}_{\tau}=v\boldsymbol{\sigma}_{\tau}\cdot(\boldsymbol{p}+\tau e\boldsymbol{A})+\Delta\sigma_{z},~{}$$

(17)

where $\tau\boldsymbol{A}$ is the strain vector potential and its details will be discussed later. Note that there might also exist strain-induced scalar potentials. We discuss such terms in Sec. S3 of Supplementary Material sup , and show that WSe₂ is an ideal candidate, in which such terms can be neglected. Importantly, profile of the pseudo-magnetic field, $\tau\boldsymbol{B}=\nabla\times(\tau\boldsymbol{A})$, can be controlled by tuning the strain distribution. Strain engineering in 2D materials Guinea et al. (2010); Tomori et al. (2011); Reserbat-Plantey et al. (2014); Jiang et al. (2017); Zhang et al. (2018a); Phong and Mele (2022); Mahmud et al. (2023) thus offers a solution to overcome the obstacle of generating patterned magnetic fields with zero flux.

The Dirac model $H^{\text{D}}_{\tau}$ describes charge carriers in both the conduction and valence bands. From experimental perspectives, for TMDs it is most natural to consider holes (h) near the valence band edge, which are described by

$\displaystyle H^{h}_{\tau}\psi=[\frac{1}{2m}(\boldsymbol{p}+\tau e\boldsymbol{A})^{2}+\boldsymbol{\mu}_{\tau}^{h}\cdot(\tau\boldsymbol{B})]\psi=\varepsilon\psi$

(18)

near the valence band edge in the limit of $\varepsilon\ll\Delta$, where $\boldsymbol{\mu}_{\tau}^{h}=-\tau\boldsymbol{\mu}_{B}$ is the valley magnetic moment Xiao et al. (2007) for holes. In this quadratic model, $\varepsilon$ is measured from the valence band edge and $m=\Delta/v^{2}$ is the effective mass. Note that a magnetic moment in the form of Bohr magneton naturally appears here, which is a characteristic of the low-energy behavior of the massive Dirac model. More importantly, the strain pseudo-magnetic field can only couple to this valley magnetic moment, but not other sources (e.g., spin and orbital Aivazian et al. (2015)). Therefore, the second challenge mentioned above for realizing SUSY is resolved.

In this section, we propose a physical system for the scenario of $N=2$ SUSY in Eq. (7). It allows us to obtain $C_{1}^{-}-C_{1}^{+}=2$ thanks to the quartic dispersion around the zero mode [Fig. 2(a)]. Such cases are useful for achieving topological bands with higher Chern numbers. The system is realized with the replacement of each monolayer in Fig. 5(a) by a bilayer, as shown schematically in Fig. 6(a). Notably, the isospectral SPs realized are still in stark contrast to the time-reversal partners, with one SP being topologically trivial ($C_{1}^{-}=0$), and another nontrivial ($C_{1}^{+}=-2$).

We first look at the basic ingredient of the setup– a strained bilayer. For illustrative purposes, it suffices to consider a simple bilayer Hamiltonian

$$\mathcal{H}_{\tau}^{\text{2L}}=\begin{pmatrix}H_{\tau}^{\text{D}}&T\\ T^{\dagger}&H_{\tau}^{\text{D}}\end{pmatrix}~{}\text{ with }~{}T=\begin{pmatrix}0&0\\ t&0\end{pmatrix},$$

(21)

where $H_{\tau}^{\text{D}}=v\boldsymbol{\sigma}_{\tau}\cdot(\boldsymbol{p}+\tau e\boldsymbol{A})+\Delta\sigma_{z}$ is the Dirac Hamiltonian for a strained monolayer, and $t$ is the interlayer coupling constant. This simple model can describe rhombohedral bilayer TMDs Zhang et al. (2018b); Tong et al. (2017), where $t\sim 0.1$ eV from a recent experiment Liang et al. (2022). By focusing on holes, one can consider a low-energy model

$\displaystyle\mathcal{H}_{\tau}^{h}\psi=\frac{v^{4}}{2\Delta t^{2}}(\tau\pi_{x}^{\tau}+i\pi_{y}^{\tau})^{2}(\tau\pi_{x}^{\tau}-i\pi_{y}^{\tau})^{2}\psi=\varepsilon\psi,$

(22)

where $(\pi_{x}^{\tau},\,\pi_{y}^{\tau})=\boldsymbol{\pi}^{\tau}=\boldsymbol{p}+\tau e\boldsymbol{A}$. This equation is valid for $\varepsilon\ll t,\,\Delta$. One may notice that $\mathcal{H}_{\tau}^{h}$ can also describe holes in Bernal bilayer graphene by introducing an interlayer displacement field into the two-band model McCann and Koshino (2013). In this case, $\Delta$ is replaced by the displacement field, and $t\sim 0.4$ eV. The Hamiltonian $\mathcal{H}_{\tau}^{h}$ resembles that in Eq. (7) with $N=2$, thus lays the basis for constructing the SUSY model.

By following a similar idea as in Fig. 5(a), the SUSY model in Eq. (7) with $N=2$ can be realized by using two copies of bilayer TMDs with layer-opposite strain pseudo-magnetic fields (i.e., $\boldsymbol{B}^{(1)}=-\boldsymbol{B}^{(2)}$). Then holes from opposite valleys of opposite bilayers form SPs, for example,

$$\begin{aligned} \mathcal{H}_{\tau=+}^{h}(\text{bilayer}~{}1)&=\frac{v^{4}}{2\Delta t^{2}}[\pi_{+}(\boldsymbol{A}^{(1)})]^{2}[\pi_{-}(\boldsymbol{A}^{(1)})]^{2}\\ \mathcal{H}_{\tau=-}^{h}(\text{bilayer}~{}2)&=\frac{v^{4}}{2\Delta t^{2}}[\pi_{-}(\boldsymbol{A}^{(1)})]^{2}[\pi_{+}(\boldsymbol{A}^{(1)})]^{2}\end{aligned},~{}$$

(23)

where the notation $\pi_{\pm}(\boldsymbol{A}^{(1)})=(p_{x}+eA_{x}^{(1)})\pm i(p_{y}+eA_{y}^{(1)})$ has been used, and $-\boldsymbol{A}^{(2)}$ has been replaced by $\boldsymbol{A}^{(1)}$ in layer 2. It now becomes clear that they correspond to $\mathbb{H}_{\mp}$ in Eq. (7) with $N=2$. In this form it is also evident that the low-energy dispersion will be quartic, thus one expects $C_{1}^{-}-C_{1}^{+}=2$ if there is one zero mode on the first band [Fig. 2(a)].

By depositing double bilayers with $\theta=15^{\circ}$ and $\theta=45^{\circ}$ respectively on the substrate [Fig. 6(a)], opposite pseudo-magnetic fields in the two bilayers can be achieved [see Fig. 4(b)], thus realizing the model of Eq. (23). Figure 6(c) shows representative energy bands for such SPs by considering double bilayers of WSe₂. The low-energy bands are obtained from Eq. (21) due to its simplicity. The band structure around the zero mode at $\gamma$ clearly shows a quartic dispersion [also compare with Fig. 5(c)]. The first bands of the SPs exhibit distinct band topology characterized by $C_{1}^{+}=-2$ vs $C_{1}^{-}=0$, which is consistent with $C_{1}^{-}-C_{1}^{+}=2$ due to the quartic dispersion around $\gamma$ [Eq. (15) and Fig. 2(a)].

One notices that the lowest bands of strained TMDs, e.g., Fig. 5(c), are much narrower than their counterparts in graphene Milovanović et al. (2020); Manesco et al. (2021); Manesco and Lado (2021); Phong and Mele (2022); Mahmud et al. (2023); Gao et al. (2023). Hence Coulomb interaction can play important roles when these bands are fractionally filled. Despite having identical energy dispersion, the SPs possess drastically different wave functions [see e.g., Fig. 1(c) vs Fig. 1(d)] that underlie their distinct topological and quantum geometrical properties, hence spontaneous SUSY breaking by Coulomb interaction can be anticipated.

Taking the system of Fig. 5(c) as an example, we project Coulomb interaction onto the lowest band of each SP and perform exact diagonalization calculations. The adopted Coulomb interaction before projection reads

$$H_{\text{int}}=\frac{1}{2S}\sum_{\boldsymbol{k},\,\boldsymbol{k}^{\prime},\,\boldsymbol{q}}V(\boldsymbol{q})c_{\boldsymbol{k}+\boldsymbol{q}}^{\dagger}c_{\boldsymbol{k}^{\prime}-\boldsymbol{q}}^{\dagger}c_{\boldsymbol{k}^{\prime}}c_{\boldsymbol{k}},$$

(24)

where $S$ is the system’s area, $c_{\boldsymbol{k}}^{\dagger}$ creates a plane wave with momentum $\boldsymbol{k}$, $V(\boldsymbol{q})=e^{2}\tanh{(qd)}/(2\epsilon_{0}\epsilon_{r}q)$ is the dual-gate screened Coulomb potential, $d$ is the distance from the material to the metallic gates, $\epsilon_{0}$ is the vacuum permittivity, and $\epsilon_{r}$ is the relative dielectric constant. With their distinct wave functions [see Fig. 1(c, d)], projection of the Coulomb potential to the lowest band leads to different form factors of the projected interactions for the twin SPs.

Figure 7(a) and (b) present the many-body spectra as a function of crystal momentum obtained at filling factor $\nu=2/3$ for the two isospectral SPs with topologically nontrivial and trivial lowest band, respectively, with exact diagonalization calculations performed on a $3\times 9$ unit-cell cluster put on a torus. The crystal momentum is assigned as $k_{1}/N_{1}\boldsymbol{b}_{1}+k_{2}/N_{2}\boldsymbol{b}_{2}$ ($N_{1}=3$ and $N_{2}=9$ here) and labeled by an integer index $k_{1}+N_{1}k_{2}$. Three nearly degenerate ground states can be identified in Fig. 7(a), which are separated from excited states by an energy gap of $\sim 3$ meV. A similar feature of ground state degeneracy is also observed from the calculation with a different cluster size of $4\times 6$ unit cells in Fig. 7(c). This approximate ground state degeneracy is consistent with that of a fractional quantum Hall state on a torus Wen and Niu (1990). Additionally, upon magnetic flux insertion, these three ground states evolve into each other after $2\pi$ flux and return to themselves with a period of $6\pi$ [Fig. 7(d)]. These features provide strong evidence that the ground state at $\nu=2/3$ filling of this SP is FQAH state in nature. In stark contrast, the many-body spectrum at the same filling for its twin SP with topologically trivial lowest band is completely distinct [e.g., Fig. 7(b) vs Fig. 7(a)]. The spectrum in Fig. 7(b) shows the opening of a larger correlation gap, and has the character of the CDW phase, where the three degenerate ground states have momenta that can be folded back to the $\gamma$ point of a tripled unit cell Reddy et al. (2023); Yu et al. (2024b). The spontaneous SUSY breaking thus takes the form of developing distinct correlated phases for the two SPs with identical single-particle spectrum.

The above results suggest that the SUSY systems are promising platforms for exploring correlation-driven effects that are sensitive to topology and quantum geometry, which are being intensively explored with the groundbreaking experimental discovery of FQAH effects in twisted bilayer MoTe₂ Park et al. (2023); Cai et al. (2023); Zeng et al. (2023); Xu et al. (2023) and in pentalayer graphene/hBN moiré system Lu et al. (2024).

Here we define the square root of $\mathcal{H}_{\text{susy}}$ in Eq. (1) as

$$H_{\text{sqrt}}=\begin{pmatrix}\Delta\mathbb{I}&M\\ M^{\dagger}&-\Delta\mathbb{I}\end{pmatrix},$$

(25)

where $\Delta$ is a constant and $\mathbb{I}$ is the identity matrix com (e). In the following we will assume $\Delta\geq 0$, and the discussions are not limited to the forms of $M$ in Eq. (6). $H_{\text{sqrt}}$ can be considered as the square root of $\mathcal{H}_{\text{susy}}$ in the sense that its square yields $\mathcal{H}_{\text{susy}}$ up to a constant:

$$H_{\text{sqrt}}^{2}=\begin{pmatrix}MM^{\dagger}&0\\ 0&M^{\dagger}M\end{pmatrix}+\Delta^{2}=\mathcal{H}_{\text{susy}}+\Delta^{2},$$

(26)

where the identity matrix has been removed for simplicity. As the SP Hamiltonian $\mathbb{H}_{+}=MM^{\dagger}$ and $\mathbb{H}_{-}=M^{\dagger}M$ have identical positive energies, which will be denoted as $\mathcal{E}>0$, Eq. (26) dictates that the spectrum of $H_{\text{sqrt}}$ consists of $E^{\text{sqrt}}=\sqrt{\mathcal{E}+\Delta^{2}}>\Delta$, $E^{\text{sqrt}}=-\sqrt{\mathcal{E}+\Delta^{2}}<-\Delta$, and $E^{\text{sqrt}}=\pm\Delta$ if $\mathbb{H}_{\pm}$ have zero modes. This points to a symmetric spectrum for $H_{\text{sqrt}}$ with the possible exception at energies $\pm\Delta$, where the eigenmodes are also dubbed as threshold modes.

When $\Delta=0$, the square-root Hamiltonian $H_{\text{sqrt}}$ has chiral symmetry, i.e., it satisfies $\{H_{\text{sqrt}},\,\Sigma_{z}\}=0$ with $\Sigma_{z}=\text{diag}(\mathbb{I},\,-\mathbb{I})$. Many lattice models consider such a case Ezawa (2020); Mizoguchi et al. (2020); Yoshida et al. (2021); Roychowdhury et al. (2022). Chiral symmetry ensures a symmetric spectrum $E\leftrightarrow-E$, however, it should be noted that symmetric spectrum can also appear from SUSY in the absence of chiral symmetry when $\Delta\neq 0$.

When $\Delta\neq 0$, the relation $\{H_{\text{sqrt}},\,[H_{\text{sqrt}},\,\Sigma_{z}]\}=0$ is satisfied, thus $[H_{\text{sqrt}},\,\Sigma_{z}]\propto Q^{\dagger}-Q$ [see Eq. (3)] transforms eigenstates of $H_{\text{sqrt}}$ with opposite energies into each other Kailasvuori (2009). An exception occurs for the threshold modes: $Q$ and $Q^{\dagger}$ annihilate the threshold modes, $\Psi_{\Delta}^{\text{sqrt}}=(\psi_{+}^{0},\,0)^{T}$ and $\Psi_{-\Delta}^{\text{sqrt}}=(0,\,\psi_{-}^{0})^{T}$, which are also the zero modes of $\mathcal{H}_{\text{susy}}$ satisfying $M^{\dagger}\psi_{+}^{0}=0$ and $M\psi_{-}^{0}=0$. Thus the energy levels at $E^{\text{sqrt}}=\Delta$ and $E^{\text{sqrt}}=-\Delta$ are not necessarily paired. However, intriguing features can emerge when $H_{\text{sqrt}}$ has symmetric threshold modes at both $\Delta$ and $-\Delta$, or equivalently, when $\mathbb{H}_{+}$ and $\mathbb{H}_{-}$ both have zero modes, which will be the focus now.

In the following we consider periodic systems. Bloch functions of the periodic $H_{\text{sqrt}}$ are also related to those of $\mathbb{H}_{\pm}$. We denote the eigenvalue of the Bloch Hamiltonian $H_{\text{sqrt}}(\boldsymbol{k})$ as $E_{n\boldsymbol{k}}^{\text{sqrt}}$, with $n=\pm 1,\pm 2,\cdots$ labeling conduction ($n>0$) and valence ($n<0$) bands in ascending order of energy magnitude. The eigenvalues of $\mathbb{H}_{\pm}(\boldsymbol{k})$, which are non-negative, will be denoted as $\mathcal{E}_{|n|\boldsymbol{k}}^{\pm}$ with $|n|=1,2,\cdots$, and they satisfy

$$\mathcal{E}_{|n|\boldsymbol{k}}^{\pm}=\left(E_{n\boldsymbol{k}}^{\text{sqrt}}\right)^{2}-\Delta^{2}\geq 0.$$

(27)

Given the SP Bloch states $u_{|n|\boldsymbol{k}}^{\pm}$ satisfying $\mathbb{H}_{\pm}(\boldsymbol{k})u_{|n|\boldsymbol{k}}^{\pm}=\mathcal{E}_{|n|\boldsymbol{k}}^{\pm}u_{|n|\boldsymbol{k}}^{\pm}$, one can obtain the normalized Bloch state of $H_{\text{sqrt}}(\boldsymbol{k})$ in terms of $u_{|n|\boldsymbol{k}}^{\pm}$ as

$$u_{n\boldsymbol{k}}^{\text{sqrt}}\equiv\begin{pmatrix}c_{n\boldsymbol{k}}^{+}u_{|n|\boldsymbol{k}}^{+}\\ c_{n\boldsymbol{k}}^{-}u_{|n|\boldsymbol{k}}^{-}\end{pmatrix}=\begin{pmatrix}\sqrt{\frac{E_{n\boldsymbol{k}}^{\text{sqrt}}+\Delta}{2E_{n\boldsymbol{k}}^{\text{sqrt}}}}u_{|n|\boldsymbol{k}}^{+}\\ \text{sgn}(E_{n\boldsymbol{k}}^{\text{sqrt}})\sqrt{\frac{E_{n\boldsymbol{k}}^{\text{sqrt}}-\Delta}{2E_{n\boldsymbol{k}}^{\text{sqrt}}}}u_{|n|\boldsymbol{k}}^{-}\end{pmatrix}.$$

(28)

In general, the two-component $u_{n\boldsymbol{k}}^{\text{sqrt}}$ consists of both SP Bloch functions $u_{|n|\boldsymbol{k}}^{\pm}$, except for the threshold modes at energy $E_{n\boldsymbol{k}}^{\text{sqrt}}=\Delta$ ($-\Delta$), which are fully polarized on $u_{|n|\boldsymbol{k}}^{+}$ ($u_{|n|\boldsymbol{k}}^{-}$) that correspond to the zero modes of $\mathbb{H}_{+}$ ($\mathbb{H}_{-}$).

When the bands are narrow, satisfying $|E_{n\boldsymbol{k}}^{\text{sqrt}}|-\Delta\ll\Delta$, one finds $c_{n\boldsymbol{k}}^{+}\approx 1$ & $c_{n\boldsymbol{k}}^{-}\approx 0$ ($c_{n\boldsymbol{k}}^{+}\approx 0$ & $c_{n\boldsymbol{k}}^{-}\approx-1$) in the conduction (valence) bands, thus $u_{n\boldsymbol{k}}^{\text{sqrt}}$ is near fully polarized on the $u_{|n|\boldsymbol{k}}^{+}$ ($u_{|n|\boldsymbol{k}}^{-}$) component. In such cases, electrons in the conduction and valence bands of $H_{\text{sqrt}}$ can be well described by $H_{\text{sqrt}}^{c}=MM^{\dagger}/(2\Delta)$ and $H_{\text{sqrt}}^{v}=-M^{\dagger}M/(2\Delta)$, respectively.

The mapping from energy $\mathcal{E}_{|n|\boldsymbol{k}}^{\pm}$ and states $u_{|n|\boldsymbol{k}}^{\pm}$ of $\mathbb{H}_{\pm}$ to those of $H_{\text{sqrt}}$ are shown schematically in Fig. 8(a). In the right column, the ${E}_{n\boldsymbol{k}}^{\text{sqrt}}$ dispersions consist of dots whose red and blue portions represent the $u_{|n|\boldsymbol{k}}^{+}$ and $u_{|n|\boldsymbol{k}}^{-}$ components of $u_{n\boldsymbol{k}}^{\text{sqrt}}$ respectively [see Eq. (28)]. Reciprocally, one can deduce the eigenenergies and eigenstates of $\mathbb{H}_{\pm}$ from $E_{n\boldsymbol{k}}^{\text{sqrt}}$ and $u_{n\boldsymbol{k}}^{\text{sqrt}}$. And usually diagonalizing $H_{\text{sqrt}}$ can be done more conveniently as $H_{\text{sqrt}}$ is linear in $M$ and $M^{\dagger}$.

Eq. (28) further implies that the band topology of $H_{\text{sqrt}}$ are related to those of $\mathbb{H}_{\pm}$. Recall that the states $u_{|n|\boldsymbol{k}}^{\pm}$ are paired by Eq. (5) when the energy is nonzero, this allows one to relate the Berry connection of $H_{\text{sqrt}}$ to those of $\mathbb{H}_{\pm}$, whose details are given in Sec. S5 of Supplementary Material sup . From such relations, we arrive at the following connections between the band Chern numbers of $H_{\text{sqrt}}$ (denoted as $C_{n}^{\text{sqrt}}$) and those of $\mathbb{H}_{\pm}$ (denoted as $C_{|n|}^{\pm}$):

$$C_{n}^{\text{sqrt}}=\begin{cases}C_{|n|}^{+}=C_{|n|}^{-}&\text{if $|n|>1$}\\ C_{1}^{+}&\text{if $n=1$}\\ C_{1}^{-}&\text{if $n=-1$}\end{cases}.$$

(29)

These relations, schematically summarized in Fig. 8(b), are not limited to the forms of $M$ in Eq. (6) and independent of the value of $\Delta\neq 0$ (with a fixed sign). Therefore, the symmetric energy bands of $H_{\text{sqrt}}$ are accompanied by symmetric Chern numbers, i.e. $C_{n}^{\text{sqrt}}=C_{-n}^{\text{sqrt}}$, except for the first conduction ($n=1$) and valence ($n=-1$) bands that contain the threshold modes.

Incidentally, in the case of $\Delta=0$, one readily notices that chiral symmetry also underlies the identical Chern numbers between the $|n|>1$ conduction and valence bands (the Chern numbers for individual $n=\pm 1$ bands become ill-defined if there is degeneracy at the zero energy, they are equal if the zero modes are absent), but in the generic case this equality relation is traced back to Eq. (12)– a consequence of SUSY. And with the chiral symmetry at $\Delta\equiv 0$, topology of the square-root Hamiltonian is also noted in some lattice models Arkinstall et al. (2017); Ezawa (2020); Mizoguchi et al. (2020); Yoshida et al. (2021).

When chiral symmetry is broken in the case of $\Delta\neq 0$, our findings give insights on how contrasted topology of conduction and valence bands arises from the SUSY of the squared Hamiltonian, as dictated by the zero modes and the dispersion in their vicinity in the SUSY spectrum (e.g., Fig. 2). With $\Delta>0$, the Chern numbers are well defined in the first conduction and valence bands of $H_{\text{sqrt}}$, and $C_{1}^{\text{sqrt}}=C_{1}^{+}$ and $C_{-1}^{\text{sqrt}}=C_{1}^{-}$ are satisfied. So $C_{1}^{\text{sqrt}}\neq C_{-1}^{\text{sqrt}}$ when the dispersion in the vicinity of threshold/zero mode determines a finite difference between $C_{1}^{+}$ and $C_{1}^{-}$ [see e.g., Eq. (15) and Fig. 2]. The contrasted topology of the SPs of the squared Hamiltonian ($\mathcal{H}_{\text{susy}}$) in bands containing their zero modes are therefore transferred to the first conduction and valence bands of $H_{\text{sqrt}}$ respectively [Fig. 8(b)].

Figure S3 in the Supplementary Material sup confirms the above results with a simple example of $H_{\text{sqrt}}$ that describes massive Dirac fermions in a periodic magnetic field [see Fig. 1(a)]. Specifically, here $M=\pi_{-}$, and the associated SUSY system $\mathbb{H}_{\pm}$ are given by Eq. (9). By comparing Fig. S3(b) and Fig. 1(b), we find $C_{1}^{\text{sqrt}}=-1=C_{1}^{+}$ versus $C_{-1}^{\text{sqrt}}=0=C_{1}^{-}$, and $C_{n}^{\text{sqrt}}=C_{-n}^{\text{sqrt}}=C_{|n|}^{+}=C_{|n|}^{-}$ otherwise, which is consistent with Eq. (29) and Fig. 8(b). This model of $H_{\text{sqrt}}$, and its TR counterpart, can be realized by using strained graphene with a staggered potential Mahmud et al. (2023); Phong and Mele (2022); Wan et al. (2023b); Gao et al. (2023).

TBG in the chiral limit Tarnopolsky et al. (2019) and strained bilayer graphene (SBG) Wan et al. (2023b) are another two important examples where the contrasted topology of the conduction and valence minibands containing the threshold modes are dictated by Eq. (15) and Eq. (16). On one hand, they serve as cases with non-Abelian gauge fields San-Jose et al. (2012), on the other hand, they both exhibit $C_{-1}^{\text{sqrt}}-C_{1}^{\text{sqrt}}=2$, but due to drastically different origins. They both can be formulated in the form of $H_{\text{sqrt}}=\sigma_{x}\otimes(p_{x}+e\mathcal{A}_{x})+\sigma_{y}\otimes(p_{y}+e\mathcal{A}_{y})$, as detailed in Supplementary Material sup . Unlike a usual magnetic field, here $\boldsymbol{\mathcal{A}}=(\mathcal{A}_{x},\,\mathcal{A}_{y})$ is a non-Abelian gauge potential due to interlayer coupling, and the associated non-Abelian gauge field reads $\boldsymbol{\mathcal{B}}=\nabla\times\boldsymbol{\mathcal{A}}+i\frac{e}{\hbar}[\mathcal{A}_{x},\,\mathcal{A}_{y}]\hat{\boldsymbol{z}}$. They are matrices due to the layer degree of freedom in a hybridized bilayer. The associated SUSY Hamiltonian is given by the square of $H_{\text{sqrt}}$, i.e., $\mathbb{H}_{\pm}=\frac{1}{2m}(\boldsymbol{p}+e\boldsymbol{\mathcal{A}})^{2}\pm\boldsymbol{\mu}_{B}\cdot\boldsymbol{\mathcal{B}}$. The results discussed earlier, e.g., contrasted band topology in $\mathbb{H}_{\pm}$, thus can be applied directly. Let us focus on band topology and highlight the similarity/difference between TBG and SBG. In TBG, one expects $C_{1}^{-}-C_{1}^{+}=2$ because each SP has two zero modes (at the Dirac points, and we assume the twist angle is not magic Bistritzer and MacDonald (2011); Tarnopolsky et al. (2019)), around which the low-energy dispersion of $\mathbb{H}_{\pm}$ is quadratic [see Fig. 2(b)]. While $C_{1}^{-}-C_{1}^{+}=2$ is also expected in SBG, it instead results from one zero mode for each SP with quartic dispersion. If a sublattice staggered potential is added to introduce a finite $\Delta$ in $H_{\text{sqrt}}$, these results will manifest in the conduction and valence bands of $H_{\text{sqrt}}$ as $C_{-1}^{\text{sqrt}}-C_{+1}^{\text{sqrt}}=2$, which are consistent with numerical results (see Fig. S4(c) and Fig. S5(b) of Supplementary Material sup ).

The strained $\Gamma$-valley model with quadratic band touching Wan et al. (2023a) is an example that differs from all previous examples. Its Hamiltonian also has the form of $H_{\text{sqrt}}$, whereas $M=(p_{x}^{2}-p_{y}^{2}+f_{1})-i(2p_{x}p_{y}+f_{2})$ with $f_{1,2}(\boldsymbol{r})$ being some space-dependent functions, thus the effects of strain cannot be understood in terms of a gauge field. Nevertheless, topological flat bands and magic angles –like those in TBG– are also discovered. Our results [Fig. 8(b)] show that its band topology are inherited from the associated $\mathcal{H}_{\text{susy}}$ (despite in absence of an effective gauge field), thus provide insights for understanding the emergence of nontrivial topology and the distinct Chern numbers in the bands that contain the zero/threshold modes.

In addition to being fundamentally intriguing, the contrasted band topology in the lowest conduction and valence bands could be important for hosting nontrivial physical phenomena. Firstly, similar correlation-driven states as those discussed in earlier sections can be expected with electron or hole doping when interaction is strong. Furthermore, the interplay of distinct topology in the conduction and valence bands could lead to other interesting states of matter. For example, it was predicted recently that topological excitons with a finite vorticity in momentum space can emerge in systems with topologically distinct conduction and valence bands, which can exhibit optical selection rules of geometrical/topological origin as well as lead to interesting nonlinear anomalous Hall effects Xie et al. (2024). The square-root systems discussed here serve as natural playgrounds for exploring topological exciton insulators and condensates.