Z2 characterization for three-dimensional multiband Hubbard models
Bernhard Irsigler, Jun-Hui Zheng, Fabian Grusdt, Walter Hofstetter
ZZ characterization for three-dimensional multiband Hubbard models Bernhard Irsigler, Jun-Hui Zheng,
1, 2
Fabian Grusdt,
3, 4 and Walter Hofstetter Institut f¨ur Theoretische Physik, Goethe-Universit¨at Frankfurt am Main, Germany Center for Quantum Spintronics, Department of Physics,Norwegian University of Science and Technology, NO-7491 Trondheim, Norway Munich Center for Quantum Science and Technology (MCQST), Schellingstraße 4, 80799 M¨unchen, Germany Fakult¨at f¨ur Physik, Ludwig-Maximilians-Universit¨at, 80799 M¨unchen, Germany
We introduce three numerical methods for characterizing the topological phases of three-dimensionalmultiband Hubbard models based on twisted boundary conditions, Wilson loops, as well as the localtopological marker. We focus on the half-filled, three-dimensional time-reversal-invariant Hofstadtermodel with finite spin-orbit coupling. Besides the weak and strong topological insulator phases wefind a nodal line semimetal in the parameter regime between the two three-dimensional topologicalinsulator phases. Using dynamical mean-field theory combined with the topological Hamiltonianapproach we find stabilization of these three-dimensional topological states due to the Hubbardinteraction. We study surface states which exhibit an asymmetry between left and right surfaceoriginating from the broken parity symmetry of the system. Our results set the stage for furtherresearch on inhomogeneous three-dimensional topological systems, proximity effects, topologicalMott insulators, non-trivially linked nodal line semimetals and circuit-based quantum simulators.
I. INTRODUCTION
Three-dimensional (3d) topological states surpass theirtwo-dimensional (2d) counterparts in terms of complexityand richness. A 2d quantum spin Hall (QSH) state, e.g.,is characterized by a single Z number ν , i.e., the systemis either in a topologically trivial state ν = 0 or in thenon-trivial QSH state ν = 1. The 3d analogue, however,is characterized by four Z numbers ( ν ; ν , ν , ν ) leadingto a total of 16 topologically distinct states [1–3]. Thestraightforward way of picturing this 3d generalization isby stacking many 2d QSH layers. If the coupling betweenthese layers is weak one finds a weak topological insulator(WTI), e.g., ( ν ; ν , ν , ν ) = (0; 0 , , ν , ν , ν ) axis. The strong topological insulator (STI),on the other hand, emerges if ν = 1 and features helicalsurface states in any direction. There is no spin-conservedbackscattering of the surface states due to spin-momentumlocking which is protected by the time-reversal invariance(TRI).A further difference to the 2d case is that in 3d alsogapless topological states can emerge, such as Dirac andWeyl semimetals [6–8]. If nonsymmorphic symmetries arepresent, however, a symmetry-protected Dirac semimetalis predicted in 2d [9]. Another prominent example in 3dare nodal-line semimetals (NLSM) which exhibit a bulkband touching along a closed line embedded in the 3dBrillouin zone (BZ). These lines are not accidental but aretopologically as well as symmetry protected and cannotsimply gap out. A particle moving on a path in the 3d BZlinking the nodal line picks up a nontrivial Berry phase[10]. Here, the nodal line acts as singularity around whichthe Berry phase is acquired. The only way to open a gap isto shrink the nodal line to a point which can then gap out.On the other hand, integrating the Berry curvature on a2d manifold enclosing the complete nodal line can yield anonzero Chern number. This corresponds to a topological FIG. 1. Schematic of Hamiltonian (1). The parameters ofthe noninteracting system are described in the text. Hubbardinteractions of strength U are treated using dynamical mean-field theory leading to local selfenergies Σ i in the unit cell. charge similar to Weyl points in Weyl semimetals [8, 11,12]. If the nodal line carries a topological charge, itcannot gap out by simply shrinking to a point but has torecombine with another nodal line carrying the oppositetopological charge [10]. Even more complex physics occursif one combines many nodal lines which are topologicallynon-trivially linked [13–17].In contrast to real materials, cold atomic gases allowto experimentally rebuild model Hamiltonians such as- in the context of topological states - the celebratedHofstadter [18, 19] and Haldane [20, 21] 2d models. 3dtopological states, however, are still on their way to beexperimentally accessible. In theory there are differentgeneralizations of the Hofstadter model [22–24]. Here,we study the TRI Hofstadter model [25] generalized to3d [26] with Hubbard interactions between two fermionicspin components. We find besides WTI and STI a NLSMin the phase diagram. In order to characterize thesetopologically nontrivial phases, we generalize three topo-logical invariants to 3d, TRI, and interacting systems.Calculating surface states confirms the bulk-boundarycorrespondence.The structure of the manuscript is as follows: in Sec. II a r X i v : . [ c ond - m a t . qu a n t - g a s ] J un we introduce the 3d TRI Hofstadter-Hubbard model, inSec. III we discuss all quantum phases of the noninteract-ing system. In Secs. IV and V we put increased emphasison the Wilson loop method and the local Z marker,respectively. We continue by discussing the interactingphases within dynamical mean-field theory in Sec. VIas well as the corresponding surface states in Sec. VII.Section VIII concludes the manuscript. II. MODEL
The Hamiltonian of the 3d TRI Hofstadter-Hubbardmodel reads [26]ˆ H = (cid:88) j (cid:34) (cid:88) µ = x,y,z ( − t µ ) (cid:16) ˆ c † j + µ e πi θ µ ˆ c j + h.c. (cid:17) +( − x λ ˆ c † j ˆ c j + U ˆ c †↑ j ˆ c ↑ j ˆ c †↓ j ˆ c ↓ j (cid:35) , (1)where ˆ c j = (ˆ c ↑ j , ˆ c ↓ j ) T is the spin-1/2 fermionic annihi-lation spinor, j = ( x, y, z ) is a 3d lattice vector, t µ isthe hopping energy and µ the unit vector in µ directionwhere we set t x = t y = 1, θ = ( γσ x , αxσ z , αxσ y ) is a vec-tor of generalized matrix Peierls phases with γ being thespin-mixing amplitude, α the flux, and σ i the i th Paulimatrix. Moreover, λ is the staggered potential amplitudein x direction, and U is the Hubbard interaction energy.The Hamiltonian is schematically depicted in Fig. 1. III. NONINTERACTING PHASES
Before studying the interacting system, let us firstunderstand the noninteracting case, i.e., U = 0. Weshow the gap of the half-filled system (1) in Fig. 2i) as afunction of t z and λ for different values of γ and α = 1 / γ = 0 and 0.1 we find a gapped phase only for large λ > .
75 depending on the value of t z . For 0 . < γ < . λ emerge which wewill characterize by their topological invariants. The 3dtopological invariants of the present system, being a stackof coupled QSH layers, can be simplified through thefollowing invariant [26]: ν = Z (0) + Z ( π ) , (2)where Z ( k z ) denotes the 2d Z topological index at fixed k z . ν can thus assume the three values 0,1, or 2. Where0 represents a trivial band insulator (BI) (0;0,0,0) sinceboth 2d invariants are zero. If ν is 1, we find an STI(1;0,0,0). If both of them assume the value 1 Eq. (2)assumes 2 which corresponds to the WTI with invariants(0;0,0,1). Note that we find only this particular WTIdue to the chosen anisotropy t x = t y (cid:54) = t z . From Eq. (2)we understand that the 3d invariant only requires the FIG. 2. Noninteracting, U = 0, phase diagrams of the systemdescribed by Eq. (1): i) band gap ∆ and topological invariant ν defined in Eq. (2) for the 3d topological insulators obtainedby: ii) twisted boundary conditions, iii) Wilson loops, and iv)local Z marker for α = 1 / computation of 2d Z numbers [1–3]. In the following, wewill develop and apply three different methods in orderto compute the 2d Z indices for k z = 0 and π . We thusdirectly obtain the invariant (2). The first method isthe generalization of Fukui’s method [27] to TRI systemsusing twisted boundary conditions (TBC) [28–30]. Wefirst Fourier transform Eq. (1) for the z direction. Forthe x and y direction we apply spin-dependent TBC, i.e.,ˆ c x + N x ,y,k z = ˆ c x,y,k z e iϑ x and ˆ c x,y + N y ,k z = ˆ c x,y,k z e iϑ y σ z where N x × N y is the size of the 2d system. Note that thespin dependence σ z appears only once, however, in whichdirection is a freedom of the gauge. After introducingthe twist angles ϑ x and ϑ y , Fukui’s method is applied in( ϑ x , ϑ y ) space. This yields the Z invariant Z ( k z ) withparameter k z . We show ν in Fig. 2ii) obtained by theTBC method if the gap i) is finite.We find an STI phase, shown in red, as well as a WTI,shown in yellow. We also observe in Fig. 2e) that formaximal spin mixing γ = 0 .
25 there are gapless transi-tion lines between the topological insulator phases. For γ = 0 .
22 as shown in Fig. 2d) these transition lines extendto gapless regions which we discuss further below.
IV. WILSON LOOP
The second method to compute ν is the Wilson looptechnique [31, 32] which is an extension of the Zak phaseto multi-band systems. We present the Wilson loop tech-nique in the 3d case and provide details for the numericalcomputation in the following.We Fourier transform the Hamiltonian, defined inEq. (1), for all three spatial dimensions. For α = 1 / k -dependent Hamiltonian matrix has Hilbertspace dimension 12 where the spin as well as the position x within the unit cell are treated as internal degrees offreedom: H ( k ) = O (1) T T † e ik x T † O (2) TT † O (3) TT † O (4) TT † O (5) TT e − ik x T † O (6) T = te πiγσ x (3) O ( x ) = − t cos( k y ) cos(2 παx ) − t sin( k y ) sin(2 παx ) σ z − t z cos( k z ) cos(2 παx ) − t z sin( k z ) sin(2 παx ) σ y + λ ( − x . We define the time-reversal-invariant, gauge-independentmultiband formulation [31] of the discretized Wilson loop D ( C k ) = (cid:89) k j ∈ C k F j , with F mnj = (cid:104) u m ( k j ) | u n ( k j +1 ) (cid:105) , (4)with | u n ( k j ) (cid:105) being the cell-periodic part of Bloch stateof the n th band and k j are discretized values of the closedcontour C k in the BZ. If we set k z = 0 , π and choose C k to go along k x , we find Eq. (4) to be a parametric functionof k y only. The eigenvalues of D ( k y ) are λ m ( k y ). Theirphases θ m ( k y ) = Im log λ m ( k y ) will perform trajectorieson a cylinder which we define through ( k y , θ m ) ∈ [0 , π ] × [0 , π ]. Here, θ m is the periodic part of the cylinder. Atthe ends of the cylinder, i.e., k y = 0 and π , the θ m willbe degenerate in pairs due to time-reversal symmetry.By tuning k y from 0 to π these pairs will split and the θ m may wind around the cylinder. At k y = π the θ m reconnect again in pairs. This integer valued windingnumber around the cylinder is directly connected to thetime-reversal polarization [33] and corresponds to the Z number.Numerically, we find this winding number by dividingthe cylinder into three regions: I where 0 < θ m < π/
3, IIwhere 2 π/ < θ m < π/
3, and III where 4 π/ < θ m < π .The winding is depicted in Fig. 3 for a trivial (blue) anda nontrivial instance (orange). We sample a sufficientset of values of k y and count the number n i of θ m valuesbeing in the region i , with i =I,II,III. This yields the data( n , n , n ) as a function of k y . We then compute thechange ∆ n i of n i with respect to ∆ k y . From this data,we only keep the ones where (∆ n , ∆ n , ∆ n ) follow somepermutation of − , ,
1. Data, where (∆ n , ∆ n , ∆ n ) allare zero, do not carry information and data where a 2appears could be removed by increased sampling of k y FIG. 3. Example for the numerical calculation of the Wilsonloops. The parameters are set to γ = 0 . t z = 0 . λ = 0 .
75. The blue data points correspond to k z = 0 and theorange data points to k z = π . and can thus be safely omitted. Finally, to each of theremaining data points a chirality can be assigned by meansof the Levi-Civita tensor. Summing these chiralities yieldsa nontrivial Z number for odd and a trivial one for evenvalues of the sum of the chiralities. The results are shownin Fig. 2iii) and they agree exactly with the method usingTBC in Fig. 2ii). V. LOCAL Z MARKER
We now turn to the generalization of the local Chernmarker [34] to the TRI case, first in 2d. A Z general-ization to Kitaev’s real-space formulation of the Chernnumber [35] has recently appeared in Ref. [36]. Here, weintroduce the spin-projected version of the local Chernmarker C µν ( x, y ) = (cid:104) x, y |P µ ˆ P ˆ x ˆ P ˆ y ˆ P P ν | x, y (cid:105) , (5)where P µ is the projector onto the states of band µ =I,IIin the { I,II } eigenbasis of time-reversed partners [33], ˆ P isthe projector onto the occupied eigenstates of the Hamil-tonian, and | x, y (cid:105) is the eigenstate of the position operatorin 2d. The eigenvalues of the 2 × C µν ( x, y ) cor-respond to time-reversed partners similar to the partialpolarizations of the time-reversed partners in Ref. [33],however, now defined in real space. The first eigenvaluethus resembles exactly the 2d local Z marker. Sincethe eigenvalues are independent of the basis in which C µν ( x, y ) is represented, we can also use the spin basis {↑ , ↓} such that we do not have to find the { I,II } basis.By Fourier transforming only the z coordinate of Eq. (1)and fixing the value of k z = 0 , π we can generalize the 2dlocal Z marker to a 3d local Z marker. FIG. 4. Interacting phase diagrams obtained from DMFTand the topological Hamiltonian Eq. (6). a) shows the bandgap of the topological Hamiltonian ∆ and b) the topologicalinvariant ν , Eq. (2), as functions of the interaction strength U and the staggered potential λ for γ = 0 .
25 and t z = 0 . ν as functions of U and t z for γ = 0 .
22 and λ = 0 .
75. The orange lines correspond to theorange lines in Fig. 2e)i) and ii). The gray regions denoteDMFT solutions which break the lattice symmetry. The whitesymbols correspond to the parameter sets used in Fig. 6.
The bulk value is presented in Fig. 2iv) showing ap-proximately the same behaviour as the aforementionedmethods in Fig. 2ii) and iii) computed on a 30 ×
30 lat-tice. However, the local Z marker suffers from finite sizeeffects when the gap is small [37]. This can be observed,e.g., in Fig. 2iv)e) for large t z and λ where the local Z marker is not quantized due to the finite system. If thegap is sufficiently large, however, the local Z marker iswell quantized. VI. INTERACTING PHASES
We study interaction effects by applying dynamicalmean-field theory (DMFT) [38] which neglects nonlocalfluctuations but covers all local fluctuations. Since theunit cell of the system (1) contains six lattice sites if thereis no spontaneous symmetry breaking, we make use ofthe real-space version of DMFT [39–41]. Here, the many-body problem of the full lattice with N sites is mappedonto N single-site quantum impurity problems, whereeach impurity problem interacts with a self-consistent,noninteracting bath. This approach nonperturbativelydescribes local quantum dynamics, in contrast to staticmean-field theory. After solving the single-impurity prob-lem for each site, for which we use exact diagonalizationwith four bath sites, the selfenergy Σ σσ (cid:48) j ( ω ) for each lat-tice site j and frequency ω is obtained. Here, σ, σ (cid:48) arespin degrees of freedom. Using the Dyson equation, theseare used to construct a new lattice Green’s function andthis procedure is repeated until self-consistency. In 2dDMFT has provided a successful description of topolog-ical systems for many aspects [30, 42–48]. Ref. [49] hasshown that nonlocal contributions are small already in2d. We therefore expect even more accurate results ofDMFT in 3d.To calculate topological invariants for the interactingsystem we follow the topological Hamiltonian approach [50]. The idea here is, that if the Green’s function can besmoothly deformed to a noninteracting Green’s function,i.e., no poles or zeros occur, the topological properties donot change. This holds since topological phase transitionscome along with a divergence or a zero [51, 52] of theGreen’s function. In this way, one can construct an ef-fective, noninteracting Hamiltonian H T = − G − ( ω = 0)which is used to compute topological invariants. In com-bination with the local selfenergy Σ σσ (cid:48) j ( ω ) from DMFTits matrix form reads:[ H T ] σσ (cid:48) jj (cid:48) = [ H ] σσ (cid:48) jj (cid:48) + Σ σσ (cid:48) j ( ω = 0) δ jj (cid:48) , (6)where H denotes the noninteracting part of the Hamil-tonian. We show the gap ∆ as well as the topologicalinvariant (2) of the topological Hamiltonian (6) in Figs. 4a)and b), respectively, as functions of λ and U for γ = 0 . t z = 0 .
2. The green lines for U = 0 correspond to thegreen lines in Figs. 2e)i) and ii), respectively. The greyregions correspond to DMFT results where the latticesymmetry is spontaneously broken.We first focus on the symmetric phases. We observethat the gap closing lines in Fig. 4a) coincide with thetopological phase transition lines in b) as expected. Fur-thermore, we find stabilization of the STI and the WTIphases against λ through U . For small U , Hubbard in-teractions effectively renormalize λ , which extends thetopological phases in the phase diagram. This is the 3danalogue of the interaction-induced topological phase tran-sition in 2d [30, 53] and can be understood through thecompetition between staggered potential and interactions.In Figs. 4c) and d) we present ∆ and ν as functionsof t z and U . The orange lines for U = 0 correspond tothe orange lines in Figs. 2d)i) and ii), respectively, for γ = 0 .
22 and λ = 0 .
75. As in the previous phase diagram,we observe again stabilization of the topological phasesthrough interactions.We compare our results to a DMFT study of a four-band model including a Hund’s coupling term [54]. Wefind qualitative agreement between the phase diagram inFig. 4d) and the one in Ref. [54, Fig. 3] even though thelatter corresponds to a finite Hund’s coupling which wedo not include here. This is because the Hund’s couplingeffectively reduces the interorbital interactions and thusmakes the Hubbard term the dominant interaction term.In contrast to Ref. [54], we do not find the (1;1,1,1) phase.It is a priori not clear what would be the unit cell of apossible spontaneous-symmetry-broken phase as a resultof the nontrivial exchange couplings between neighboringspins due to the Peierls phases in Eq. (1). The effectivespin Hamiltonian [42] for the 3d system reads:ˆ H spin = (cid:88) j (cid:88) µ,ν,ρ cyclic t µ U (cid:110) ˆ S µ j ˆ S µ j + µ + cos(4 πθ µ ) (cid:104) ˆ S ν j ˆ S ν j + µ + ˆ S ρ j ˆ S ρ j + µ (cid:105) + sin(4 πθ µ ) (cid:104) ˆ S ν j ˆ S ρ j + µ − ˆ S ρ j ˆ S ν j + µ (cid:105)(cid:111) , (7) FIG. 5. Classical Monte Carlo results for the groundstateenergy of Eq. (7) for a) γ = 0 .
25 and t z = 0 . γ = 0 .
22 and t z = 1 as function of the size of the unit cell N x × N y × N z . b) and d) show the spin state of the unit cellwith the smallest energy marked by a blue circle in a) and c),respectively. where we defined the spin operator ˆ S µ j = ˆ c † j σ µ ˆ c j and θ = ( γ, αx, αx ). In the spin population balanced, 2d case[43, 47] one can argue that the spins will always orderantiferromagnetically in y direction. However, in the 3dcase we cannot make this argument and the unit cellmight in fact be very large. Results of a classical MonteCarlo study to find the classical ground state E of Eq. (7)are shown in Fig. 5 for unit cells up to N x × N y × N z = 6 lattice sites. Examplarily, for γ = 0 .
25 and t z = 0 . × × x and anti-ferromagnetic ordering in y and z direction similar to the 2d collinear order in Ref. [43,Fig. 6]. For γ = 0 .
22 and t z = 1 the unit cell with smallestenergy is found with dimensions 4 × × γ and t z and can lead to complexmagnetic orders. VII. SURFACE STATES
We now study the surface states of the present sys-tem. To this end we put the system on a 3d cylindergeometry, i.e., k y , k z are good quantum numbers but in x direction we now apply open boundary conditions. Forthis geometry, we can define the single-particle Green’s function: G σσ (cid:48) xx (cid:48) ( ω, k y , k z ) = (cid:104) { ω − Σ( ω ) − H ( k y , k z ) } − (cid:105) σσ (cid:48) xx (cid:48) , (8)where x, x (cid:48) are the spatial degrees of freedom in x direc-tion. The spectral density of a spatial region X is definedas ρ X ( ω, k x , k y ) = − π (cid:88) σ,x ∈ X Im G σσxx ( ω, k y , k z ) . (9)We show the surface states of a system with 60 sites in x direction by plotting ρ X ( ω = 0 , k y , k z ) in Fig. 6 forthe left surface L where 1 ≤ x ≤
3, the bulk B where25 ≤ x ≤
36, and the right surface R where 58 ≤ x ≤ k y = 0and k z = 0. The red dots denote TRI momenta. Theparameters are chosen according to the white symbolsin Fig. 4. For U = 1, Fig. 6a) shows the Fermi surfaceenclosing only one TRI momentum which corresponds tothe surface state of an STI. Figure 6b) shows the Fermisurface crossing the BZ almost parallel to the k z axis.Thus it encloses two TRI momenta which corresponds tothe surface state of a WTI. We do not show the resultsfor the bulk because it is gapped.Let us now turn to the case of the NLSM. Figures 6c)to e) show the surface states as well as the x projectionof the bulk state for different U . First, we notice that thetwo surfaces are asymmetric. This arises from the brokenparity symmetry since upon the transformation λ →− λ the spectral density behavior of the two surfaces isexchanged. The right surface shows a state correspondingto a WTI, whereas the left surface rather shows onecorresponding to the STI. However, the surface state ofthe left surface does not fully enclose the TRI momentumbut rather stops within the BZ, reminiscent of a Fermi arc.The missing part is recovered in the bulk and identifiesas a nodal line. It is shown for the full 3d BZ in Fig. 6f).In order to justify that the nodal line is not an acci-dental band touching but is topologically protected wecompute the Berry phase on a closed path in the 3dBZ linking the nodal line. To this end, we make use ofthe multiband formulation of Ref. [31] to find the Berryphase ImTr log (cid:81) j F j , where F mnj = (cid:104) u m ( k j ) | u n ( k j +1 ) (cid:105) and | u n ( k j ) (cid:105) is the cell-periodic part of Bloch state ofthe n th band. If the path k j is (is not) linked with thenodal line the Berry phase yields π (0). This shows thatthe nodal line is topologically protected and thus cannotgap out. On the other hand, when computing the Berrycurvature on a 2d box surface enclosing the nodal line wefind a vanishing Chern number. Thus the nodal line doesnot carry a topological charge.Refs. [55–57] studied the interacting NLSM with renor-malization group calculation, static mean-field, and clus-ter perturbation theory, respectively. Since the density of FIG. 6. Surface states of different 3d, interacting, topologicallynontrivial phases: a) strong topological insulator, b) weaktopological insulator, and c)-e) nodal-line semimetals. Thewhite symbols correspond to the parameter sets marked bythe symbols in Fig. 4. See definitions of L, B, and R belowEq. (9). The bulk nodal line in the full 3d BZ is shown inf) in blue corresponding to the projected nodal line in c)B;projections onto the k i - k j planes for i, j = x, y, z are shownin red. states of a NLSM vanishes at the band touching the nodalline is robust against interactions and will only gap outfor strong interactions. On the other hand, interactionscan change the size and the shape of the nodal line.The surface states of NLSMs have attracted a lot ofattention. If they are drumhead states, they constituteflat bands with a diverging density of states which islocalized at the surfaces of the system. Ref. [58] havestudied emergent surface antiferromagnetic order in thiscontext already for small critical interaction strength,which is, however, increased by spin-orbit coupling. We do not observe any symmetry breaking in our parame-ter range and do not find this surface magnetism. Weattribute this to the strong spin-orbit coupling in oursystem which curves the flat bands and thus decreasesthe surface density of states. A recent study found aninversion of the Berry curvature of one K point driven byspin-orbit coupling and two-particle interactions in theHaldane model, which leads to a surface Chern insulator[59]. Bulk antiferromagnetism in NLSMs was investigatedin Ref. [60].From the experimental point of view, the TRI Hofs-tadter Hamiltonian has been realized in 2d using laser-assisted tunneling [18]. Theoretically, this approach hasbeen generalized to 3d [15, 24, 61]. A generic way to im-plement spin-orbit coupling proposed in Ref. [62] might begeneralizable to three dimensions. For detection, Bloch-Zener-St¨uckelberg interferometry [61], anomalous velocitymeasurement, or state tomography [24] for bulk statesas well as Bragg spectroscopy for surface states [15] havebeen proposed. Also a 3d version of a topological interface[46] could be used for the detection of the surface states.The here introduced local Z marker could be used todistinguish the topological phase at the interface. Veryrecently, a NLSM has been realized in a fermionic coldatom experiment with Yb atoms by mapping the k z component to a Zeeman field and reading out 2d layersfor each value of k z [63]. VIII. CONCLUSION
We develop three numerical techniques for the charac-terization of three-dimensional topological states of mat-ter. Based on twisted boundary conditions, Wilson loops,and the local topological marker, these techniques canbe used to compute weak and strong topological indiceseven in interacting systems. We apply these to the three-dimensional time-reversal-invariant Hofstadter-Hubbardmodel and find a topological nodal-line semimetal betweenphases of weak and strong topological insulators. Usingdynamical mean-field theory we observe stabilization ofthe three-dimensional topological states through Hub-bard interactions. The numerical methods presented hereenable the study of interacting, three-dimensional topo-logical matter in inhomogeneous systems, which will be ofgreat interest for cold atomic implementations. Moreover,we think that our results could contribute to benchmarkcircuit-based quantum simulators where is has been pos-sible to engineer artificial gauge fields as well as stronginteractions [64–66].
ACKNOWLEDGMENTS
The authors acknowledge useful discussions withMohsen Hafez-Torbati. This work was supported bythe Deutsche Forschungsgemeinschaft (DFG, German Re-search Foundation) under Project No. 277974659 via Re-search Unit FOR 2414 and Germany’s Excellence Strategy- EXC2111 - 390814868. This work was also supported by the DFG via the high performance computing centerLOEWE-CSC. [1] L. Fu, C. L. Kane, and E. J. Mele, Phys. Rev. Lett. ,106803 (2007).[2] J. E. Moore and L. Balents, Phys. Rev. B , 121306(R)(2007).[3] R. Roy, Phys. Rev. B , 195322 (2009).[4] Z. Ringel, Y. E. Kraus, and A. Stern, Phys. Rev. B ,045102 (2012).[5] B. Sbierski, M. Schneider, and P. W. Brouwer, Phys. Rev.B , 161105(R) (2016).[6] S. Murakami and S.-i. Kuga, Phys. Rev. B , 165313(2008).[7] S. M. Young, S. Zaheer, J. C. Y. Teo, C. L. Kane, E. J.Mele, and A. M. Rappe, Phys. Rev. Lett. , 140405(2012).[8] N. P. Armitage, E. J. Mele, and A. Vishwanath, Rev.Mod. Phys. , 015001 (2018).[9] S. M. Young and C. L. Kane, Phys. Rev. Lett. , 126803(2015).[10] C. Fang, H. Weng, X. Dai, and Z. Fang, Chinese Phys.B , 117106 (2016).[11] A. A. Burkov, M. D. Hook, and L. Balents, Phys. Rev.B , 235126 (2011).[12] M. Hirayama, R. Okugawa, and S. Murakami, J. Phys.Soc. Jpn. , 041002 (2018).[13] T. Bzdusek, Q. Wu, A. R¨uegg, M. Sigrist, and A. A.Soluyanov, Nature , 75 EP (2016).[14] G. Chang, S.-Y. Xu, X. Zhou, S.-M. Huang, B. Singh,B. Wang, I. Belopolski, J. Yin, S. Zhang, A. Bansil, H. Lin,and M. Z. Hasan, Phys. Rev. Lett. , 156401 (2017).[15] W. Chen, H.-Z. Lu, and J.-M. Hou, Phys. Rev. B ,041102(R) (2017).[16] Z. Yan, R. Bi, H. Shen, L. Lu, S.-C. Zhang, and Z. Wang,Phys. Rev. B , 041103(R) (2017).[17] L. Li, C. H. Lee, and J. Gong, Phys. Rev. Lett. ,036401 (2018).[18] M. Aidelsburger, M. Atala, M. Lohse, J. T. Barreiro,B. Paredes, and I. Bloch, Phys. Rev. Lett. , 185301(2013).[19] H. Miyake, G. A. Siviloglou, C. J. Kennedy, W. C. Burton,and W. Ketterle, Phys. Rev. Lett. , 185302 (2013).[20] G. Jotzu, M. Messer, R. Desbuquois, M. Lebrat,T. Uehlinger, D. Greif, and T. Esslinger, Nature ,237 (2014).[21] N. Fl¨aschner, B. S. Rem, M. Tarnowski, D. Vogel, D.-S.L¨uhmann, K. Sengstock, and C. Weitenberg, Science , 1091 (2016).[22] T. Kimura, Prog. Theor. Exp. Phys. , 103B05 (2014).[23] Y. Li, Phys. Rev. B , 195133 (2015).[24] D.-W. Zhang, R.-B. Liu, and S.-L. Zhu, Phys. Rev. A , 043619 (2017).[25] N. Goldman, I. Satija, P. Nikolic, A. Bermudez, M. A.Martin-Delgado, M. Lewenstein, and I. B. Spielman,Phys. Rev. Lett. , 255302 (2010).[26] M. S. Scheurer, S. Rachel, and P. P. Orth, Sci. Rep. ,8386 (2015).[27] T. Fukui, Y. Hatsugai, and H. Suzuki, J. Phys. Soc. Jpn. , 1674 (2005). [28] D. N. Sheng, Z. Y. Weng, L. Sheng, and F. D. M. Haldane,Phys. Rev. Lett. , 036808 (2006).[29] T. Fukui and Y. Hatsugai, Phys. Rev. B , 121403(R)(2007).[30] P. Kumar, T. Mertz, and W. Hofstetter, Phys. Rev. B , 115161 (2016).[31] R. Yu, X. L. Qi, A. Bernevig, Z. Fang, and X. Dai, Phys.Rev. B , 075119 (2011).[32] F. Grusdt, D. Abanin, and E. Demler, Phys. Rev. A ,043621 (2014).[33] L. Fu and C. L. Kane, Phys. Rev. B , 195312 (2006).[34] R. Bianco and R. Resta, Phys. Rev. B , 241106(R)(2011).[35] A. Kitaev, Ann. Phys. , 2 (2006).[36] Z. Li and R. S. Mong, arXiv:1905.12649 (2019).[37] B. Irsigler, J.-H. Zheng, and W. Hofstetter, Phys. Rev.A , 023610 (2019).[38] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg,Rev. Mod. Phys. , 13 (1996).[39] S. Okamoto and A. J. Millis, Phys. Rev. B , 241104(R)(2004).[40] R. W. Helmes, T. A. Costi, and A. Rosch, Phys. Rev.Lett. , 056403 (2008).[41] M. Snoek, I. Titvinidze, C. T˝oke, K. Byczuk, and W. Hof-stetter, New J. Phys. , 093008 (2008).[42] D. Cocks, P. P. Orth, S. Rachel, M. Buchhold, K. Le Hur,and W. Hofstetter, Phys. Rev. Lett. , 205303 (2012).[43] P. P. Orth, D. Cocks, S. Rachel, M. Buchhold, K. Le Hur,and W. Hofstetter, J. Phys. B , 134004 (2013).[44] T. I. Vanhala, T. Siro, L. Liang, M. Troyer, A. Harju,and P. T¨orm¨a, Phys. Rev. Lett. , 225305 (2016).[45] A. Amaricci, L. Privitera, F. Petocchi, M. Capone, G. San-giovanni, and B. Trauzettel, Phys. Rev. B , 205120(2017).[46] B. Irsigler, J.-H. Zheng, and W. Hofstetter, Phys. Rev.Lett. , 010406 (2019).[47] B. Irsigler, J.-H. Zheng, M. Hafez-Torbati, and W. Hof-stetter, Phys. Rev. A , 043628 (2019).[48] U. Gebert, B. Irsigler, and W. Hofstetter,arXiv:1906.11164 (2019).[49] T. Mertz, K. Zantout, and R. Valent´ı, Phys. Rev. B ,125111 (2019).[50] Z. Wang and S.-C. Zhang, Phys. Rev. X , 031008 (2012).[51] V. Gurarie, Phys. Rev. B , 085426 (2011).[52] J.-H. Zheng and W. Hofstetter, Phys Rev. B , 195434(2018).[53] J.-H. Zheng, B. Irsigler, L. Jiang, C. Weitenberg, andW. Hofstetter, arXiv:1812.01991 (2018).[54] A. Amaricci, J. C. Budich, M. Capone, B. Trauzettel,and G. Sangiovanni, Phys. Rev. B , 235112 (2016).[55] S. Sur and R. Nandkishore, New J. Phys. , 115006(2016).[56] B. Roy, Phys. Rev. B , 041113(R) (2017).[57] J. Kang, J. Zou, K. Li, S.-L. Yu, and L.-B. Shao, Sci.Rep. , 2824 (2019).[58] J. Liu and L. Balents, Phys. Rev. B , 075426 (2017). [59] W. Chen and J. L. Lado, Phys. Rev. Lett. , 016803(2019).[60] J. Wang, Phys. Rev. B , 081107(R) (2017).[61] D.-W. Zhang, Y. X. Zhao, R.-B. Liu, Z.-Y. Xue, S.-L.Zhu, and Z. D. Wang, Phys. Rev. A , 043617 (2016).[62] F. Grusdt, T. Li, I. Bloch, and E. Demler, Phys. Rev. A , 063617 (2017).[63] B. Song, C. He, S. Niu, L. Zhang, Z. Ren, X.-J. Liu, andG.-B. Jo, Nat. Phys. , 911 (2019). [64] J. Koch, A. A. Houck, K. L. Hur, and S. M. Girvin, Phys.Rev. A , 043811 (2010).[65] P. Roushan, C. Neill, A. Megrant, Y. Chen, R. Bab-bush, R. Barends, B. Campbell, Z. Chen, B. Chiaro,A. Dunsworth, A. Fowler, E. Jeffrey, J. Kelly, E. Lucero,J. Mutus, P. J. J. O’Malley, M. Neeley, C. Quintana,D. Sank, A. Vainsencher, J. Wenner, T. White, E. Kapit,H. Neven, and J. Martinis, Nat. Phys. , 146 (2017).[66] C. Owens, A. LaChapelle, B. Saxberg, B. M. Anderson,R. Ma, J. Simon, and D. I. Schuster, Phys. Rev. A97