Topology of SmB$_6$ determined by dynamical mean field theory

Whether SmB6 is a topological insulator remains hotly debated. Our density functional theory plus dynamical mean field theory calculations are in excellent agreement with a large range of experiments, from the 4f5.5 intermediate valency to x-ray and photoemission spectra. Using the pole extended (PE) Hamiltonian, which fully captures the self-energy, we show that SmB6 is a strongly correlated topological insulator, albeit not a Kondo insulator. The PE Hamiltonian is proved to be topological (in)equivalent to the topological Hamiltonian for (non-)local self-energies. The topological surface states are analyzed, addressing conflicting interpretations of photoemission data.

On the experimental side, there is clear evidence of robust metallic surface states, such as a surface-dependent plateau in the low-temperature resistance [16,17] and angular resolved photoemission spectroscopy (ARPES) data [18][19][20][21] including its spin texture [22]. However, their topological origin has been questioned [21,23]. Indeed, whether the low-energy electronic properties even originate from the surface or the bulk remains hotly debated. De Haas-van Alphen oscillations have been interpreted as stemming from both, the bulk [24] and the surface [25,26]. The same holds for the main ARPES features: [27] vs. [11,26]. The low-temperature linear specific heat has been shown to be predominantly a bulk effect [28]; and the same holds for the optical conductivity within the band gap [29]. On the other hand, it is was argued [30] that these effects are not intrinsic but stem from 154 Gd impurities which eludes a mass purification against 154 Sm.
This controversy clearly calls for a better theoretical understanding. Because of the strong electronic correlations and well-localized 4f orbitals DFT + Dynamical Mean Field Theory (DMFT) [31][32][33] is the method of choice. There have been earlier DFT+DMFT [34][35][36], DFT+Impurity [37,38], and DFT+Gutzwiller[39] calculations which however did not capture the bulk band gap, the flatness of the f -bands, and the intermediate valence of SmB 6 all at the same time [11], due to various additional approximations. Indeed, the combination of these properties pose a hard theoretical challenge that requires an accurate many-body treatment of the Sm 4f orbitals.
In this letter, we present charge self-consistent DFT+DMFT calculations for SmB 6 . Our results provide an all encompassing picture of the bulk and surface properties of SmB 6 , in excellent agreement with many experimental observations. Analyzing the symmetry properties of the corresponding PE Hamiltonian [14,15] we show that SmB 6 is a strongly correlated topological insulator. Furthermore, we prove that the "topological Hamiltonian" [12] is equivalent to the PE Hamiltonian if the self-energy is local, and pin-point the reason the former may fail for momentum-dependent self-energies.
DFT+DMFT Method. All calculations have been performed using the relativistic spin polarized toolkit (RSPt) [40][41][42][43], which is based on linearized muffin tin orbitals. This method allows the correlated Sm 4f orbitals to be readily identified and projected upon in the DMFT calculation [41]. The full Coulomb interaction, spin-orbit interaction, and local Hamiltonian for the Sm 4f orbitals are included in the DMFT impurity problem, which is solved using the exact diagonalization (ED) approach [42,44,45]. The calculations were performed at  [11]. Sm 4f 6 → 4f 5 transitions with the J and L quantum numbers of the final (4f 5 ) state are indicated. Lower panel: The trace of the local self-energy (real part) shows distinct poles, but not within the bulk band gap. T = 100 K and iterated until charge self-consistency [46]. For additional details, such as the final bath state parameters and the double counting procedure [47,48], see the Supplemental Material (SM) [49].
Bulk. The hybridization between the strongly localized Sm 4f orbitals and the surrounding orbitals is too weak to form a coherently screened (Kondo) ground state at 100 K. Instead we find an intermediate valence state with a thermal mixture of both Sm 4f 6 and 4f 5 configurations, with an average 4f occupation of n f ≈ 5.5, in good agreement with experiment [50]. The 4f 6 contribution is predominately of Γ 1 character and the 4f 5 of Γ 8 in agreement with non-resonant inelastic x-ray (NIXS) data [51]. A detailed characterization of the thermal ground state is given in Table SII in the SM [49].
The k-integrated spectral function (DOS) reflects the intermediate valence and displays distinct 4f 6 → 4f 5 and 4f 5 → 4f 4 multiplet transitions, as shown in Fig. 1 and the SM [49], respectively. The first peak below E F at -11 meV is to a 6 H 5/2 Γ 8 final state [49]. Upon closer inspection, we notice a Γ 7 subpeak with a maximum at -25 meV, too small to be well resolved in the experiment. The next peak around -170 meV corresponds to J = 7/2 final states ( 6 H 7/2 ), and the on-resonance spectrum ( ν = 140 eV) displays an additional J = 9/2 peak around -300 meV. The latter is hardly discernible in theory and at off-resonance ( ν = 70 eV) as the direct transition is largely forbidden.
These sharp multiplet peaks are generated from a many-body self-energy which has distinct poles, as seen in Fig. 1 (lower panel). However, there is no pole within the narrow bulk band gap, as one would have for a Mott insulator. We will come back to the poles of the selfenergy when discussing the topological properties.
Let us now turn to the momentum resolved spectral function in Fig. 2  on the other hand hybridize strongly with the B 2p orbitals and form a parabolic band centered at the X-point. When this dispersive band crosses the flat f -bands close to the E F their hybridization leads to a small band gap of about 9 meV (16 meV peak-to-peak), close to the experimental value of ∼10 -20 meV [11,52]. The spectral function agrees well with the bulk sensitive ARPES experiments of Denlinger et al. [11] reproduced in Fig. 2.
Altogether, we find that the local DFT+DMFT selfenergy gives an accurate description of the experimentally established bulk properties of SmB 6 . Hence we can now turn to its topological properties with confidence.
Proof of non-trivial topology. The topological Z 2 invariant of Kane and Fu [10] can be determined, for an inversion-symmetric non-interacting system, from the parity P ki;m of the (pairs of Kramers degenerate) bands m below E F at the time reversal invariant moment (TRIM) momenta k i . If im P ki;m = −1, the system has a non-trivial topology. However, this procedure can not be directly applied to an interacting system, as the one-particle self-energy can split, smear out, and reduce the weights of the bands, and even make them fade away and reappear at shifted energies. Several suggestions of how to generalize the Z 2 invariant to the interacting case have been made [12,15], but as common in the fast moving field of topological materials, the theory still needs to be developed in full.
The self-energy consists in general of a set of poles located on the real frequency axis [49,53], clearly shown in Fig. 1 (bottom). The local self-energy can hence be written as Here, E l and V † ml V ln are the pole position and weight, respectively; and mn are the local spin-orbital indices. This pole structure of the self-energy is identical to that of a hybridization function where the physical orbitals m hybridize via V lm with some local auxiliary orbitals l. That is, we can make an exact mapping, akin to a purification of a mixed state, of the interacting bulk system with the self-energy Σ mn to a non-interacting "pole extended" (PE) Hamiltonian with additional orbitals l [14,49]. As in the case of a purification, we recover the physical band structure by simply projecting the PE band structure onto the physical orbitals. As this PE Hamiltonian is non-interacting, the topological invariant Z 2 of [10] can be straight forwardly applied [15]. If this Z 2 is non-trivial, the PE Hamiltonian has topological surface states, which manifest themselves in the physical system through the projection onto the physical orbitals.
In the following, we show that SmB 6 is a topological insulator by proving that (1) its PE Hamiltonian is topologically non-trivial and that (2) the resulting topological surface states have a finite physical weight.
As for (1), the construction of the PE Hamiltonian requires the full pole structure of the self-energy, which can be difficult to obtain. Nevertheless, we can still assess its topological properties from some general considerations: (i): As shown in Fig. 1, the self-energy has no pole at the Fermi energy.
(ii): The DMFT self-energy is k-independent in its local basis. The PE bands must therefore have a finite physical weight to be dispersive.
From (i) and (ii) it follows that since there is a finite band gap in the known physical spectral function, the PE Hamiltonian is insulating as well [54].
(iii): Only the correlated Sm 4f states carry a DMFT self-energy, the Sm 5d orbitals and B 2p orbitals forming the dispersive bands have Σ = 0. This is important since the auxiliary orbitals given by a local self-energy inherit the irreducible representation of the local correlated orbitals they derive from [55]. In our case, the auxiliary orbitals have the same odd parity [56] as the Sm 4f orbitals. It is therefore enough to count the pairs of the non-interacting bands of opposite (even) parity at the TRIM points to evaluate the Z 2 invariant! These are 6, 6, 5, and 4 at Γ, R, X, and M , respectively [57]. The physical reason for the odd number at the X point is that the dispersive band of Sm 5d and B 2p character is below E F at X. This implies im P ki;m = −1, i.e., the PE Hamiltonian is topologically non-trivial.
This non-trivial topology is very robust against potential imprecisions in our DFT+DMFT calculation: they would need to shift an even parity orbital to the other side of E F at an odd number of TRIM points, which requires a shift by several eV, while maintaining the mixed valency. Furthermore, this shift must originate from the DFT potential and hence the electron density since Σ = 0 for these orbitals.
As for (2): In general one expects that a non-trivial topology in the bulk implies robust metallic surface states. For a k-dependent surface self-energy this need not be the case: In a gedanken experiment we can simply craft a surface self-energy using auxiliary orbitals that mimics a large number of additional layers to the vacuum side. Such a self-energy makes the surface layer indistinguishable from a bulk layer. Hence the topological surface states are completely absorbed by the surface self-energy. However, for a local DMFT self-energy, which is a good approximation for SmB 6 , a complete absorption is not possible. The local self-energy may vary from layer-to-layer (perpendicular to the surface), but it still fulfills (ii). Hence, the dispersive topological surface states must carry a finite physical weight as they cross the bulk band gap.
Let us also address the concerns of a possible break down [13] of the "topological Hamiltonian" [12]. For a kindependent self-energy such a breakdown is not possible as long as it is defined, i.e., (i) Σ(0) < ∞. We show this in the SM [49] by adiabatically detaching the auxiliary orbitals from the PE Hamiltonian. The resulting Hamiltonian is identical to the "topological Hamiltonian", apart from the detached auxiliary orbitals, which due to their purely local character do not contribute to the non-trivial topology. Hence, the "topological Hamiltonian" and the PE Hamiltonian are topologically equivalent when the self-energy is local. For a k-dependent self-energy, on the other hand, the auxiliary orbitals can contribute in a non-trivial way to the elementary band representations [6,49]. This can manifest itself as additional topological surface states that gap out the topological surface states predicted by the topological Hamiltonian, which is consistent with the break down seen in [13].
Surface states. To address the conflicting interpretations of the experimental photoemission data [11,20,21], we performed additional DFT+DMFT slab calculations. The [001] SmB 6 supercell, shown to the right in Fig. 3, was used with two different terminations to represent the surface patches of sputtered films[58] and cleaved samples [21]. The former has a terminating B 6 layer (depicted without bonds) and the outermost Sm atom is Sm 3+ , while the latter lacks the B 6 layer and is instead terminated by Sm 2+ surface atoms. The B 6 termination has trivial surface bands associated with a B 6 dangling bond [23] that disappears when the surface reconstructs [11,21,23]. To effectively mimic this partial surface reconstruction in our computationally expensive DMFT calculation, we simply apply an additional 0.5 eV potential to the 2p orbitals of the outermost B atom in the final step. This potential shifts the trivial surface band, but does not directly affect the topological surface states which live, as we will see, in the subsurface layer.
The resulting DFT+DMFT spectra are presented in Fig. 3 (left). We recover the flat Sm 4f bands of the bulk at -14meV and +40meV as well as the m J = ±1/2, J = 5/2 band directly above E F at the X point. On top of these bulk states surface bands emerge within the bulk band gap. These metallic surface bands are mainly  4 shows cuts of the spectral function at E F and -5 meV below. The surface states form large and intermediate sized pockets around both the X and Γ points for the B 6 and Sm termination, respectively. The X pockets show a strong spin polarization in agreement with spin-resolved ARPES experiments [20,49]. Ref. [21] reports both the large X and Γ pockets for the B 6 terminated surface, although they interpret the large Γ pocket as an umklapp state. The trivial "Rachba split" surface states reported close to the Γ point directly after cleaving [21] seem instead to correspond to the shifted B 6 derived trivial surface bands. Ref. [11] shows the intermediate sized X pockets and the H features at -5 meV associated with the Sm terminated surface, reproduced in Fig. 4, as does [18,20], but the Γ pocket is seemingly missing. However, on closer inspection of the data, in particular in the second Brillouin zone, there is clear evidence of a matching surface derived Γ pocket, which again has been interpreted as an umklapp state [18,20]. Hence, to finally settle the apparent discrepancy between theory and experiment, we suggest a critical experimental reexamination of this umklapp assignment. For example, our data for the Sm termination suggest that the Γ pocket has a much weaker spin-polarization than the X pocket.
Conclusion. The appropriateness of DFT+DMFT for describing the strong correlations in SmB 6 and the excellent description of various experimental properties gives us confidence that we have achieved an accurate theoretical description of SmB 6 . We determined the topological nature in a rigorous way from the symmetry properties of the PE Hamiltonian, formed by mapping the pole structure of the self-energy onto a set of auxiliary orbitals [14,15]. We made the key observation that also the auxiliary orbitals contribute to the "elementary band representations" [6] of the PE Hamiltonian, which clarifies the fundamental topological role of the poles of the selfenergy. Nevertheless, we prove that for a local DMFT self-energy it is enough to count the bands of the opposite parity channel at the TRIM points, under the condition that the self-energy does not have a pole at E F . This band analysis shows that SmB 6 constitutes a topological insulator.
In general, if the positions of the bands in one parity channel is well-described by DFT, then the addition of a local self-energy in the other parity channel will not change the topology of the system as long as it remains insulating. If the local self-energy has poles within the bulk band gap topological surface bands will still appear in the spectral function, but they may gradually flatten out and vanish as they approach the poles.
We have also shown, through a simple qualitative ar-gument, that a non-local symmetry preserving surface self-energy can completely absorb any topological surface bands. However, our data suggests that the "missing" Γ pocket in SmB 6 instead is simply disguised as umklapp states [18,20,21].
We are highly indebted to J. D. Denlinger and J. W. Allen for making available their ARPES data. We would also like thank Annica Black-Schaffer and Dushko Kuzmanovski for useful discussions. This work has been financially supported in part by European Research  [54] A purely auxiliary band can not cross the band gap as this requires a finite dispersion.
[55] The reason is simple, the symmetries of the system enforce that only local orbitals that share the same irreducible representation may hybridize.
[56] When the inversion center is placed at the Sm atom.
[57] Excluding the semi-core and core states. [59] The physical and auxiliary orbitals form two separated non-trivial subspaces when V = 0. However, since V = 0 is not considered a protected symmetry, this separation is not topological and the full system is still trivial.
The topology of SmB 6 determined by dynamical mean field theory P. Thunström 1,2 and K. Held 1 (Dated: July 10, 2019) In this supplemental material, we first show that the pole expansion Hamiltonian is the proper description of the topological properties in many-body systems and that it is equivalent (not equivalent) to the so-called topological Hamiltonian for a local (non-local) self-energy. Second, we discuss computational details of our density functional theory plus dynamical mean field theory (DFT+DMFT) calculations, using an adapted exact diagonalization as an impurity solver. Third, additional results of the DFT+DMFT electronic structure are presented. Fourth, we demonstrate the robustness of the topological surface states, unless time-reversal symmetry is broken by a magnetic field.

