A search for correlation-induced adiabatic paths between distinct topological insulators

Correlations in topological states of matter provide a rich phenomenology, including a reduction in the topological classification of the interacting system compared to its non-interacting counterpart. This happens when two phases that are topologically distinct on the non-interacting level become adiabatically connected once interactions are included. We use a quantum Monte Carlo method to study such a reduction. We consider a 2D charge-conserving analog of the Levin-Gu superconductor whose classification is reduced from $\mathbb{Z}$ to $\mathbb{Z}_4$. We may expect any symmetry-preserving interaction that leads to a symmetric gapped ground state at strong coupling, and consequently a gapped symmetric surface, to be sufficient for such reduction. Here, we provide a counter example by considering an interaction which (i) leads to a symmetric gapped ground state at sufficient strength and (ii) does not allow for any adiabatic path connecting the trivial phase to the topological phase with $w=4$. The latter is established by numerically mapping the phase diagram as a function of the interaction strength and a parameter tuning the topological invariant. Instead of the adiabatic connection, the system exhibits an extended region of spontaneous symmetry breaking separating the topological sectors. Frustration reduces the size of this long-range ordered region until it gives way to a first order phase transition. Within the investigated range of parameters, there is no adiabatic path deforming the formerly distinct free fermion states into each other. We conclude that an interaction which trivializes the surface of a gapped topological phase is necessary but not sufficient to establish an adiabatic path within the reduced classification. In other words, the class of interactions which trivializes the surface is different from the class which establishes an adiabatic connection in the bulk.


I. INTRODUCTION
Recent years have witnessed an intense research effort to understand topological phases of matter [1][2][3][4] . Symmetryprotected topological (SPT) phases are described by equivalence classes of phases under symmetric adiabatic deformation. This means that two SPTs belonging to the same phase can be deformed to each other without closing the bulk gap or breaking the protecting symmetries, whereas two distinct SPT phases cannot. SPTs protected by internal symmetries such as time-reversal or particle-hole symmetry have been extensively studied for free-fermion systems 5,6 . These are all characterized by the existence of gapless anomalous surface states whose existence is a direct consequence of the bulk topology, a phenomenon known as bulk-boundary correspondence.
The inclusion of interactions can modify the topological character of free-fermion SPTs in at least three different ways: (i) The spontaneous breaking of protecting symmetry can lead to the disappearance of surface states and consequently alter the topological classification. (ii) Correlations can induce topological order that is characterized by long-range entanglement as in fractional quantum Hall states 7,8 , fractional topological insulators 9 , or quantum spin liquids [10][11][12][13] . These states do not have a non-interacting analog. (iii) The free-fermion topological classification may be reduced due to the existence of symmetric adiabatic paths in the space of interacting Hamiltonians connecting states that are disconnected at the non-interacting level [14][15][16][17][18][19][20][21][22][23][24][25] .
The first example of a single-particle topology reduction was considered by Fidkowski and Kitaev in Ref. 14. In this work, the authors study a spinless superconductor in one dimension (Kitaev chain), with spinless time-reversal symmetry T 2 = +1, representing class BDI. They have constructed an explicit interaction which preserves the symmetry and has a unique and symmetric ground state, that yet gaps out 8 topological Majorana boundary modes and allows an adiabatic connection between bulk states whose winding number differ by 8. This implies a reduction of the noninteraction classification from Z to Z 8 . Later on, this result has been generalized to different symmetry classes and higher dimensions [15][16][17][18][19][20][21][22][23][24][25] . Whereas the early work by Fidkowski and Kitaev 14 used the simple properties of the 1D model to show explicitly that two phases differing by winding of 8 can be adiabatically connected, most of the investigations in higher dimensional systems relied on bulk-boundary correspondence to argue that the existence of an interaction which symmetrically gaps out the topological boundary modes is sufficient to establish the collapse of the non-interacting classification. These approaches employed several arguments such as studying the possibility of gapping surface states [15][16][17] or 0D defects that follow from dimensional reduction 18,19 as well as investigating the signatures of these boundary states in the entanglement spectrum 20,21 . Other bulk-based approches include studying the braiding statistics arising from gauging the symmetry 22 or group cohomology 23 .
In spite of this body of work, which provided a comprehensive answer to the question of the general interactionreduced classification for topological phases protected by internal symmetries, only a few works 26 addressed the fate of a given topological phase in the presence of a specific symmetric interaction. In particular, we would like to investigate the question whether an interaction which symmetrically gaps arXiv:1912.07614v1 [cond-mat.str-el] 16 Dec 2019 out surface states at sufficient strength also enables an adiabatic symmetric path connecting two distinct non-interacting SPTs. In this work, we show that this is not generally true. As we will illustrate in detail, the caveat is that the interaction results in transitions into symmetry-broken phases along the paths, preventing adiabaticity despite the fact that at strong enough interaction strength the ground state is symmetric.
To show this, we consider a simple model of four identical layers of topological insulators protected by charge conservation and an internal Z 2 symmetry (a non-superconducting analog of the one considered by Levin and Gu 22 ) whose non-interacting Z classification is reduced to Z 4 . The noninteracting theory has an emergent SU (4) symmetry corresponding to rotations among the four different flavors and is characterized by the topological winding number w = 4. However, we consider an interaction of much lower symmetry which reduces this flavor rotation symmetry to U (1) × U (1) (while still preserving the symmetries protecting the topological classification). It is worth mentioning here that our approach differs from earlier works 21,26,27 which considered highly symmetric interactions to avoid possible symmetry breaking. Here, we consider an interaction with very low symmetry precisely to show that symmetry breaking is unavoidable along the adiabatic path although the ground state is still symmetric at large enough interaction strength.
Starting from a microscopic model, we employ the projective auxiliary-field Quantum Monte Carlo method [28][29][30][31] to study its ground state phase diagram. We find that adiabatic connection between the w = 4 and w = 0 phases of the model is not possible due to the appearance of an extended region of spontaneous symmetry breaking that additionally separates the single particle topological sectors. To overcome this problem, we add extra terms to the Hamiltonian to frustrate the long range order. At weak levels of frustration, the region of spontaneous symmetry breaking is reduced in size until it gives way to a first order phase transition at strong frustration that still blocks the adiabatic connection. As a result, we conclude that the interaction we considered, while sufficient for gapping out surface states, is insufficient for the existence of a symmetric adiabatic deformation between the trivial and nontrivial phases. This poses a counter example to the criteria derived by some of the authors in Ref. 19 which is shown to be necessary but not sufficient.
This article is organized as follows. In Sec. II, we define the two-dimensional microscopic model and discuss its symmetries.In Sec. III, we analyze possible mean-field scenarios with spontaneous symmetry breaking and identify the most dangerous channel. Additionally, we present the analytic solution in the limit of infinite interaction strength and find a unique, symmetric and gapped ground state. In Sec. IV, we briefly discuss the projective auxiliary-field Quantum Monte Carlo (QMC) method. We present the numerical extracted phase diagram in Sec. V that exhibits a region of spontaneous symmetry breaking that gives way to a first order phase transition at strong frustration. In Sec. VI, we conclude with a discussion of the implications of the phase diagram and suggest future avenues.

