A monolayer transition metal dichalcogenide as a topological excitonic insulator

TL;DR

This paper predicts that monolayer T' phase transition metal dichalcogenides can form a topological excitonic insulator state, where strong electron-hole interactions induce a new ground state with enhanced topological and excitonic properties.

Contribution

It demonstrates, through first-principles calculations, the emergence of a topological excitonic insulator phase in monolayer T' transition metal dichalcogenides, revealing a novel interplay of topology and excitonic order.

Findings

Electron-hole binding energy exceeds the band gap in MoS2.

The excitonic insulator state enhances the bulk gap and breaks inversion symmetry.

Strain induces additional ferroelectric excitonic phases with distinct topological properties.

Abstract

Monolayer transition metal dichalcogenides in the T' phase promise to realize the quantum spin Hall (QSH) effect at room temperature, because they exhibit a prominent spin-orbit gap between inverted bands in the bulk. Here we show that the binding energy of electron-hole pairs excited through this gap is larger than the gap itself in MoS2, a paradigmatic material that we investigate from first principles by many-body perturbation theory (MBPT). This paradoxical result hints at the instability of the T' phase against the spontaneous generation of excitons, and indeed we find that it gives rise to a recostructed `excitonic insulator' ground state. Importantly, we show that in this system topological and excitonic order cooperatively enhance the bulk gap by breaking the crystal inversion symmetry, in contrast to the case of bilayers where the frustration between the two orders is relieved by breaking time reversal symmetry. The excitonic topological insulator departs distinctively from the bare topological phase as it lifts the band spin degeneracy, which results in circular dichroism. A moderate biaxial strain applied to the system leads to two additional excitonic phases, different in their topological character but both ferroelectric as an effect of electron-electron interactions.

Equations

$${\hat{\mathcal{H}}}_{\text{QSH}}=\frac{{\epsilon }_p+{\epsilon }_d}{2}{\hat{\mathbb{I}}}_{\tau }\otimes {\hat{\mathbb{I}}}_{\sigma }+\frac{{\epsilon }_p-{\epsilon }_d}{2}{\hat{\tau }}_z\otimes {\hat{\mathbb{I}}}_{\sigma }-\text{Im}\{{\Delta }_{\text{SO}}\}{\hat{\tau }}_y\otimes {\hat{\mathbb{I}}}_{\sigma }+\text{Re}\{{\Delta }_{\text{SO}}\}{\hat{\tau }}_x\otimes {\hat{\sigma }}_x.$$

$${\hat{\mathcal{H}}}_{\text{X}}=\text{Re}\{{\Delta }_{\text{X}}\}{\hat{\tau }}_x\otimes {\hat{\mathbb{I}}}_{\sigma }-\text{Im}\{{\Delta }_{\text{X}}\}{\hat{\tau }}_y\otimes {\hat{\sigma }}_x,$$

$$\eta (\text{k})=\frac{\mathcal{A}({\sigma }^{+})-\mathcal{A}({\sigma }^{-})}{\mathcal{A}({\sigma }^{+})+\mathcal{A}({\sigma }^{-})},$$

$$\text{Re}\{{\Delta }_{\text{X}}(\text{k})\}=-\sum _{{\text{k}}^{\prime }}W(\text{k}-{\text{k}}^{\prime })⟨p_{{\text{k}}^{\prime }\uparrow }^{†}{d}_{{\text{k}}^{\prime }\uparrow }⟩,$$

$$\text{Im}\{{\Delta }_{\text{X}}(\text{k})\}=-\sum _{{\text{k}}^{\prime }}W(\text{k}-{\text{k}}^{\prime })⟨p_{{\text{k}}^{\prime }\uparrow }^{†}{d}_{{\text{k}}^{\prime }\downarrow }⟩,$$

$$\text{Re}\{{\Delta }_{\text{X}}(\text{k})\}$$

$$\text{Im}\{{\Delta }_{\text{X}}(\text{k})\}$$

$$P(\text{k})=⟨\text{k},v,\lambda |\hat{\Theta }{|\text{k},v,-\lambda ⟩}_{\text{sp}}=u_0{u}_0^{\prime }-v_0{v}_0^{\prime }{e}^{i(\phi -{\phi }^{\prime })},$$

$$\mathcal{A}({\sigma }^{\pm })=\frac{2\pi }{\hslash }\sum _{\lambda }{\left|⟨{\Psi }_0|\mathcal{D}({\sigma }^{\pm }){\gamma }_{\text{k}c\lambda }^{†}{\gamma }_{\text{k}v\lambda }|{\Psi }_0⟩\right|}^2\delta \left(\hslash \omega -E_c(\text{k})+E_v(\text{k})\right),$$

$$\mathcal{D}({\sigma }^{\pm })=e\left({\mathcal{E}}_x{x}_{\Gamma }\pm i{\mathcal{E}}_z{z}_{\Gamma }\right)\sum _{\text{k}\sigma }\left(p_{\text{k}\sigma }^{†}{d}_{\text{k}\sigma }+d_{\text{k}\sigma }^{†}{p}_{\text{k}\sigma }\right),$$

$$⟨\text{P}⟩=\frac{e\left(x_{\Gamma }\text{i}+z_{\Gamma }\text{k}\right)}{N_x{a}_x{N}_y{a}_y}\sum _{\text{k}\sigma }⟨p_{\text{k}\sigma }^{†}{d}_{\text{k}\sigma }+d_{\text{k}\sigma }^{†}{p}_{\text{k}\sigma }⟩,$$

Full text

{affiliations}

CNR-NANO, Via Campi 213a, 41125 Modena, Italy.

INFN, Dipartimento di Fisica, Università degli Studi di Roma Tor Vergata, Via Della Ricerca Scientifica 1, 00133 Roma, Italy.

Dipartimento di Scienze Fisiche, Informatiche e Matematiche, Università degli Studi di Modena e Reggio Emilia, 41125 Modena, Italy.

The monolayer transition metal dichalcogenides that were recently proposed as candidates for the QSH effect all have overlapping metal-$d$ conduction and chalcogenide-$p$ valence bands[2]. Such ‘band inversion’ makes the system either a narrow-gap semiconductor, due to $p-d$ spin-orbit hybridization (Fig. 1c), or a semimetal, whose band edges are displaced in momentum space. In both cases long-range Coulomb attraction, which is poorly screened in two dimensions, tends to bind electrons ($e$) at the bottom of conduction band with holes ($h$) at the top of valence band, thus giving rise to excitons. If the e-h binding energy is larger than the semiconductor gap (or if it is non vanishing in the semimetal), then excitons will spontaneoulsy form and condense, until a correlated gapped phase is built at thermodynamic equilibrium: the excitonic insulator proposed in the sixties[5, 6, 7, 8]. This paradigm has been recently invoked for the QSH insulator WTe2, as its bulk gap is strongly sensitive to temperature[21, 22, 23] and doping[24], whereas in the absence of excitonic effects MBPT predicts semimetallic behaviour[2]. Mounting evidence of the excitonic insulator has been accumulating in the last two years in low-dimensional materials[25, 26, 27, 16]—noticeably transition metal dichalcogenides[28, 29, 30].

The relation between topological and excitonic order is intriguing, as the former emerges in a noninteracting picture whereas the latter is driven by e-h interactions. So far the problem has been discussed for bilayers[14, 15, 16, 17, 18]: the spin-orbit gap, associated with topological order, depends on interlayer tunneling but the exciton binding does not, as $e$ and $h$ in separate layers remain bound by long-range attraction. This leads to scenarioes of frustration between the QSH phase and the topologically trivial excitonic phase. The monolayer case is different, since screening is suppressed as the gap opens: here both spin-orbit gap and exciton binding are affected by interband hybridization.

Here we study an archetypical member of the $T^{\prime}$ family of QSH candidates, monolayer MoS2, through first-principles calculations by means of MBPT, and take into account e-h binding by solving the Bethe-Salpeter equation[31]. This allows us to assess the inherent excitonic instability of the material. We then demonstrate the coexistence of topological and excitonic orders through a self-consistent approach, which predicts a chiral ground state wave function with unique fingerprints.

Results

Exciton binding and instability. In spite of the small point symmetry group, $T^{\prime}$-MoS2 has an inversion center, located at the midpoint between two neighbour Mo atoms (black dot in Fig. 1a): these form a zig-zag chain, parallel to the $y$ axis, which is characteristic of the $T^{\prime}$ phase. We compute the energy bands from first principles, including spin-orbit interaction at the density functional theory (DFT) level, and then evaluating many-body corrections within the GW approximation (Methods). The resulting band structure is highly anisotropic (Fig. 1b), the conduction band being almost flat in the direction perpendicular to the zig-zag chain. The cut of energy surfaces along the $\Gamma\Lambda$ direction (dashed red line in the Brillouin zone domain of Fig. 1b) clarifies why $T^{\prime}$-MoS2 is a QSH insulator. As shown in Fig. 1c, the overlap of $p$ and $d$ bands forms two valleys located at k = $\pm\Lambda$, similar to Dirac valleys of graphene. Like in the Kane-Mele model[1], spin-orbit interaction opens a gap, acting as a valley-dependent magnetic field for each spin species: hence the valley contributions to the topological invariant $\mathbb{Z}_{2}$ do not cancel out[1, 2]. Importantly, the GW renormalization of the DFT gap is gigantic (respectively solid and dashed lines in Fig. 1c), which points to the relevance of electronic interactions.

From the solution of Bethe-Salpeter equation (Methods) we find that the binding energy of the lowest exciton exceeds the GW gap by 32 meV, hence $T^{\prime}$-MoS2 is unstable against the spontaneous generation of excitons. The exciton probability weight in momentum space is localized in the two $\Lambda$ valleys (Fig. 2b): the weight is stretched along the $k_{x}$ direction, following the energy profile of uncorrelated e-h pairs excited with zero total momentum. This exciton is eight-fold degenerate within numerical accuracy, as multiplet states include both bonding and antibonding combinations of the wave functions in the two valleys, for all possible $e$ and $h$ spin projections along the two-fold screw axis $y$ (cf. Fig. 1a); the spin full rotational symmetry is reduced to a rank-two representation by spin-orbit interaction. The exciton wave function in real space is shown in Fig. 2a as the conditional probability of finding a bound electron (red colour), provided the hole position is fixed (black dot): this density is substantially delocalized in the direction of the Mo zig-zag chain, consistently with Fig. 2b.

Topological excitonic insulator. The excitonic instability reminds us of the Cooper problem of two bound electrons filling a Fermi sea: whereas the gluing of $e-e$ pairs heralds the transition to the Bardeen-Cooper-Schrieffer superconductor, the strong binding of $e-h$ pairs signals the formation of the excitonic insulator: both Cooper pairs and excitons collectively enforce the many-body gap, which is the self-consistent order parameter[32]. Were $T^{\prime}$-MoS2 a semimetal (Fig. 3a), the analogy would be complete at the formal level, as the excitonic order parameter, $\Delta_{\text{X}}$, would then hybridize conduction and valence bands, opening a gap[7]. However, the actual bands are effectively hybridized by spin-orbit interaction, $\Delta_{\text{SO}}$ (Fig. 3b), hence the role of $\Delta_{\text{X}}$ requires further clarification (Fig. 3c).

Following the seminal work by Volkov and Kopaev[9], a key observation is that $\Delta_{\text{SO}}$ and $\Delta_{\text{X}}$ have opposite parity in k space. In order to preserve the inversion symmetry of the crystal, $\Delta_{\text{SO}}$ must be odd, $\Delta_{\text{SO}}(\text{\bf k})=-\Delta_{\text{SO}}(-\text{\bf k})$, as the periodic parts of conduction and valence states transform like $p_{y}$ and $d_{yz}$ orbitals, respectively, as illustrated by their wave functions at $\Gamma$ (insets in Fig. 1c). Since $\Delta_{\text{X}}$ is associated to the lowest-exciton wave function in k space[32], it must have $s$-wave symmetry, $\Delta_{\text{X}}(\text{\bf k})=\Delta_{\text{X}}(-\text{\bf k})$. Besides, the spin degeneracy associated with the exciton that drives the instability rules out the breaking of time reversal symmetry. Together, these conditions provide us with a tractable two-band model, by reducing the number of independent order parameters from eight (the spin-resolved, complex interband hybridizations) to two, i.e., the real and imaginary parts of $\Delta_{\text{X}}(\text{\bf k})$ (see Methods).

The mean-field Hamiltonian in k space, $\hat{\cal{H}}(\text{\bf k})=\hat{\cal{H}}_{\text{QSH}}(\text{\bf k})+\hat{\cal{H}}_{\text{X}}(\text{\bf k})$, is therefore a 4 $\times$ 4 matrix that acts on vectors spanned by spin-resolved $p$ and $d$ Bloch states at $\Gamma$. It adds the excitonic term, $\hat{\cal{H}}_{\text{X}}$, to the single particle term, $\hat{\cal{H}}_{\text{QSH}}$, based on an effective-mass model[2] of GW bands:

[TABLE]

Here the model parameters are the band inversion and effective masses of $p$ and $d$ energy bands, $\varepsilon_{p}(\text{\bf k})$ and $\varepsilon_{d}(\text{\bf k})$ [plotted in Fig. 3d as $\xi(\text{\bf k})=\left(\varepsilon_{p}-\varepsilon_{d}\right)/2$ ], as well as the velocities, $v_{1}$ and $v_{2}$, associated with the complex spin-orbit interaction, $\Delta_{\text{SO}}(\text{\bf k})=\hbar v_{2}k_{y}-i\hbar v_{1}k_{x}$ (Figs. 3e and 3f). $\hat{\sigma}_{\alpha}$ and $\hat{\mathbb{I}}_{\sigma}$ are the 2 $\times$ 2 Pauli matrices and identity in spin space, whereas $\hat{\tau}_{\alpha}$ and $\hat{\mathbb{I}}_{\tau}$ act on the pseudospin space of $p$ and $d$ orbital components ($\alpha=x,y,z$). The chosen parameters provide good matching between model and first-principles GW bands (respectively dashed and dotted lines in Figs. 4a and 4b) in the k-space region of interest. The excitonic hybridization, $\Delta_{\text{X}}(\text{\bf k})$, which appears in

[TABLE]

is obtained numerically (Figs. 3g and 3h) by solving two self-consistent coupled equations [see (6) and (7) below], ruled by the screened Coulomb interaction, $W(\text{\bf k})$, which we extract from first principles (Supplementary Figure 1). In the QSH phase these equations only admit the trivial solution, $\Delta_{\text{X}}=0$. The agreement between effective-mass (Supplementary Figure 2) and first-principles (Fig. 2b) exciton wave functions points to the reliability of the two-band model.

The conduction and valence bands of the excitonic insulator are respectively $E_{c\lambda}(\text{\bf k})=(\varepsilon_{p}+\varepsilon_{d})/2+E_{\text{\bf k}\lambda}$ and $E_{v\lambda}(\text{\bf k})=(\varepsilon_{p}+\varepsilon_{d})/2-E_{\text{\bf k}\lambda}$, where $E_{\text{\bf k}\lambda}=\left[\xi^{2}+\left|\Delta_{\text{SO}}\;\lambda\;\Delta_{\text{X}}\right|^{2}\right]^{1/2}$ and $\lambda=\pm$ is a quantum label that reduces to the spin projection in the absence of spin-orbit interaction. Low-lying $e-h$ excitations close to the gap (Fig. 4c), of energies $E_{c\lambda}(\text{\bf k})-E_{v\lambda^{\prime}}(\text{\bf k})$, occur at the pristine semimetal Fermi surface, which is an ellipse in k space obeying $\xi=0$ (white colour in Fig. 3d). In particular, as shown in Figs. 3e-h, at $\Lambda$ points both $\text{Im}\{\Delta_{\text{SO}}\}$ and $\text{Re}\{\Delta_{\text{X}}\}$ vanish, hence $\text{Re}\{\Delta_{\text{SO}}\}$ and $\text{Im}\{\Delta_{\text{X}}\}$ add quadratically without interfering, resulting in the approximate gap value $E_{c\lambda}(\Lambda)-E_{v\lambda^{\prime}}(\Lambda)=2\left[(\text{Re}\{\Delta_{\text{SO}}\})^{2}+(\text{Im}\{\Delta_{\text{X}}\})^{2}\right]^{1/2}$. Therefore, the insulator is simultaneously excitonic and topological (QSHX), as an effect of the self-organization of $\Delta_{\text{X}}$(k). This may be checked most easily by adiabatic continuation, since the pristine QSH gap smoothly increases as the temperature is lowered below the critical value of the QSHX phase, around 700 K (blue curve in Fig. 5b). We explicitly calculate $\mathbb{Z}_{2}$ in Supplementary Note 1.

Fingerprints The excitonic hybridization spontaneously breaks the inversion symmetry of the crystal ($\hat{\cal{H}}_{\text{X}}$ anticommutes with the inversion operator, $\hat{I}=-\hat{\tau}_{z}\otimes\hat{\mathbb{I}}_{\sigma}$), which leads to clear-cut observable properties. The Kramers degeneracy of bands is lifted, as shown in Fig. 4b (black and red solid lines point to the two spin projections), with band splittings reaching a maximum of about 20 meV for $k_{x}\sim 0.05-0.1$ [units of $(2\pi)/a_{x}$] and vanishing on the $\Gamma\Lambda$ line (solid black line in Fig. 4a), as a consequence of the two-fold rotational symmetry that restores degeneracy. The QSHX phase exhibits circular dichroism, like monolayer $H$-MoS2 (Refs. 33, 34): the optical absorption of a photon whose wave vector lies in the $yz$ plane (Fig. 4d), $\cal{A}$, depends on the circular polarization of the photon itself, $\sigma^{+}$ or $\sigma^{-}$. Figure 4e shows the k-resolved degree of optical polarization[33],

[TABLE]

evaluated for the optical transition from the highest valence to the lowest conduction band, whose excitation energy, $E_{c\lambda}(\text{\bf k})-E_{v\lambda}(\text{\bf k})$, is shown in Fig. 4c. Remarkably, as $\eta$(k) circles around the semimetal Fermi surface, it modulates from its lower to its upper bound, respectively $-1$ (blue color in Fig. 4e) and $+1$ (red colour). These limit values are effectively optical selection rules coupling orbital and spin degrees of freedom, like in the case of monolayer $H$-MoS2. On the contrary, in the QSH ground state $\eta$(k) $=0$ at each k point, since the microscopic transverse currents, which are responsible of the net angular momentum, exactly cancel out due to Kramers degeneracy[35]. We expect the dichroic signal to survive to final-state interactions, not included in our calculation (Methods), and hence disclose the intrinsic chirality of the QSHX ground state.

Excitonic phases and ferroelectricity The application of biaxial strain to $T^{\prime}$-MoS2 is a practical handle to tune the band inversion[2] and hence the energy scales ruling the ground state, as shown in Fig. 5a. Excitonic correlations tend to destroy topological order as the band overlap is suppressed, since $\Delta_{\text{SO}}$ decreases while $\Delta_{\text{X}}$ increases. Such balance eventually leads to a topologically trivial excitonic insulator (region in yellow colour, X-FE), after crossing a phase allowing for destructive interference between $\Delta_{\text{SO}}$ and $\Delta_{\text{X}}$ (cyan colour, QSHX-FE). This tiny intermediate region, located around the value of 180 meV, is broadened by temperature up to the frontier with the QSH phase. The topological character of each phase is made evident by the cuts of the phase diagram along the temperature axis (violet and blue dashed lines in Fig. 5a) displayed in Fig. 5b. Starting from the QSH phase and lowering the temperature, the gap increases when entering the QSHX phase (blue curve) but decreases while crossing the QSHX-FE region (violet curve), until the gap closes and opens again in the, now topologically trivial, X-FE phase (cf. Supplementary Note 1).

Both QSHX-FE and X-FE phases totally distort the pristine $C_{2h}$ symmetry by breaking the screw axis symmetry along $y$ (Supplementary Figure 4), in addition to inversion. This is related to the macroscopic condensation of the exciton electric dipole in the $xz$ plane, which makes the system ferroelectric[19]. Contrary to usual ferroelectrics, like BaTiO3, here the permanent electric dipole, P, is not due to the displacement of anions and cations but to the modulation of the electronic charge density associated with the exciton polarization (Fig. 5c). This may open fascinating new routes, like the realization of ultrafast switches between conductive and insulating (ferroelectric) beahvior, locally controllable by strain or screening, or the exploration of exotic electronic collective modes[19], coherently radiating in the THz range.

In conclusion, we have demonstrated that a paradigmatic two-dimensional topological insulator is also an excitonic insulator, by combining calculations from first principles with a self-consistent mean-field model. We expect our results to be relevant to other $T^{\prime}$ polytypes, including those that are semimetal according to MBPT, like WTe2, which are harder to simulate in view of the enhanced screening of $e-h$ interactions. In particular, WTe2 owns the record temperature of 100 K for the QSH effect[3] and exhibits unexpected but intriguing properties that might be related to the excitonic insulator, like ferroelectricity[20] and gate-induced superconductivity[36, 24] in close proximity to the insulating phase.

{methods}

0.1 Computational details of many-body

perturbation theory from first principles.

The ground-state structure and Kohn-Sham states were calculated for a single layer of $T^{\prime}$-MoS2 by using a DFT approach, as implemented in the Quantum ESPRESSO package[37]. The generalized gradient approximation (GGA) PBE parametrization[38] was adopted together with plane wave basis set and norm-conserving pseudopotentials to model the electron-ion interaction. Fully relativistic pseudopotentials treating the $sp$ semicore states of the transition metal atoms as valence electrons were employed. The kinetic energy cutoff for the wave functions was set to 90 Ry. The Brillouin zone was sampled by using a 16 $\times$ 16 $\times$ 1 k-point grid. The supercell size perpendicular to the $T^{\prime}$-MoS2 layer was set to $a_{z}=15.98$ Å and checked to be large enough to avoid spurious interactions with its replica.

The equilibrium atomic lattice parameters were obtained by performing a full relaxation of the cell and atomic positions. The obtained equilibrium lattice parameters, $a_{x}=5.74$ Å, $a_{y}=3.19$ Å, as well as the Kohn-Sham electronic gap were in very good agreement with previous literature[2].

Many-body perturbation theory[4] calculations were performed using the Yambo code[39, 40]. Many-body corrections to the Kohn-Sham eigenvalues were calculated within the $G0W0$ approximation to the self-energy operator, where the dynamic dielectric function was obtained within the plasmon-pole approximation[41]. The spectrum of excited states was then computed by solving the Bethe-Salpeter equation (BSE). The static screening in the direct term was calculated within the random-phase approximation with inclusion of local field effects; the Tamm-Dancoff approximation for the BSE Hamiltonian was employed after having verified that the correction introduced by coupling the resonant and antiresonant part was negligible. Converged excitation energies were obtained considering respectively 2 valence and 2 conduction bands in the BSE matrix. For the calculations of the GW band structure and the Bethe-Salpeter matrix, the Brillouin zone was sampled with a 70 $\times$ 35 $\times$ 1 k-point grid. A kinetic energy cutoff of 60 Ry was used for the evaluation of the exchange part of the self energy and 10 Ry for the screening matrix size. 248 unoccupied bands were used to build the polarizability and integrate the self-energy. The Coulomb interaction was truncated[42] in the layer-normal direction to avoid spurious interactions with the image systems. Note that the GW gap obtained in this work (0.26 eV) is similar to the value reported in Ref. 43 calculated at HSE06 level (0.23 eV). A smaller value of the $G0W0$ gap (0.08 eV) is instead reported in Ref. 2. The observed discrepancy can be probably explained by the use of a non-truncated Coulomb potential in that work. The approach used in the GW and BSE calculations takes into account the full spinorial nature of the electronic Kohn-Sham states, providing superior accuracy than perturbative treatments of spin-orbit coupling.

Two-band model and self-consistent mean-field theory. The effective-mass Hamiltonian, $\hat{\cal{H}}(\text{\bf k})=\hat{\cal{H}}_{\text{QSH}}(\text{\bf k})+\hat{\cal{H}}_{\text{X}}(\text{\bf k})$, acts on four-component vectors, $\left(u_{\uparrow}(\text{\bf k}),v_{\uparrow}(\text{\bf k}),u_{\downarrow}(\text{\bf k}),v_{\downarrow}(\text{\bf k})\right)$, with $u_{\sigma}$(k) and $v_{\sigma}$(k) being the spin- and k-resolved envelope functions[44] of the $p_{y}$ and $d_{yz}$ Bloch states at $\Gamma$, respectively, and $\sigma=\uparrow,\downarrow$. The vector normalization is such that $\sum_{\sigma}(\left|u_{\sigma}(\text{\bf k})\right|^{2}+\left|v_{\sigma}(\text{\bf k})\right|^{2})=1$.

QSH Hamiltonian. The non-interacting term of Eq. (1), $\hat{\cal{H}}_{\text{QSH}}$, is modelled after Ref. 2 to provide the bands of the QSH insulator—renormalized within the GW approximation—by accounting for the spin-orbit interaction. It complies with the $C_{2h}$ point symmetry group, which includes the inversion, $\hat{I}=-\hat{\tau}_{z}\otimes\hat{\mathbb{I}}_{\sigma}$, and the two-fold screw axis rotation along $y$, $\hat{C}_{2y}=i\,\hat{\tau}_{z}\otimes\hat{\sigma}_{y}$. The band dispersions are $\varepsilon_{p}(\text{\bf k})=-\delta\;-\hbar^{2}k_{x}^{2}/(2m_{px})\;-\hbar^{2}k_{y}^{2}/(2m_{py})$ and $\varepsilon_{d}(\text{\bf k})=\delta\;+\hbar^{2}k_{x}^{2}/(2m_{dx})\;+\hbar^{2}k_{y}^{2}/(2m_{dy})$, with $-\delta$ being the band inversion. We optimize the matching with first-principles GW bands by choosing $\delta=-250$ meV, $m_{px}/m_{e}=0.6$, $m_{py}/m_{e}=0.31$, $m_{dx}/m_{e}=2.48$, $m_{dy}/m_{e}=0.5$, $v_{1}=$ 1.5$\cdot 10^{15}$ Å/s, $v_{2}=$ 1.3$\cdot 10^{15}$ Å/s ($m_{e}$ is the free electron mass).

Excitonic hybridization. The mean-field excitonic Hamiltonian, $\hat{\cal{H}}_{\text{X}}$(k), in principle accounts for eight independent order parameters, i.e., the four complex hybridization terms coupling all possible pairs of conduction- and valence-band spin projections. The form (2) complies with the constrain that $\hat{\cal{H}}_{\text{X}}$ is even in k, $\hat{\cal{H}}_{\text{X}}$(k) $=\hat{\cal{H}}_{\text{X}}(-\text{\bf k})$, and invariant under time reversal, $\hat{\Theta}\hat{\cal{H}}_{\text{X}}(-\text{\bf k})=\hat{\cal{H}}_{\text{X}}(\text{\bf k})\hat{\Theta}$, where $\hat{\Theta}=i\,\hat{\mathbb{I}}_{\tau}\otimes\hat{\sigma}_{y}\hat{K}$ and $\hat{K}$ are the time-reversal and the complex conjugation operator, respectively. The excitonic hybridization terms are derived from the interband Coulomb interaction through the usual mean-field decoupling procedure, as

[TABLE]

and

[TABLE]

where the Fermi field operators $p_{\text{\bf k}\sigma}$ and $d_{\text{\bf k}\sigma}$ destroy electrons occupying $p$ and $d$ Bloch states, respectively, and the symbol $\left<\ldots\right>$ stands for the quantum statistical average [the expressions (4) and (5) do not change as one reverses the spins of $p^{\dagger}_{\text{\bf k}\sigma}$ and $d_{\text{\bf k}\sigma^{\prime}}$ simultaneously]. Whereas quantum and thermal fluctuations are believed to disrupt long range order in low dimensions, here we expect that excitonic correlations are stabilized by the long range of the Coulomb interaction, as verified in the one-dimensional case through extensive quantum Monte Carlo simulations[25].

Dressed interaction. The key quantity above is the q-resolved screened Coulomb interaction[45], $W(\text{\bf q})=V_{0}(\text{\bf q})/(1+2\pi\alpha_{\text{2D}}\left|\text{\bf q}\right|)$, which we fit to the one obtained from first principles within the random phase approximation (Supplementary Figure 1). The long-wavelength term, $V_{0}(\text{\bf q})\sim 1/\left|\text{\bf q}\right|$, is unscreened in two dimensions[45] and hence the dominant contribution to $e-h$ attraction. Since the prefactor of $1/\left|\text{\bf q}\right|$ is determined by the dimensions of the supercell used in the first-principles calculation[42], i.e., $a_{x}=5.743$ Å, $a_{y}=3.191$ Å, $a_{z}=15.98$ Å, the only free parameter is the two-dimensional polarizability[45], $\alpha_{\text{2D}}$. We adjust $\alpha_{\text{2D}}$ to match effective-mass and first-principles exciton binding energies, with $\alpha_{\text{2D}}=10.75$ corresponding to the three-dimensional dielectric constant of 9.5 for the stack of $T^{\prime}$-MoS2 monolayers—the actual bulk of the supercell calculation[45]. This latter figure reasonably compares with the values—between 10 and 20—assumed by the first-principles dielectric function in the $q$-range of $0.02\div 0.05$ [units of $(2\pi)/a_{x}$].

Eigenvectors. The eigenvectors of $\hat{\cal{H}}(\text{\bf k})$ are the conduction and valence band states of the correlated insulator. The first of the two valence bands, $(v,\lambda=+)$, has as wave function envelopes $u_{\sigma}(\text{\bf k},v,+)=u_{0}/\sqrt{2}$ and $v_{\sigma}(\text{\bf k},v,+)=iv_{0}\exp(-i\varphi)/\sqrt{2}$, with $u_{0}$ and $v_{0}$ positive numbers whose magnitudes are given by $u_{0}^{2}=1/2[1-\xi(\text{\bf k})/E_{\text{\bf k}+}]$ and $u_{0}^{2}+v_{0}^{2}=1$; the phase is provided by $iu_{0}v_{0}\exp(i\varphi)=(\Delta_{\text{SO}}+\Delta_{\text{X}})/2E_{\text{\bf k}+}$. The second band, $(v,\lambda=-)$, has $u_{\sigma}(\text{\bf k},v,-)=-\sigma u_{0}^{\prime}/\sqrt{2}$ and $v_{\sigma}(\text{\bf k},v,-)=-\sigma iv_{0}^{\prime}\exp(i\varphi^{\prime})/\sqrt{2}$ (the spin index takes the values $\sigma=\pm$ when occurring in the body of a formula), with $u_{0}^{\prime 2}=1/2[1-\xi(\text{\bf k})/E_{\text{\bf k}-}]$, $u_{0}^{\prime 2}+v_{0}^{\prime 2}=1$, and $iu_{0}^{\prime}v_{0}^{\prime}\exp(i\varphi^{\prime})=(\Delta_{\text{SO}}-\Delta_{\text{X}})/2E_{\text{\bf k}-}$. The envelopes of conduction bands $(c,\lambda)$ have similar expressions, which are obtained from $u_{\sigma}(\text{\bf k},v,\lambda)$ and $v_{\sigma}(\text{\bf k},v,\lambda)$ by swapping $u_{0}$ for $v_{0}$ ($u_{0}^{\prime}$ for $v_{0}^{\prime}$) in the formulae and simultaneously adding $\pi$ to $\varphi$ ($\varphi^{\prime}$).

Self-consistent equations. The many-body ground state, $\left|\Psi_{0}\right>$, is the Slater determinant with all valence band states filled, i.e., $\left|\Psi_{0}\right>=\prod_{\text{\bf k}}\gamma^{\dagger}_{\text{\bf k},v,+}\gamma^{\dagger}_{\text{\bf k},v,-}\left|0\right>$, where $\left|0\right>$ is the vacuum and $\gamma_{\text{\bf k},\alpha,\lambda}=\sum_{\sigma}u_{\sigma}(\text{\bf k},\alpha,\lambda)\,p_{\text{\bf k}\sigma}+v_{\sigma}(\text{\bf k},\alpha,\lambda)\,d_{\text{\bf k}\sigma}$ is the Fermi operator destroying an electron occupying the eigenstate of $\hat{\cal{H}}(\text{\bf k})$ of wave vector k and band index $(\alpha,\lambda)$, with $\alpha=c,v$. The knowledge of many-body states, together with Eqs. (4) and (5), allows us to write explicitly the self-consistent equations for $\Delta_{\text{X}}$,

[TABLE]

and

[TABLE]

Here $f_{F}(E)=(\exp{\beta(E-\mu)}+1)^{-1}$ is the Fermi distribution function, $\beta=1/k_{B}T$ is inversely proportional to the temperature $T$, $k_{B}$ is Boltzmann’s constant, and $\mu$ is the equilibrium chemical potential, which we place at midgap neglecting small deviations expected at finite temperature. Note that the divergence of $W$ at long wavelength is harmless since the interaction is integrated over the Brillouin zone. Equations (6) and (7), which are coupled together essentially through the denominators $E_{\text{\bf k}^{\prime}\lambda}$, allow for the trivial solution $\Delta_{\text{X}}(\text{\bf k})\equiv 0$ in the QSH phase. The non trivial solution is obtained numerically through recursion, using the k-resolved exciton wave function as a first-iteration seed[25], which allows for quick and robust convergence.

Topological invariant $\mathbb{Z}_{2}$. We derive the topological invariant $\mathbb{Z}_{2}$ by applying the test originally developed by Kane and Mele for graphene[46] to $\hat{\cal{H}}(\text{\bf k})$. This relies on the evaluation of the overlap, $P(\text{\bf k})$, between the single-particle state $(\text{\bf k},v,\lambda)$ and the time reversed of $(\text{\bf k},v,-\lambda)$,

[TABLE]

with $\left<...\right>_{\text{sp}}$ being the scalar product between four-component vectors. In Supplementary Note 1 we assess the number of pairs of zeroes of $P(\text{\bf k})$ occurring at ${\text{\bf k}}^{*}$ and $-{\text{\bf k}}^{*}$, which at once characterizes the topology of the many-body state.

Circular dichroism and permanent electric dipole. The degree of optical polarization, $\eta$(k), is defined by Eq. (3) through the optical absorption, ${\cal A}(\sigma^{\pm})$, of a photon having circular polarization, $\sigma^{+}$ or $\sigma^{-}$, and wave vector lying in the $yz$ plane. We evaluate ${\cal A}(\sigma^{\pm})$ through Fermi golden rule,

[TABLE]

by selecting the mimimum $e-h$ pair excitation energy available at a given point in k space, $E_{c}(\text{\bf k})-E_{v}(\text{\bf k})$, i.e., the photon energy $\hbar\omega$ depends on k, as shown in Fig. 4c. The light-matter interaction within the dipole approximation,

[TABLE]

depends on the (real) coordinate interband matrix elements evaluated at $\Gamma$, $x_{\Gamma}$ and $z_{\Gamma}$, which we extract from first principles ($x_{\Gamma}=1.94$ Å, $z_{\Gamma}=0.308$ Å, $y_{\Gamma}=0$ due to symmetry). Here ${\cal{E}}_{x}={\cal{E}}_{z}={\cal{E}}$ is the electric field, and the sign $\pm$ picks up the helicity of the field that is circularly polarized in the $xz$ plane. The permanent electric dipole, $\left<\text{\bf P}\right>$, is the equilibrium average of the periodic part of the interband dipole operator over the many-body ground state[19],

[TABLE]

where $N_{x}N_{y}$ is the number of unit cells used in the simulation.