Emergent dual topology in the three-dimensional Kane-Mele Pt$_2$HgSe$_3$

Recently, the very first large-gap Kane-Mele quantum spin Hall insulator was predicted to be monolayer jacutingaite (Pt$_2$HgSe$_3$), a naturally-occurring exfoliable mineral discovered in Brazil in 2008. The stacking of quantum spin Hall monolayers into a van-der-Waals layered crystal typically leads to a (0;001) weak topological phase, which does not protect the existence of surface states on the (001) surface. Unexpectedly, recent angle-resolved photoemission spectroscopy experiments revealed the presence of surface states dispersing over large areas of the 001-surface Brillouin zone of jacutingaite single crystals. The 001-surface states have been shown to be topologically protected by a mirror Chern number $C_M=-2$, associated with a nodal line gapped by spin-orbit interactions. Here, we extend the two-dimensional Kane-Mele model to bulk jacutingaite and unveil the microscopic origin of the gapped nodal line and the emerging crystalline topological order. By using maximally-localized Wannier functions, we identify a large non-trivial second nearest-layer hopping term that breaks the standard paradigm of weak topological insulators. Complemented by this term, the predictions of the Kane-Mele model are in remarkable agreement with recent experiments and first-principles simulations, providing an appealing conceptual framework also relevant for other layered materials made of stacked honeycomb lattices.