II. MODEL & SYMMETRIES
Here, we design a two-dimensional microscopic model that obeys an anti-unitary time-reversal symmetry (TRS) and a unitary Z 2 symmetry 18,19 . The latter can be implemented, e.g., as the conservation of the z component of the spin S z modulo 2. The topology of these free fermion model is given by a Z valued winding number w that is related to the number of helical Dirac cones at the edge of the sample. In the presence of correlations, the classification is expected to be reduced from Z to Z 4 . We define a specific interaction term that should allow for adiabatic deformations of free fermion states whose winding number differ by (multiples of) ∆w = ±4.
We begin by introducing the free fermion part of the model H 0 = k Ψ † k H(k)Ψ k which represents two copies of a Quantum-Hall system with opposite Chern number representing a topological insulator [32][33][34][35] which preserves the spin projection σ z . We have where Ψ † k is the creation operator of a four-component spinor with momentum k. The above Hamiltonian is block diagonal and we denote the sub-blocks by H ± (k). They represent a Chern insulator with a ±1 Chern number 36 . The energy scale is set by t which will be used as the unit of energy (t = 1) throughout the rest of this manuscript. Eq. (1), corresponds to a gapped Dirac cone at each of the four time-reversal invariant momenta, with the gap dictated by m(k). For the choice of parameters λ = 0, λ = −2 or λ = −4 at least one of the Dirac cones remains gapless (compare with Fig. 1) whereas any other value describes an insulating state. These points correspond to topological phase transitions. The last term in Eq. (1a), ∆ = 0.25 is used to break the particle-hole symmetry within H ± (k). We note that the dimensional reduction arguments used in Ref. 19 relied on the existence of chiral symmetry in the continuum model. The symmetry is explicitly broken here following the aforementioned scheme of having as little symmetries as possible. Also in real materials, the particle-hole symmetry is an approximate low energy symmetry, unless one considers superconductors.
The Hamiltonian H(k) obeys a spurious U (1) symmetry exp{iφσ z } instead of the required Z 2 Ising symmetry R = σ z . It separates the H ± (k) sectors, where ± is given by the eigenvalues of R. Additionally, there is one independent, anti-unitary time-reversal symmetry T = σ y τ y K. Here, K refers to the complex conjugation. Hence, the model is very closely related to the well know topological insulators (TIs). There, one may introduce spin-orbit coupling which breaks the spin conservation as long as the time-reversal symmetry is respected.
In order to discuss the topology of the gapped states, it is useful to first focus on H + . Defining the three-component vector d = (sin(k x ), sin(k y ), m(k)), the Chern number of withd = d/|d| 33 . As it can be seen in Fig. 1, the parameter λ tunes the system from a trivial insulator (λ > 0) through a semi-metal with a Dirac cone at k = (π, π) (λ = 0) to a Chern insulator (−2 < λ < 0) with Chern number +1. At λ = −2, the system exhibits one Dirac cone at each k = (0, π) and k = (π, 0). The winding number of the full Hamiltonian is then given as w = (C + −C − )/2 where C − = −C + due to the time-reversal symmetry connecting the two sectors. Note that small values of ∆ only modify the energy of the bands and, as long as the band gap does not close, the wave functions do not change as it represents a (momentum dependent) chemical potential. Hence the topology is insensitive to (small) ∆.
To study the topological reduction from Z to Z 4 , we introduce four copies of H 0 labelled by an orbital index o ∈ {A, B, C, D}. The dimensional reduction scheme 19 can be used to derive the form of the interactions that allow the adiabatic connection by introducing a lattice of zerodimensional defects. Such defects inherit the bulk topology in the sense that they exhibit n topologically degenerate zero modes where n matches the topological invariant of the bulk 37 . For two-dimensional models, this is done by first realizing one-dimensional edge modes at a domain wall and secondly adding an oscillating mass term along this domain wall with appropriately chosen symmetries, such that each node of the mass term localizes zero modes. This construction in turn allows to derive an explicit interaction term which gaps those defects without breaking any symmetry. Using this recipe, we design the following interaction with the projectors P α = 1 2 (1 + iαγ 3 γ 4 ) satisfying P = P 2 , and the γ matrices acting on the original Dirac components as γ 1,2 = σ 0 τ x,y , and γ 3,4,5 = σ z,y,x τ z . Note that the spurious U (1) symmetry is R = iγ 4 γ 5 and that the operators R, γ 4 and γ 5 form an SU (2) algebra. Hence, the symmetry generates continuous rotations in the γ 4 γ 5 plane. Using only γ 5 in the interaction terms breaks the U (1) symmetry down to the required discrete Z 2 symmetry that transforms γ 5 → −γ 5 . Physically, this interaction introduces correlated pair hopping of electrons between layers A → B and D → C while flipping the R charge. This term does allow, e.g., two R = + electrons being scattered into two R = − ones such that the R charge is only conserved modulo 2 which illustrates the U (1) → Z 2 reduction 38 . We note that this interaction term is not invariant under any rotation between the flavors since it singles out a very specific channel of coupling the flavors. Thus, out of the SU(4) flavor rotation of the non-interacting theory, only a U(1)×U(1) symmetry remains corresponding to simultaneous phase rotations of Ψ A and Ψ B or Ψ C and Ψ D .