POLE EXPANDED HAMILTONIAN
Local self-energy In the following we will show that for a local (e.g. DMFT) self-energy the pole extended (PE) Hamiltonian [14] and the topological Hamilton [12] are topologically equivalent. The proof employs some DMFT concepts explicitly, which can also be used whenever we have a local self-energy (we can just define an impurity model from that local self-energy and the non-interacting Green's function). In the end of this Section, we will then proof that they are in general not equivalent if non-local selfenergies are permitted.
In DMFT each correlated lattice site is mapped to a single impurity Anderson model and the resulting selfenergy is mapped back to the lattice. The impurity model consists formally of a set of local correlated orbitals, e.g. the Sm 4f orbitals (F ), and a set of noninteracting bath orbitals (B). The one-particle term H of the impurity Hamiltonian can hence be partitioned as where the subscripts denote the orthogonal projections upon F and B, i.e. H F B ≡ F HB where F and B are the corresponding projection operators. Together, the full projection P = F + B spans the entire system at hand, i.e. H P P = H. The Lehmann representation of the interacting local Green's function shows that it can be written as a (large) sum of simple poles where E l and U ml U † ln give the pole positions and their weights, respectively, and mn runs over all orbitals in the system (P). The poles E l originate from energy differences between two many-body states, and the (rectangular) matrix U fulfills l U ml U † ln = δ mn because of the normalization of the Green's function.
The latter property allows us to decompose the matrix U as U = P W , where W contains the eigenvectors of a large auxiliary HamiltonianHW = W E and P is the projection onto the (physical) system as before, with the eigenvalues E l forming the diagonal matrix E. Substituting U = P W into Eq. (S.2) yields the following expression for the local Green's function Eq Substituting the operator identity ("downfolding") Hence, the dynamical part of the self-energy,H P A (ν − H AA ) −1H AP , corresponds to a hybridization function to the auxiliary orbitals. The projection in Eq. (S.7) is onto both the impurity orbitals and the bath P = F + B. However, the uncorrelated bath orbitals do not have a self-energy, i.e. BΣ(ν) = Σ(ν)B = 0, which can be shown e.g. from the equation of motion or by integrating out the non-interacting bath orbitals. This implies thatH P P = H P P + Σ F F (∞) andH BA =H AB = 0, which in turn yields Eq. (S.8) is equivalent to Eq. (1) in the main text. In DMFT, the lattice self-energy is replaced with periodically repeated copies of the local self-energy Σ(ν). The bulk lattice Green's function G k (ν) is hence given by where k belongs to the first Brillouin zone; and the projection operator F corresponds to the Wannier projection of the lattice basis onto the local correlated orbitals. We may now use Eq. (S.6) again, but this time in the other direction ("upfolding"), to extend the Hilbert space with the auxiliary orbitals in A, where L projects upon the physical orbitals of the lattice, just as P projected on the physical orbitals (impurity + bath) in Eq. (S.3) in the context of the impurity model. The pole extended (PE) lattice Hamiltonian is hence given by While our derivation assumed a local, i.e., kindependent self-energy, Eq. (S.2) and Eq. (S.6) can also be formulated for each k-vector. The important difference is that this makes the auxiliary Hamiltonian kdependent,H k , but with this modification Eq. (S.10) still holds. Hence, both, for a k-independent and a kdependent self-energy, the physical Green's function, and thus the spectral function, is completely determined by the eigenvalues of H P E k and the projection onto the correlated orbitals.
If the non-interacting PE Hamiltonian H P E k is topologically non-trivial then there will be topological surface states. The physical weights of these surface states will depend on the hybridization between the physical and the auxiliary orbitals at the surface. Only a vanishing hybridization can potentially confine a topological surface state of H P E k to the auxiliary sector, and thus prevent it from appearing in the physical spectrum of G k (ν). However, the topological surface states must be dispersive to cross the bulk band gap of H P E k , which implies that if the self-energy is k-independent the topological surface states have to carry a finite physical weight.
so that H k (λ = 0) = H T k and H k (λ = 1) = H P E k . The corresponding physical lattice Green's function becomes

(S.13)
A key point is that interpolation is chosen in such a way that the physical Green's function evaluated at the Fermi energy, G λ k (0), is independent of λ, v such that v = Av and H k (λ)v = 0 for some k and 1 > λ > 0. However, since v has no physical (F ) com- Further-more, the auxiliary part of H k (λ) is identical to that of the PE Hamiltonian, H k (λ) AA = (H P E k ) AA , which together withH F A v = 0 gives H P E k v = 0. That is, the purely auxiliary vector v must also be an eigenvector of H P E with a zero eigenvalue. This eigenvector can not exist since the PE Hamiltonian was assumed to be insulating. This implies that no band can cross the Fermi energy during the continuous interpolation in Eq. (S.12).
The PE Hamiltonian H P E k and the combined topological auxiliary Hamiltonian H T A k are hence topologically equivalent. The topological Hamiltonian H T k is by itself only topologically equivalent to H P E k if the detached auxiliary HamiltonianH k AA is topologically trivial.H k AA is automatically trivial if the self-energy is local (k-independent), but for a non-local self-energy it becomes necessary to consider how the auxiliary orbitals contribute to the elementary band representations [6], as shown by two concrete examples in the next section.

Non-local self-energy
Let us illustrate the topological role of the non-local self-energy by constructing two simple examples for which PE and topological Hamiltonian are not topologically equivalent. From the definition of the topological Hamiltonian H T k ≡ H k + Σ k (0) it is clear that if Σ k (0) is much smaller than the band gap of H k it cannot change the topology of H T k . On the other hand, the auxiliary orbitals in the PE Hamiltonian may change its topology even with an infinitesimal coupling strength. Let us therefore start with a system governed by the topologically non-trivial Bernevig-Hughes-Zhang (BHZ) Hamiltonian [3] and define a non-local self-energy as Σ k (ν) = |V | 2 [ν + H k BHZ ] −1 , where V is the coupling strength. Such a self-energy can be obtained if we have a non-local and frequency-dependent interaction. The corresponding PE and topological Hamiltonians can be easily identified as where 1 is the identity matrix, and the second row and column in the matrix in Eq. S.16 corresponds to the auxiliary orbitals (c.f. Fig. S.1 Another example is given by a system with the same non-local self-energy but with the constant Hamiltonian

The corresponding PE and topological Hamiltonians become
The PE Hamiltonian of Eq. (S. 19) is non-trivial if |V | 2 < 1, while it becomes trivial when |V | 2 > 3. The topological Hamiltonian has just the opposite classification.
These two examples clearly demonstrates the topological role of the auxiliary orbitals and the poles of the self-energy. Even when it is possible to detach the auxiliary orbitals from the physical system, their contribution to the topology cannot be neglected.
One may make an analogy to a system with spinorbit coupling, and let the spin-up and the spin-down channel take the role of the physical and the auxiliary orbitals, respectively. Even if it would be possible to smoothly switch off the spin-orbit coupling, and hence detach the two spin channels, one can still not determine the topology of the original spin-full system (PE Hamiltonian) from the spin-up channel alone (topological Hamiltonian). Just as the topological invariants must trace over both spin channels, they must also at least formally trace over the auxiliary orbitals.
To sum up, we have proven that the PE and topological Hamiltonian are topologically equivalent for a local (e.g. DMFT) self-energy but may be topologically distinct for a non-local self-energy. Only the PE Hamiltonian properly reflects the topology of the system, whereas for a non-local self-energy the topological Hamiltonian TABLE S.I. Screened Sm 4f Slater integrals (F 0 -F 6 ) and the double counting potential (µDC) of the Sm atoms in the bulk and the Sm terminated (Sm) and B6 terminated (B6) supercells in Fig. 3. The labels start from the middle Sm atom (Sm1) and go toward the surface (Sm2 -Sm4) may yield the wrong topology.

COMPUTATIONAL DETAILS
All the SmB 6 DFT+DMFT calculations were performed using the experimental cubic crystal structure (space group: 221, prototype: CaB 6 ) with the lattice parameter a = 4.13Å. The DFT exchange-correlation functional was set to the local density approximation [8], and the Brillouin zone was sampled through a conventional Monkhorst-Pack mesh of 16 x 16 x 16 k-points in the bulk and 16 x 16 x 1 k-points in the slab calculations. The local Coulomb interaction was parameterized in terms of the Slater parameters F 0 , F 2 , F 4 , and F 6 . The Hubbard U parameter F 0 is heavily screened by the valence electrons and was set to 8.0 eV. The less screened Slater parameters F 2 , F 4 , and F 6 were instead calculated through Slater integrals at the beginning of each DFT iteration, and then scaled by 0.95, 0.97, and 1.00, respectively [43].The final self-consistent values are given in Table S

Exact diagonalization
The Sm 4f orbitals are strongly contracted compared to the Sm 5d and 6s orbitals, which makes the Sm 4f orbitals hybridize very weakly with the orbitals on neighboring atoms. The hybridization is completely captured by the local hybridization function ∆(ν) [31][32][33] that describes how the electrons propagate in the material once they leave the Sm 4f orbitals of a given atom. To put the weak Sm 4f hybridization in perspective, the Ni 3d hybridization function in NiO is about 20 times larger than the Sm 4f hybridization function in SmB 6 for comparable computational setups. The electronic structure of SmB 6 is therefore already well-described on the eV energy scale by completely neglecting the hybridization [37]. However, this approach is not accurate enough to describe the important meV energy scale close to the Fermi energy. It is particularly important that the hybridization function at the Fermi energy ∆(0) is accurately captured, as well as the Sm 4f occupation, to ensure that the Fermi energy remains in the hybridization band gap during the DFT+DMFT self-consistency cycle. The Exact Diagonalization (ED) impurity solver [42, 44, 45] takes most of the hybridization between the correlated Sm 4f orbitals and the rest of the material into account by including a limited number of effective bath orbitals in the impurity problem. However, the total number of bath orbitals that can be included in the impurity problem is severely limited by the growth of the many-body Hilbert space. In order to still get an accurate representation of ∆(0) and the Sm 4f occupation we need to go beyond the standard hybridization fitting scheme as well as modifying the double counting correction. The details thereof are described in the next two subsections.

Bath discretization
In the self-consistent DMFT scheme the lattice Green's function in Eq. (S.9) is projected onto the local Sm 4f orbitals (F ) and integrated over the Brillouin zone to yield the local Green's function G F F (ν). The hybridization function is then extracted from the inverse of G F F (ν), where H F F is the local impurity Hamiltonian in Eq. (S.1) which includes the spin-orbit coupling, the doublecounting, and all crystal field terms. Σ F F (ν) is the local self-energy given by Eq. (S.8). The ED method proceeds by fitting a few bath state parameters H F B and H BB to ∆(ν) via the model function We fit the hybridization function on the Matsubara axis with 6 bath states per Sm orbital using a conjugate gradient scheme which also takes the off-diagonal terms in ∆(ν) into account. The high-energy bath states (B ′ ), i.e. the eigenstates of H BB with energies E b ′ relatively far from the Fermi energy ( give only a perturbative contribution to the low energy physics due to the large energy cost of exciting these bath states. Their effect can be estimated by introducing the scaling H B ′ B ′ → λH B ′ B ′ and H F B ′ → √ λH F B ′ in ∆ ED (ν), which keeps ∆ ED (0) invariant, and let the scaling parameter λ → ∞. The scaling shows that the high energy bath states can be replaced to zeroth order by a static contribution where B ′′ is the low-energy complement to B ′ , i.e. B = B ′′ + B ′ . In our calculations we put the high-energy cut-off at |E b ′ | > 10w b ′ . The remaining low-energy bath states were ordered according to their weight w b , and included in B ′′ to the extent allowed by the computational resources, with a minimum of 4 bath states in total. The converged one-particle term of the ED Hamiltonian H ED 0 = H F F + H F B + H BF + H BB for the bulk calculation is presented in Fig. S.4.

Double counting correction
In the DFT+DMFT scheme, the explicit addition of a local Coulomb interaction term to the DFT Hamiltonian introduces a double counting (DC) of the electronelectron interaction. The unknown form of the screening processes in the exchange correlation functional prevents the implementation of an exact double counting correction. In particular, the screening of the local interaction between the Sm 4f electrons, implicitly described within the local density approximation, may be different to the effective (static) screening implied by the renormalized Slater parameters. Due to this ambiguity several different double counting corrections schemes have been suggested over the years, such as the Fully Localized Limit (FLL) [47] and Around Mean Field (AMF) [48]. The FLL correction removes the spherically averaged Hartree-Fock contribution of an effective atomiclike system with integer orbital occupations. The AMF correction considers instead an effective itinerant system with uniform (non-integer) orbital occupations. However, the thermal ground state of an intermediate valence compound such as Sm 6 falls outside these two scenarios as it contains several thermally occupied almost atomic-like many-body eigenstates having different number of electrons, as shown in Table S.II. The occupation of an intermediate valence ground state responds to the DC potential in Fermi-Dirac-like steps. However, the underlying assumptions of both FLL and AMF make their DC correction linear in the occupation. The two different behaviors always lead to a feed-back loop away  from intermediate valence towards integer valence. A second less sever issue is that the finite discretization of the bath states in ED causes a small mismatch between the impurity green's function and the local green's function G F F (ν) projected from the lattice. To minimize these two problems we automatically adjusted the double counting potential at each DMFT iteration to obtain the same number of Sm 4f electrons in the impurity as in the lattice. In the bulk calculation the number of Sm 4f electrons stabilizes at 5.48 at a temperature of 100 K, remarkable close to the experimental value of 5.45 [50]. In the slab calculations we chose to enforce the bulk Sm 4f occupation of 5.48 for all the Sm atoms except at the surface, to make the center Sm atom as bulk-like as possible and minimize the hybridization between the surface states located at the top and bottom of the slab.

MANY-BODY ELECTRONIC STRUCTURE Thermal ground state
The lowest energy many-body eigenstates of the selfconsistent impurity problem Hamiltonian H ED is given in Table S.II. The states have almost atomic-like Sm 4f occupations (N f ) and total angular momenta (J), but the crystal fields and the hybridization with the bath mix the different J z configurations. The thermal groundstate is an incoherent mixture of the eigenstates according to their Boltzmann weights (GS) e −βE /T r[e −βH ] at 100 K.

Extended photoemission spectrum
In Fig. 1 of the main paper, we have already shown the spectral function in a small energy range from -1.2 eV to FIG. S.5. The Sm 4f -projected DFT+DMFT spectral density of SmB6 compared to experimental on-resonance ( ν = 140 eV) photoemission data [11]. The approximate J and L quantum numbers of the final states are indicated. 0.1 eV, which resolves the J and j z multiplets for L = 5 and L = 3. In Fig. S.5 we show the same k-integrated spectral function in a larger energy window. Clearly all peaks agree well with the experimentally measured (onresonance) photoemission data [11]. The relative weight of the 4f 6 → 4f 5 transitions are accurately captured. The relative weight of the different 4f 5 → 4f 4 transitions are in reasonable agreement with the experimental data, even though we do not include the resonance effect which strongly enhances their total weight. This gives us further assurance that our DFT+DMFT calculation faithfully describe bulk SmB 6 .

ROBUSTNESS OF THE TOPOLOGICAL SURFACE STATES
Adding a surface potential To numerically confirm the topological protection of the surface bands we applied time-reversal symmetric potentials to the (sub)surface layer of the Sm terminated supercell displayed to the right in Fig. S.8. The topological surface states should survive under these perturbing potentials, which can occur for example due to surface reconstructions. Our results are shown in Figs. S.6 and S.7, where we applied the potential to the j z = ±5/2 and j z = ±1/2 spin-orbitals of the J = 5/2 Sm 4f manifold, respectively. These spin-orbitals were chosen as the surface states around the Γ-point has primarily j z = ±5/2 character, while the surface states around the X-point have mainly j z = ±1/2 character. Technically, we simply added the potential term to the self-energy of these states after DMFT convergence (a self-consistency including this potential is beyond our illustrative purposes). Both figures show the Sm 4f projection on the sub-surface (left panel) and subsub-surface states (right panel). Please remember that, as discussed in the main text, the surface layer itself has another Sm 4f valence and is insulating. The topological surface states hence appear in the sub-surface layer in the unperturbed systems.
Fig. S.6 (left) clearly shows that one of the flat 4f orbitals (i.e., the j z = ±5/2) is shifted above the Fermi energy upon increasing the j z = ±5/2 potential V 5/2 . At the same time, the topological surface states around the Γ-point shift from the sub-surface layer for V 5/2 = 0 to the subsub-surface layer at V 5/2 = 0.73 eV. There is a crossover in-between with the surface states being extended to both layers.
At the same time the topological surface states around the X-point remain on the sub-surface layer. If we instead apply a potential to the j z = ±1/2 spin-orbitals as in Fig. S.7, it is these topological surface states which start to shift to the next layer below. The X-pocket is however more robustly anchored to the sub-surface atom compared to the Γ-pocket, and does not shift completely away even for V = 0.73 eV.
The larger "mobility" of the Γ pocket in the j z = ±5/2 J = 5/2 spin-orbitals is interesting, as it is seen in some experiments [18,20,21] but not in others [11,18,20]. In the latter experiments, the Γ pocket might simply hide a few layers deeper in the bulk, escaping its detection in surface-sensitive photoemission experiments.
Another possible explanation of the experimental discrepancies reveals itself at V 5/2 = 0.36 eV in Fig. S.6. At this potential strength the topologically derived Xpocket and Γ-pocket get very close to each other around the Fermi energy, and the latter can easily be misidentified as an umklapp-state [20,21]. Interestingly, at the very same potential there are trivial bands which grace the Fermi energy close to the Γ-point. These trivial bands may be connected to the weak α-pocket seen in some experiments [18,20,21]. However, to put these observations on firmer ground several fully self-consistent slab calculations with various realistic surface reconstructions are needed, which is beyond the scope of current study. In the middle panel, the field is applied to the jz = ±1/2 states, and a smaller gap emerges around the X-point, while the Γ-pocket remains ungapped.
Magnetic field perpendicular to the surface Next we apply a spin-polarizing field instead of a potential term to the surface. The field (B z ) is directed perpendicular to the surface and it breaks the time reversal symmetry of the system. The perturbation is hence expected to destroy the topological protection of the surface states. We target the same two different pairs of states as for the potential term. That is, we apply the field to the subsurface Sm 4f J = 5/2, j z = ±5/2 and j z = ±1/2 spin-orbitals. As shown in Fig. S.8, the perturbing fields indeed allow the surface states around the Fermi level to hybridize, including the otherwise protected topological surface states.
Let us start with the field applied to the j z = ±5/2 spin-orbitals in Fig. S.8 (left). The large coupling be-tween the field and the jz = ±5/2 component of the surface states makes the surface bands very susceptible to the perturbing (time reversal symmetry breaking) field. The topological surface states start to hybridize and an indirect band gap opens. Quantitatively, the spectrum is altered more at the Γ point than at the X-point as here the j z = ±5/2 character dominate. But also at the X-point a small gap opens.
It is exactly vice versa if we apply the perturbing field to the subsurface Sm 4f J = 5/2, j z = ±1/2 spin-orbitals, see Fig. S.8 (right). Here, actually only the topological surface state around the X-point are gapped out, whereas the topological surface bands around Γ remains intact. While applying a field to only part of the Sm 4f states is, as a matter of course, a theoretical construct, it still demonstrates that the Γ and X pockets can clearly shift independently.