Graphene's crystal and electronic structure has been fundamental for the development of the theory of topological insulators. The very first model of a topological insulator ever proposed, namely the Chern (or quantum anomalous Hall) insulator by Haldane, is essentially a two-band tight-binding model for graphene in the presence of a staggered magnetic field with zero flux over the unit cell [1]. The experimental isolation of graphene [2] inspired Kane and Mele to assert that by doubling Haldane's model and introducing spins one could describe intrinsic spin-orbit coupling (SOC) in graphene, leading to a novel gapped topological phase [3,4]. Such phase, identified by a Z 2 topological invariant, is named quantum spin Hall insulator (QSHI) and it is protected by time-reversal symmetry. Nowadays, graphene and the Kane-Mele (KM) model stand as one of the archetypal time-reversal invariant topological insulators, although negligible relativistic effects in carbon open only a vanishingly small band gap in graphene [5].
Notably, the KM model still applies to all Xenes [6], i.e. the two-dimensional (2D) unary honeycomb materials made of group IV elements (e.g. silicene, germanene and stanene), where the presence of heavier atoms should lead to sizable band gaps driven by KM SOC [6][7][8]. Recently, monolayers of jacutingaite (Pt 2 HgSe 3 ) have also been proposed as novel QSHIs that display the Kane-Mele physics, and at a much larger energy scale than in Xenes, with a band gap estimated to be around ∼ 0.5 eV [9]. Jacutingaite is a naturally-occurring layered mineral first discovered in 2008 [10] in a Brazilian mine and then synthesized in 2012 [11]. Although jacutingaite is a ternary material with several differences with respect to the Xenes, it shares a (buckled) honeycomb structure of mercury atoms (see Fig. 1a), which is ultimately responsible for the KM physics in monolayers [9] that sparked experimental [12] and theoretical [13][14][15][16] interest in this material.
Being a stacking of 2D QSHIs, bulk jacutingaite is expected to be a 3D weak topological insulator with indices (0; 001), and thus with no surface states on the (001) surface [17]. Recent first-principles simulations confirmed this weak topological classification [14][15][16], but at the same time surprisingly predicted the presence of basal surface states associated with a non-trivial mirror Chern number, thus promoting bulk jacutingaite to a dual topological material [18,19] with both weak and crystalline topological properties. Such (001) surface states have now been demonstrated independently through angleresolved photoemission spectroscopy (ARPES) experiments on synthetic jacutingaite single crystals [16]. The unexpected dual topology of bulk jacutingaite cannot be understood through the standard paradigm of weak topological insulators [17] and its interpretation opens interesting perspectives on non-trivial extensions of the KM model to describe 3D stacks of honeycomb layers.
In this Letter, we show that the non-trivial topology in bulk jacutingaite emerges from a strong interlayer hybridization that leads to a 3D generalization of the KM model including a large peculiar second nearest-layer hopping term, while nearest layers are almost decoupled. Within this picture, even and odd layers are approximately independent and can be separately described by a 3D KM model where the novel hopping term drives a band inversion, giving rise to a nodal line that is gapped by SOC and a non-zero Chern number. Remarkably, when a coupling between even and odd layers is restored, the Chern numbers add up to a non-trivial value while the Z 2 classification becomes weak, thus providing a microscopic understanding for the emergent dual topology of this material. Mercury, platinum, and selenium atoms are represented in green, grey, and yellow, respectively. b-c) Band structure of monolayer (b) and bulk (c) jacutingaite along a high-symmetry path obtained using density-functional theory with (green) and without (orange) spin-orbit coupling (SOC). The inset shows the bulk and surface Brillouin zone. In monolayer jacutingaite the low-energy physics is described by graphene's Kane-Mele model (KM), with a substantial gap (0.15 eV with density-functional theory and ∼ 0.5 eV with many-body perturbation theory at the G0W0 level) that opens at the Dirac point due to SOC. In 3D jacutingaite, Dirac cones are clearly visible at K and H and they are weakly gapped by KM-type SOC; the band structure shows a strong dispersion along the stacking axis.
We first shortly review the band structure of monolayer Pt 2 HgSe 3 (shown in Fig. 1b) and its relation to the KM model, as this will be instrumental to unveil the emergent topological properties of bulk jacutingaite crystals. When SOC is neglected, the valence and conduction bands touch at the Fermi energy, forming Dirac cones located at the corners K/K of the 2D Brillouin zone (BZ). In a basis of maximally-localized Wannier functions [20] centered on mercury atoms [9], the linear dispersions arise from a hopping term between nearest neighbors on the (buckled) honeycomb Hg (sub)lattice, similarly to what happens in graphene. The inclusion of SOC opens a substantial gap at K/K , turning the system into a QSHI. As discussed in Ref. 9, gap opening mainly stems from a complex-valued second-nearest neighbor hopping, exactly as proposed by Kane and Mele, so that all qualitative features of the band structure of monolayer jacutingaite can be understood in terms of the KM model [3,4,21]: where the sums are restricted to pairs ij ( ij ) of first (second) nearest neighbor sites i and j, v ij , u ij , d 0 ij are geometrical parameters [22], and s = (s x , s y , s z ) are spin Pauli matrices. Here, in addition to the original KM hopping amplitudes t and ∆ associated respectively with the nearest-neighbor hopping and the KM SOC, we are adding an "in-plane" SOC term with amplitude ∆ , which is not present in planar honeycomb lattices such as graphene, but that appears when in-plane mirror symmetry is broken [21] as is the case in monolayer jacutingaite. Extending naively the analogy between graphene's KM model and monolayer jacutingaite to 3D, one would expect the electronic properties of bulk jacutingaite to be almost identical to its 2D form, as it is the case between graphene and graphite. In Fig. 1c we report the band structure of bulk jacutingaite computed along a highsymmetry path by density-functional theory (DFT) [23]. Without SOC, a linear dispersion is observed close to the K and H points, which is indeed reminiscent of the 2D Dirac cones. Still, the linear behavior extends over a much larger energy range than in 2D and a remarkable energy dispersion (compared to the overall band width) appears along the vertical direction between H and K. Even more compelling, SOC opens a gap between valence and conduction bands at K (and H) as in 2D, but the magnitude of the splitting is 1-2 orders of magnitude smaller than in the monolayer limit. Overall, these features suggest a significant coupling between layers, which is consistent with the non-negligible interlayer binding energy reported in Refs. 9 and 24 (∼ 60 meV/Å 2 , compared to 20 meV/Å 2 for graphene) that sets jacutingaite as "potentially exfoliable" [24].
To show that this significant interlayer coupling is responsible for the unexpected emergence of the bulk topology from the properties of the monolayer, we develop a minimal tight-binding model that captures all the relevant physics by extracting the most important hopping terms from the complexity of the full electronic structure. We first build a 4-band model (including spin) us- IG. 2. Wannier functions and even/odd layers decoupling. a) Top and lateral views of the maximally-localized Wannier functions (MLWFs) underlying the minimal model for bulk jacutingaite (see text). A reference MLWF is shown together with two additional ones (with lighter colors) to highlight the initial and final states of two relevant hopping terms in the model Hamiltonian: a second-nearest neighbor hopping process in the same layer (with two possible paths marked with solid or dashed lines) and a second-nearest layer hopping term (dash-dotted line) that drives bulk jacutingaite into a crystalline topological phase. b) Schematic of the most relevant hopping terms in the effective 1D model describing bulk jacutingaite at fixed parallel momentum k = (kx, ky). The strongest terms are the intralayert(k ) and second-nearest layert2(k ), whilet1(k ) is almost negligible. This effectively leads to a decoupling between even and odd layers, which behave as a k -dependent Su-Schrieffer-Heeger-like chain with a doubled unit cell (orange dashed line) with respect to the primitive one (black solid line).
ing Hg-centered maximally-localized Wannier functions (MLWFs) [20] that reproduces the main features of the band structure around the Fermi level [23], such as the presence of SOC-gapped Dirac cones and the band dispersion between H and K. The corresponding Wannier functions are plotted in Fig. 2a: they clearly resemble the ones obtained for the monolayer in Ref. [9], but a notable difference is due to the spatial extension of the MLWFs; those of bulk jacutingaite spread over a neighboring layer, so that they are effectively localized on 2 layers.
We now examine the strength of the hopping terms that appear in the MLWF Hamiltonian. The strongest term is the in-plane nearest-neighbor hopping t that is responsible for the linear dispersion close to K and H, with a small splitting between valance and conduction bands which is opened by a very weak KM SOC. The fact that for bulk jacutingaite the effective KM SOC, as well as the in-plane SOC [9,21], are strongly renormal-ized with respect to the monolayer can be understood by looking at the geometry of the MWLFs. Indeed, as shown in Fig. 2a, in bulk Pt 2 HgSe 3 the overlap between MLWFs allows two alternative paths to hop to secondnearest neighbors within the same layer: one identical to monolayer jacutingaite (solid line) and one extending through the closest layer (dashed line). Owing to the different sign of the geometrical parameters v ij and u ij (see Eq. (1)), the two paths give opposite contributions and, being very similar in magnitude, result in a very weak KM (and in-plane) SOC compared to the monolayer.
Remarkably, the second strongest hopping, which is comparable in magnitude to t, is a second nearest-layer hopping term t 2 that connects two MLWFs as shown in Fig. 2a (dash-dotted line). Although the MLWFs that are involved are relatively far apart, the strong hopping amplitude stems from the partial delocalization of the MLWFs over the neighboring layer, which gives rise to a large overlap between MLWF that are two layers apart through the intermediate layer. Owing to the3m symmetry and the geometrical arrangement of Wannier functions, the strong overlap takes place between a reference MLWF and three others located two layers above (or below) in the opposite sublattice, thus mimicking nearestneighbor hopping but with a doubled in-plane separation [23]. As we are going to see, the inclusion of this single term in the KM model is sufficient to understand the appearance of the nodal line and the surface states.
After having identified the most important hopping terms from the MLWF Hamiltonian, we build the following tight-binding model: where H KM is the KM model of Eq. (1),H 2 nd N L is the second nearest-layer hopping just mentioned, and λ is a dimensionless coupling constant that interpolates between graphene or monolayer jacutingaite (λ = 0) and bulk jacutingaite (λ = 1). This model is a generalization of the Kane-Mele model to 3D jacutingaite and hereafter we will call it J3KM model. As shown in Fig. 2b, at this level even and odd layers are completely decoupled, so it is natural to choose a unit cell that is doubled along the stacking axis and to describe separately even and odd layers through Eq. (2). In this way,H 2 nd N L becomes a first-nearest layer term in the even/odd subspace, and the BZ is halved along k z . For the sake of simplicity, we start by considering the J3KM model without SOC, corresponding to an Hamiltonian with only two terms: first nearest-neighbor (as in graphene) and inter-layer hopping. The band structure for a 30-layer slab is reported in Fig. 3a for different values of λ, where states localized at the (001) surface are colored in red. For λ = 1 (corresponding to bulk jacutingaite), in addition to the graphene-like Dirac states at K, such simple model accounts for the presence of additional linear crossings between valence and conduction bands (e.g. at a low-symmetry point between Γ and M), in perfect agreement with current ARPES measure- ments [16]. Notably, the model also exhibits 001-surface states that roughly span the same BZ region as observed in experiments [16], thus providing a remarkably realistic qualitative description of the system, despite its simplicity. An interesting feature of the J3KM model that helps to rationalize the presence of surface states is that, at a given parallel momentum k = (k x , k y ), it is equivalent to a k -dependent 1D tight-binding Hamiltonian analogous to the Su-Schrieffer-Heeger (SSH) one, where the alternating hopping energies,t(k ) andt 2 (k ), between intralayer and interlayer neighboring sites of a Hg chain depend parametrically on k (see Fig. 2b and [23]). The polarization of this effective 1D chain can be written as P (k ) = eγ(k )/(2π) in terms of the Zak phase [26,27]: where u(k , k z ) is the periodic part of the occupied oneelectron eigenstate of the J3KM model at zero SOC. The combination of time-reversal and inversion symmetry dictates that the Zak phase can only assume two topologically distinct values: γ(k ) = π and γ(k ) = 0, depending on the relative strength oft(k ) andt 2 (k ) as in the SSH model. According to the surface charge theorem [27,28], the chain has an end charge whenever γ(k ) = π (polarization e/2), while it is topologically trivial when γ(k ) = 0. We thus expect surface states in the J3KM for all values of k for which γ(k ) = π. Indeed, as shown in Fig. 3a and b, surface states in bulk jacutingaite (λ = 1) appear precisely in regions of the BZ where the Zak phase is non-trivial. The analogy with the SSH model also reveals that, along the lines that separate topologically distinct regions of the BZ we need to have |t(k )| = |t 2 (k )|, so that the gap closes for some value of k z . In other words, those lines can be considered as projections on the (k x , k y ) plane of nodal lines where the gap between valence and conduction bands vanishes. Indeed, by computing the energy bands of the J3KM model without SOC we uncover the presence of a nodal line (see Fig. 3c) dispersing across the border of the reduced BZ -also predicted by firstprinciples simulations [15,16]-whose projection on the parallel plane is consistent with the boundary between regions with topologically distinct Zak phases [29].
We can now understand the emergence of surface states and nodal lines in the J3KM model by studying the evolution as a function of the interlayer coupling λ. Fig. 3a and b show the slab band structure and the Zak phase computed at different values of λ. For small λ, the band structure is essentially graphene-like, with no surface states and a trivial Zak phase over the full 2D BZ.
With increasing λ, the occupied and empty bands get closer and closer to each other, until a band inversion occurs at the time-reversal-invariant point L of the reduced BZ (corresponding to M in 2D). The band inversion creates three inequivalent nodal lines, whose projections separate regions with different Zak phases. Correspondingly, surface states appear in the slab calculation wherever γ(k ) = π. With a further increase in λ, the three nodal lines merge into a single one, as shown in Fig. 3c for bulk jacutingaite (λ = 1). The interlayer coupling thus plays a crucial role in driving the essentially trivial electronic structure of weakly coupled layers into the rich physics of bulk jacutingaite, with surface states associated with a nodal line in the absence of SOC.
We now include SOC to show the robustness of surface states and their topological protection within the J3KM model. First, we consider the KM SOC [3,4] only, which gaps the nodal line almost everywhere, but not on the intersection with the vertical plane containing the Γ-M line (and its 3-fold rotation symmetric partners). The inclusion of also the in-plane SOC [6,21] fully gaps the residual Dirac points and the system becomes a topological crystalline insulator, as supported by calculations of the mirror Chern number (see Fig. 2d) providing C M = −1.
As shown in Fig. 3e, the non-trivial Chern number protects the presence of 001-surface states even when SOC is included, with a Dirac-like dispersion close to the M point.
Further calculations of the strong Z 2 invariant ν show that the J3KM model for the even/odd subspace actually describes a strong Z 2 topological insulator, in agreement with the fact that ν ≡ C M mod 2 in this space group. When considering together even and odd layers, the Z 2 invariant of the two subspaces is summed and becomes trivial (1 + 1 ≡ 0 mod 2), while the mirror Chern numbers add up to C M = −2. On one side, this means that the weak Z 2 topology of bulk jacutingaite does not fit the standard paradigm of weakly coupled 2D QSHI, but it is intimately related to the even and non-zero Chern number and to the double band inversion (one for each even/odd subspace) driven by the strong interlayer coupling. On the other, we expect the mirror Chern number to protect the surface states also when both subspaces are considered together.
In order to support this conclusion and at the same time to provide further evidence that the above predictions are not related to some extra symmetries (e.g. particle-hole) of the simple J3KM model, we finally consider the full MLWF Hamiltonian, which includes in particular additional terms to Eq. (2) that: (i) restore the coupling between even and odd layers; (ii) introduce a finite dispersion along the K-H line as in the first-principles results of Fig. 1c; and (iii) break particle-hole symmetry, making the system a compensated semimetal. Still, they do not affect the topological classification of bulk jacutingaite. The corresponding band structure for a 60-layer slab is reported in Fig. 3f. Consistently with C M = −2, two pairs of surface states (degenerate at M) are present, slightly split by the coupling between even and odd layers. Remarkably, these bands are very similar to what is observed in ARPES experiments [16].
In conclusion, we provide a microscopic insight on how symmetry-protected topological order in layered jacutingaite emerges from a non-trivial coupling between Kane-Mele-type QSHI monolayers. The essential physical features can be captured by a simple generalization of the Kane-Mele model to account for interlayer hopping. This J3KM model predicts the presence of surface states and nodal lines gapped by spin-orbit interactions, in remarkable agreement with recent ARPES measurements and first-principles simulations, providing an appealing strategy to break the standard paradigm of weak topological insulators that becomes relevant for all other layered materials made of stacked honeycomb lattices.