III. MEAN-FIELD THEORY & ATOMIC LIMIT
Before presenting the method and the numerical results, let us develop an intuition on symmetric and symmetry breaking phases that the model presented in the previous section admits. Here, we discuss various mean-field scenarios and identify the channel which is most likely to exhibit symmetry breaking long-range order. We also argue that this phase may only be realized when the coupling strength is comparable or larger then the band gap U c,1 < U . Additionally, we will solve the atomic limit analytically for U → ∞. We find a unique and symmetric ground state that is gapped from the rest of the spectrum. Hence this state is also stable with respect to finite values of U . Accordingly, the symmetry broken phase is bounded both from below and above with U c,1 < U < U c,2 . Interestingly, the states of the limiting cases all share the same symmetry such that an intermediate symmetry broken phase is not required. This provides another argument, complementary to the edge state analysis of Ref. 19, on the existence of an adiabatic path and the according reduction of the topological classification.

A. Weak interaction limit and Mean-field scenarios
To discuss mean-field scenarios relevant at weak interactions, it is very useful to divide the four layers into two pairs, namely (A, B) and (C, D) and introduce the two pseudo-spin operators with β = x, y where the Pauli matrices µ β act on layer index within each pair. This allows us to rewrite the interaction as and shows the possibility to minimize the energy for ζ = −, given that U > 0, by the generation of pseudo-magnetic order in the xy-plane of Note that the operator M β i contain terms like Ψ † i,A γ 5 Ψ i,B . The required R symmetry transforms γ 5 to −γ 5 and is therefore broken by the order parameter. Physically, M β i introduces single electron hopping processes, e.g., from layer B to A, that flip the R charge. We have illustrated this for the helical edge states in Fig. 2(b). This also points out that the order parameter generates a gapped edge spectrum.
Additionally, this order parameter anti-commutes with Eq. (1) such that it also introduces a gap for the bulk semimetals. As Eq. (5) points out, the interaction is symmetric under rotation around the z-axis of the pseudo-spins, and hence the orientation of M within the xy-plane is arbitrary and also breaks this symmetry spontaneously.
Without loss of generality, we choose the x direction for the magnetization and introduce H MF = m i M x i . In Fig. 2(a) we present the energy spectrum for H 0 + H MF with open boundary condition in y direction. The solid black lines represent a non-zero mean-field expectation value m > 0 and the red lines overlay the symmetric version (m = 0). We can clearly observe the introduced gap of the former massless Dirac edge state.
To make the connection between this mean-field scenario and the phase diagram, let us discuss the limiting cases. Keeping the interaction strength small, we expect stable Dirac cones for λ ∼ 0 and λ ∼ −2 as the density of states at the Fermi level vanishes at half filling [39][40][41][42][43][44] . The insulating states provide an intrinsic scale of energy, namely the band gap, such that the correlation should reach comparable strength before it leads to significant changes. Hence we expect that a symmetry broken phase, if at all, occurs at finite interaction strength U > U c,1 .

B. Strong interaction limit
The more interesting limiting case is the strongly interacting one with U/t 1. Starting from the limit t = 0, we can solve H int analytically, as the lattice sites completely decouple and we are left with a zero-dimensional problem. In the following, we calculate the full spectrum and show that there is a unique ground state. For readability, let us drop the position index i for the remainder of this analytic derivation. Note that P ± act as projectors such that the Fock space can be decomposed into two separate blocks which have an identical spectrum. Hence it is sufficient to focus on one subspace, say H + . Let i be an eigenvalue of H + with degeneracy g i . The full spectrum is then given by i + j with degeneracy g i g j .
In the previous Sec. II, we have chosen a basis for the γ-matrices in which R is diagonal to remind the reader of popular QSH models. Here, however, it is more convenient to choose a different basis in which both γ 5 = σ x τ z and iγ 3 γ 4 = σ x τ 0 are diagonal 45 . Note that the local fermion degrees of freedom H + within one layer o is then fully classified by the eigenvalues s 5 = ± of γ 5 such that P + Ψ o,s5 = Ψ o,s5 and γ 5 Ψ o,s5 = s 5 Ψ o,s5 . This leads to the definition of four new spin operators: such that the Hamiltonian can be written as H + = h ac − h ad − h bc + h bd where we used the shorthand notation h ij = U (S + i S − j + h.c.). Note that this has mapped Eq. (3) to a model with four sites on a ring (a-c-b-d) of spinful fermions. Each site can be empty, double occupied, or have a single fermion with up or down spin. The total Hilbert space is then 2 8 = 256 dimensional. The Hamiltonian conserves the parity at each site 46 as well as the total S z spin component. Using the local parities, we can decompose the Hilbert space into 16 16-dimensional Hilbert spaces which can be studied as shown below.
We start by considering the subspaces where at least both a and b or both c and d sites are parity even. There are 7 16dimensional subspaces with this property. Observe that every second site represents a spin singlet such that spin-flip processes on all bonds are prohibited. Hence, the Hamiltonian and all eigenvalues vanish.
Next, we discuss the cases where only the a or b as well as only the c or d site is single occupied which occurs for 4 16-dimensional subspaces. Then, only one bond operator, say h ac , is non-zero. From the spin conservation follows that the sectors with maximal value of |S z | have a vanishing Hamiltonian as the spin-flip operators cannot act on those states. The S z = 0 subspace contains only two eigenstates with nonzero eigenvalues ±U given by | ↑ a ↓ c ± | ↓ a ↑ c .
Third, we consider the subspace where exactly one parity is even (there are 4 such subspaces), e.g., site b, such that H + = h ac − h ad . Once more, the states with maximal |S z | = 3/2 are eigenstates of energy zero. In the S z = ±1/2 subspace we find the Hamiltonian Here the eigenvalues are given by ± √ 2U and zero. The last and most interesting subspace has only odd parity sites. Like before, the two states with S z = 2 are eigenstates with vanishing energy. In the S z = ±1 sector we find eigenvalues of 0 and ±2U . The wave function of the −2U state is (from here on out, we drop the indicies for readability) Finally, for S z = 0 there are multiple states with vanishing energy, but only one state is associated with the eigenvalues ±2 √ 2U , respectively, where the ground state wave function is In summary, the full spectrum of H + consists of 186 states of vanishing energy, 16 states with energies ±U and ± √ 2U each, 2 modes for ±2U and a unique state at energy ±2 √ 2U . It is also interesting to notice that lowest excitation has an energy of ω = 2( √ 2 − 1)U and changes the spin by ∆S z = ±1. It is straight forward to show that the lowest states φ 0 and φ 1 are related as Note that the b and d sites have negative γ 5 eigenvalues. This, combined with the second relative minus sign between (a, b) and (c, d) indicates the relation to the operator M β i defined in Eq. (6) which has already been identified as the most dangerous mean-field channel. Here, we learn that this operator also exhibits the lowest excitation of the strongly interacting limit.
We have shown that the block H + exhibit a gapped and unique ground state, which is therefore also symmetric. The overall ground state of the full lattice model is then a product state that is also gapped, unique and respects all symmetries. Allowing finite, but small values of t will change details of the ground-state wave-function. Especially it will no longer be a product state of purely local states. However, it will remain unique and symmetric until the hopping energy-scale set by t is comparable to the many-body gap. This is quite opposite to local interactions with a degenerate ground state, e.g., the regular Hubbard model on the two-dimensional square lattice, where small perturbations can generate symmetry broken states, e.g., anti-ferromagnetic order due to a finite hopping amplitude in the Hubbard model.
To summarize, the limiting cases all exhibit unique ground states that share the same symmetries. As the strongly interacting state is a local product state, it is a representation of an atomic limit and very likely adiabatically connected to the trivial band insulators. The topological insulator is also stable against small interactions that were specifically chosen to FIG. 3. Comparison of the single-particle spectrum A(ω) at the Dirac point k = (π, π) extracted by analytic continuation using the stochastic maximum entropy method, with the band dispersion (solid black line) extracted from fitting the tail of time-displace Greens function. This proves that the assumption of a single low-energy excitation is justified. allow a connection of the topological state to the strong interacting limit as the topologically protected defect states are symmetrically gapped 19 . However, especially for intermediate coupling strength along the path from the topological insulator to the Mott phase at strong interacting, the energy scales mix and spontaneous symmetry breaking might occur. The dangerous channel here is given by pseudo-magnetic order.