S1. FIRST-PRINCIPLES SIMULATIONS
Density-functional theory calculations are performed with the Quantum ESPRESSO distribution [35,36], Wannier functions are obtained using WANNIER90 [37]. Structural optimization is performed by using the experimental lattice parameters as obtained from X-ray diffraction and relaxing the atomic coordinates with a non-local van der Waals functional, namely the vdW-DF2 functional [38] with C09 exchange (DF2-C09) [39], and the SSSP precision pseudopotential library v1.0 [40] with 100 Ry of wavefunction cutoff and a dual of 8. Further calculations on the optimized crystal structure (band structures and Wannier functions) are performed using the PBE functional [41] and ONCV [42] scalar and fully relativistic pseudopotentials from the PseudoDojo library [43] with 80 Ry of wavefunction cutoff and a dual of 4. All calculation are perfomed with k-point density of 0.09Å −1 , that corresponds to a k−point grid of 12 × 12 × 14, and Marzari-Vanderbilt smearing [44] of 0.015 Ry. Wannier functions are constructed from a k−point grid of 6 × 6 × 6. Part of the calculations are powered by the AiiDA [45] materials' informatics infrastructure.

S2. 4-BAND MODEL FOR BULK JACUTINGAITE
In order to construct a minimal tight-binding model for bulk jacutingaite, we follow a similar strategy as for the monolayer [9] by mapping first-principles calculations onto a set of maximally-localized Wannier functions (ML-WFs) [20], constructed from an initial projection on Hg-centred s-like orbitals. This results in a 4-band (including spin and relativistic effects) tight-binding model on a lattice made of buckled honeycomb layers directly stacked on top of each other (AA stacking). We denote the lattice vectors of the unit cell with a 1 = (a, 0, 0), a 2 = (−a/2, √ 3a/2, 0), and a 3 = (0, 0, c), with the two inequivalent sublattices centered at τ A = 2/3a 1 +1/3a 2 +δa 3 and τ B = 1/3a 1 +2/3a 2 −δa 3 . The band structure obtained from this 4-band model is compared to the original first-principles results in Fig. S1. Although a perfect quantitative agreement is not attainable, the model reproduces all qualitative features around the Fermi energy, including in particular the extended linear dispersion and the very small band gap between valence and conduction bands at K and H.