IV. METHOD
Why should we go beyond mean field? The considerations of Sec. III allow two scenarios for intermediate coupling strength, both are in agreement with the theorems on non-interacting systems. Either the order parameters vanishes and the bulk gap closes at the topological phase transition, or, the non-zero order parameters ensures a finite bulk gap at the expense of a broken symmetry. However the interaction constructed following the rules of Ref. 19 aims at a different path, namely an adiabatic connection of two topological distinct phases. Such a setup has to keep the band gap finite while preserving all protecting symmetries. Hence this connection cannot be made on a mean-field level and requires an intermediate state without a quasi particle description. In other words, one replaces poles of the Greens function that cross the Fermi surface (band gap closing) by zeros (no spectral weight) in order to change the topological invariant.
To solve the interacting system, we use the ALF package 31 , a general implementation of the auxiliary field Quantum Monte Carlo 28 . The zero-temperature version of this algorithm provides access to ground-state properties by using a trial wave function |ψ T , we take the non-interacting ground state, and project it to the correlated one by applying the exponentiated Hamiltonian exp(−ΘH)|ψ T 47,48 . Here Θ controls the projection length, the result converges exponentially fast and we choose Θ = 20 for the remainder of this work. The implemented auxiliary-field QMC algorithm uses the Trotter decomposition of the partition sum Z = Tr e −β(H0+Hint) = Tr [e −∆τ H0 e −∆τ Hint ] NTrotter + O(∆τ 2 ) with ∆τ = β/N Trotter as well as a discrete Hubbard-Stratonovich transformation for the interaction of perfect squares e ∆τÂ 2 = l=±1,±2 γ(l)e √ ∆τ η(l)Â + O(∆τ 4 ). The Monte Carlo weight of each configuration {l τ,i,ζ,α,β }, where τ labels the imaginary time-slice and i, ζ, α, β are the indices of Eq. (5), is determined by integrating out the fermions and is given by the determinate of a matrix W ({l τ,i,ζ,α,β }).
The simulation of this model in the above formulation does not suffer a sign problem due an anti-unitary symmetry T defined as, We note that γ 1 γ 5 commutes with the projection P α and anticommutes with γ 5 such that T −1 S 1,α,β . This symmetry is satisfied for each field configuration and squares to −1, such that the eigenvalues of W ({l τ,i,ζ,α,β }) come in complex conjugated pairs which guarantees the positivity of each configurations weight ∼ det[W ({l τ,i,ζ,α,β })]. 49 Three observables are the main focus of this study. The first derivative of the free energy ∂F/∂U = −β/U H int signals a first-order phase transition for increasing interaction strength U . To detect second-order phase transitions, we define a correlation ratio r = 1 − S(q=δq) S(q=0) where S(q) is the correlation function in momentum space and δq the smallest but finite momentum available on the given lattice size. Observe that q = 0 assumes a homogeneous instability which is justified according to the mean-field analysis. The correlation function is given as In case of long-range homogeneous order S(0) diverges linearly with system size L 2 whereas the correlation function remains finite for any other value of q, hence r = 1 for the thermodynamic limit. In systems without long-range order, S(q) is a smooth function such that r converges to zero for large lattices. As r is an RG-invariant quantity it exhibits a crossing point for different system sizes at the phase transition. The third observable of interest is the single-particle gap that allows us to track the semi-metallic Dirac cones that separate the insulators of different topology for U = 0.
The single-particle greens functions are given by for hole excitations. If we assume a single quasi-particle mode at low energies gapped from higher energy excitations, then both greens function behave as G k (τ ) ∼ a k exp(−τ k ) for large values of τ where k is the excitation energy and a k its spectral weight. As a sanity check of this assumption, we compare the extracted energies with the full spectrum (see Fig. 3) determined by MaxEnt 50 which proves that the assumption is justified.

V. RESULTS
We have shown above that the limit of strong interaction generates a gapped and symmetric ground state. Those two properties also apply to both non-interacting ground states, representing −2 < λ < 0 and 0 < λ, such that the adiabatic connection seems to be plausible, at least in principle. To test this hypothesis, we will track the semi-metallic phase and analyze the most dangerous correlation function identified in the mean-field considerations.
A. Tracking the semi-metal Note that the two Hamiltonians with ±(λ + 2) can be mapped onto each other. First, the mapping has to shift the momenta as k → k + (π, π) such that m (λ+2) (k) = −m −(λ+2) (k + (π, π)). To absorb the sign changes in the first three terms of Eq. (1a) the Dirac spinor has to transform as Ψ k,o → γ 4 Ψ k,o . Therefore, the position of the semi-metal with two Dirac cones at (π, 0) and (0, π) has to remain at fixed λ = −2.0 whereas the Dirac cone at (π, π) generically occurs at renormalized values λ ∼ 0. It is interesting to notice that the γ 4 anti-commutes with R such that the winding w → −w is inverted and the two Hamiltonian represent opposite topologies.
One might also be concerned that the interaction might lead to a meandering of the Dirac cone within the Brillouin zone. If we had kept the PHS with ∆ = 0, then the cones are symmetry constraint to the time-reversal invariant momenta. On one hand, we fine-tuned the symmetry breaking such that the Dirac cones remain gapless in the free fermion system, and on the other hand, the numerical results show that the Dirac cone remain where they are.
Let us introduce λ c (U ) as the critical value at which the semi-metal marks the topological phase transition between the TI with winding w = +4 and the trivial insulator. For the free fermion system, we have λ c (0) = 0. To locate the phase transition, we set a fixed interaction strength, e.g., U = 1.0, and scan the single particle spectrum for various values of λ. The resulting excitation energies (π,π) of the Dirac cone are presented in Fig. 4(a). The results depend on the lattice size and a visual extrapolation suggests a semi-metal at λ c (1.0) = −0.04 ± 0.02. We repeat this analysis for various values of U and also confirm the symmetry constrained position λ c (U ) = −2.0 for the Dirac semi-metal with cones at (0, π) and (π, 0) which separates the w = +4 TI from the one with winding w = −4 at −4.0 λ < −2.0. The results are summarized in panel (c).