S3. SECOND-NEAREST-LAYER HOPPING TERM
Within the 4-band model, the largest hopping term involves MLWFs that are centered on neighboring sites in the same layer, with amplitude t = 0.27 eV. As mentioned in the main text, the next largest contribution is a secondnearest layer hopping term with amplitude t 2 = −0.18 eV −0.7 t. This hopping process involves a reference Wannier function on the A (B) sublattice and one of three B (A) sites that lie 2 layers above (below) in unit cells identified by the lattice vectors In Fig. S2 we sketch the sites involved in this hopping process, starting either from the A sublattice (red solid arrows, upper signs in Eq. (S1)) or the B sublattice (blue dashed arrows, lower signs in Eq. (S1)). This hopping term has been adopted in the main text to construct a minimal extension of the 2D Kane-Mele model that is able to describe the emergent topology in bulk jacutingaite. The shaded yellow area shows the primitive unit cell, with lattice vectors a1, a2, and a3. Arrows connect sites involved in the second-nearest layer hopping process responsible for the non-trivial topological properties of bulk jacutingaite, starting from either the A sublattice (red solid arrows, corresponding to the vectors τB + R1,2,3 − τA with the upper sign in Eq. (S1)) or the B sublattice (blue dashed arrows, corresponding to the vectors τA + R1,2,3 − τB with the lower sign in Eq. (S1)).
The condition for the gap to close reads |t(k )| = |t 2 (k )|, which has non-trivial solutions only when |t 2 /t| > 1/3. At the values of k for which this condition is satisfied there must exist a value of k z at which the gap closes, thus giving rise to a nodal line. We have verified that, as expected, the projection of the nodal line on the (k x , k y )-plane (as identified by the solution of |t(k )| = |t 2 (k )| coincides with the line separating regions where the Zak phase (computed from the eigenstates of (S2)) is trivial and regions where it is non-trivial.