B. Symmetry-broken phase
Here, we focus on the intermediate region of the phase diagram where the energy scales of the correlation and the kinetic energy compete with each other. During the meanfield analysis (see Sec. III), we have identified this regime between the TI and the Mott insulator at strong interactions as FIG. 4. On the left hand side, we present the extracted energies k=(π,π) for the lowest particle/hole excitation. The system size scaling suggest a gap closing for λc = −0.04 ± 0.02. In the central panel, the correlation ratio r is presented and the size scaling is consistent with a symmetry broken phase between Uc = 1.65 ± 0.02 and Uc = 3.2 ± 0.2. The right hand side summarizes various scans in the phase diagram. the one most prone to spontaneous symmetry breaking with long-range pseudo-magnetic order.
Let us start with a fixed value of λ = −0.5 and analyze the correlation ratio r with increasing interaction strength U . The resulting data is depicted in Fig. 4(b) for various lattice sizes. We clearly see that the ratio systematically decreases with increasing L if the correlation strength is smaller than U c = 1.65 ± 0.02 or larger than U c = 3.2 ± 0.2 such that there is no long-range order here. In the intermediate region, the correlation ratio increases with system size that indicates spontaneous symmetry breaking due to a finite pseudomagnetic order in xy-plane of M β i . The critical values stated above are extracted from the crossing point where the ratio coincides for all lattices. The second phase transition from the ordered phase to the Mott insulator requires larger interaction strength such that the QMC simulation become more challenging, hence the smaller lattice sizes and larger error estimate. Again, we repeat this calculation for various values of λ and summarize the phase boundary in panel (c).

C. Phase diagram
In panel (c), we plot the full phase diagram and confirm the expected stability of the Dirac semi-metals as well as the insulators at weak coupling strength. The simulations also detect the symmetric state with strong correlation. In the middle of the phase diagram, where kinetic and potential energies are of similar order, we find long-range order in exactly the dangerous channel that we have identified in Sec. III. This phase breaks the protecting R symmetry and therefore allows a hybridization of counter-propagating edge modes as shown in Fig. 2. As a result, we cannot find an adiabatic path between the two non-interaction topological insulators. Instead, any path in this phase diagram either contains a semi-metallic state or a symmetry-breaking phase.

D. Can frustration remedy the problem?
The main idea is to add the z-component of the pseudospins defined in Eq. (4) and use this to frustrate the in-plane order without changing the wave functions of the limiting cases. To form a proper SU (2) algebra, we have to drop the γ 5 matrix acting on the Dirac components as (γ 5 ) 2 = 1 such that Observe that this z-component generates rotations within the xy-plane of the pseudo-spins that leave the Hamiltonian invariant. Additionally, the transformation (A ↔ B) combined with (C ↔ D) also is a symmetry operation under which S 1/2,α,z i → −S 1/2,α,z i is inverted. Hence, any unique ground state has to be an eigenstate of i,α=± S 1,α,z i + S 2,α,z i with a vanishing expectation value. In the large U limit, the sites and α-sub-blocks decouple such that we introduce an additional interaction term H frust = V i,α=± S 1,α,z i + S 2,α,z i 2 which minimizes S z locally without changing the ground state.
As depicted in Fig. 5, weak frustration does reduce the size of the symmetry broken phase. In panel (a), we present the phase diagram for V = 0.75. The symmetry broken phase now extends only to λ ∼ 0.5 whereas in the unfrustrated model it reaches λ ∼ 0. Additionally, we find that the long-range ordered region is shifted towards weaker coupling strength U while also the range in U has been reduced. This trend is also clearly visible in panel (b) where we kept λ = −2.0 fixed. The Dirac cones at (π, 0) and (0, π) persist for weak coupling strength U and V . Increasing U generates the long-ranged ordered state before the appearance of the symmetric Mott insulator at large U . With higher level of frustration the symmetry broken phase is replaced by a direct first order phase transition between the Dirac semi-metal and the Mott insulator. In Fig. 5(c), we show that the first order phase transition extends also to λ > −2.0 and connects to FIG. 5. The left hand side show the phase diagram for weak frustration V = 0.75 with a smaller region of long-range order compared to V = 0 that seems to be most stable for λ = −2. In the middle, we present the phase diagram for fixed λ = −2.0 for various frustration strength V . For V > 2.5, the symmetry broken phase is replaced by a first order phase transition. On the right hand side, the phase diagram is presented for high levels of frustration.

VI. DISCUSSION
In this study, we found an interaction which at sufficient strength trivializes the topological phase w = 4. At the same time, surprisingly, it does not allow for an adiabatic path between this topological phase and the trivial phase w = 0. The semi-metal separating the non-interacting insulators persists for small coupling strength in U . It is terminated either by a first order transition to a symmetric Mott insulator related to the large U limit or a second-order phase transition to a long-range ordered phase. This in turn is separated by another second-order transition at larger U to the symmetric Mott insulator. Similarly, the topological insulator either undergoes a direct first order transition to the Mott phase or a second order transition into a symmetry-broken phase followed by another transition into the symmetric Mott phase.
We emphasize that our results do not contradict the statement that the non-interacting classification is reduced from Z to Z 4 with interaction. Instead, our results show that the existence of a specific interactions which symmetrically gaps out the surface states of a topological insulator is not sufficient to establish that the corresponding bulk phase can be adiabatically connected to a trivial phase. This contradicts a popular line of reasoning which focuses on gapping out the surface states as a criterion for establishing the classification reduction. The underlying reason for the failure of such arguments is that the stability of surface states can be essentially reduced to a zero-dimensional problem by studying the zero-energy states within defects constructed in a specific way 19 . However, this decreased dimensionality does have strong implications on the possibility of spontaneous symmetry breaking in the ground state that may only occur in one (two) or higher dimensions for discrete (continuous) symmetries according to Mermin-Wagner theorem 51 . And it is exactly this mechanism which blocks the path we were looking for.
From the two-dimensional bulk perspective, the phase diagram exhibits various critical points. Even though the focus of this study was to establish the phases themselves, it is interesting to discuss those critical theories briefly. There are Wilson-Fisher transitions between the topological insulator and the ordered phase as well as between the ordered phase and the Mott insulator. Additionally, we expect the critical point where the semi-metal is gapped by symmetry breaking to be described by a Gross-Neveu theory 44,52 .
The results of this work raise a few questions, mainly focused on the missing pieces required to find the adiabatic path. In Fig. 6, we sketch two alternative scenarios for the bulk phase diagram, (a) the symmetric mass generation for the Dirac cone as well as (b) a separated region of symmetry breaking that terminates the semi-metal line. Several studies have reported the formation of a correlated single-particle gap of SU (4)-symmetric Dirac system without the generation of long-range order 26,27 . Most surprisingly, it is claimed to be a second order transition. In Ref. 53, the authors propose a theory that involves fractionalization in order to explain this exotic phase transition. It would be very promising to include the same kind of bulk criticality in our setup in order to find the adiabatic connection and then investigate the details of how this affects the topological aspects. In contrast, the scenario (b) shows that this symmetric mass generation is not required and that there also exist more conventional options in which the symmetry broken region we find is split into two separate ones.
However, there is no obvious approach to engineer this phase diagram. One possibility is to consider the symmetry of the interaction. Upon using highly symmetric interactions in the flavor space, the previous studies 54 found a direct secondorder transition to the symmetric Mott phase suggesting the possibility of an adiabatic path between the topological and trivial phases. In contrast, the interaction we use has very low symmetry in the flavor space, and we always found a symmetry-broken phase or a first order transition which completely blocks such adiabatic path. It would be interesting to There exist at least two alternative bulk scenarios with (a) symmetric mass generation (green cross) or (b) a symmetry breaking region (red circle) with long-range order (LRO) terminating the semimetal between topologically distinct non-interacting states. Observe that the region with LRO connects to only one semi-metal and not two both as in our numerical phase diagram in Fig. 4(c). Both scenarios allow an adiabatic path and are not realized in the range of parameters investigated in this study.
see if interactions which are more symmetric than ours but less symmetric than those considered previously could yield a phase diagram similar to Fig. 6a).
Another possibility yet again involves the bulk-boundary correspondence, similar to previous studies, in order to finetune the energy scales of the model. The energy scale of the helical edge state is set by its Fermi velocity v edge . Interestingly, it is possible to reduce v edge without a significant change of the bulk gap 55 . This is achieved by localizing the Berry curvature in momentum space. Hence, we can control and therefore separate the bulk energy scale ∆ bulk from the edge v edge without introducing open boundary condition by using bulk parameters. This constitutes a promising avenue towards the aforementioned hierarchy of energy scales v edge < U < ∆ bulk , and this hierarchy is also often used in the arguments on the reduced topological classifications. There, the interaction is supposed to be strong enough to gap out the edge spectrum and therefore trivializes the topology, but still weak enough to avoid the broken symmetry in the bulk. Here, it is usually believed that the phase transition of the edge coincides with the topological phase transition of the bulk and thus enables the adiabatic connection.
However, this also demonstrates an important subtlety of bulk boundary correspondence in such systems. In particular, what would happen if we consider an interaction which is very weak compared to the bulk energy scale U/∆ bulk 1 so that we expect the phase to be adiabatically connected a non-interacting bulk topological phase while at the same time being rather strong compared to the characteristic energy scale at the edge U/v edge > U c edge /v edge (where U c edge the critical interaction strength). Thus, the bulk and edge phase transition do not have to coincide, at least based on the involved energy scales. There are at least two distinct scenarios in this case: (i) the edge may gap out by spontanuously breaking the protecting symmetry, (ii) it may gap-out symmetrically and become topologically trivial 56 . Whereas the latter scenario violates bulk boundary correspondence, the former does not. In this case, we expect the symmetry at the edge to be restored as the strength of the interaction is increased beyond a characteristic scale associated with the bulk rather than the edge. Investigating which of these scenarios is realized represents a very interesting extension of the current work.