Topological Invariant and Quantum Spin Models from Magnetic \pi\ Fluxes in Correlated Topological Insulators

The adiabatic insertion of a \pi flux into a quantum spin Hall insulator gives rise to localized spin and charge fluxon states. We demonstrate that \pi fluxes can be used in exact quantum Monte Carlo simulations to identify a correlated Z_2 topological insulator using the example of the Kane-Mele-Hubbard model. In the presence of repulsive interactions, a \pi flux gives rise to a Kramers doublet of spinon states with a Curie law signature in the magnetic susceptibility. Electronic correlations also provide a bosonic mode of magnetic excitons with tunable energy that act as exchange particles and mediate a dynamical interaction of adjustable range and strength between spinons. \pi fluxes can therefore be used to build models of interacting spins. This idea is applied to a three-spin ring and to one-dimensional spin chains. Due to the freedom to create almost arbitrary spin lattices, correlated topological insulators with \pi fluxes represent a novel kind of quantum simulator potentially useful for numerical simulations and experiments.


I. INTRODUCTION
A topological insulator represents a novel state of matter characterized by a special band structure that can result, e.g., from strong spin-orbit interaction [1,2]. In two dimensions, this state is called a quantum spin Hall insulator, and has deep connections with the quantum Hall effect, including the coexistence of a bulk band gap and metallic edge states, the absence of symmetry breaking, and the possibility of a mathematical classification [3,4]. Importantly, because of the absence of a magnetic field, the quantum spin Hall insulator preserves time-reversal symmetry which provides protection against interactions and disorder [5][6][7]. The quantum spin Hall insulator has been realized in HgTe quantum wells [8,9].
Correlated topological insulators with strong electronelectron interactions are the focus of current research [10]. Intriguing concepts include electron fractionalization in the presence of time-reversal symmetry [11][12][13][14], spin liquids [14][15][16], and topological Mott insulators [17,18]. Remarkably, some of the theoretical models can be studied using exact numerical methods. A central problem in this context is the question of how to detect a topological state directly from bulk properties, for example in cases where the bulk-boundary correspondence breaks down. Experimentally, this issue also arises in the absence of sharp edges in proposed cold atom realizations as a result of the trapping potential [19,20]. The classification in terms of a Z 2 Chern-Simons index relies on Bloch wavefunctions and is therefore only valid for noninteracting systems. Generalizations involve twisted boundary conditions [21] or Green functions [22][23][24][25][26][27], and are challenging to use in experiments or exact simulations. Indirect signatures such as the closing of gaps [16] or the crossing of energy levels [28] require, among other difficulties, experimental tuning of microscopic parameters.
Topological insulators show a unique response to topo-logical defects such as dislocations [29,30] or π fluxes [12,30,31]. Upon adiabatic insertion of a π flux, Faraday's law together with the quantized transverse conductivity gives rise to midgap charge and spin fluxon states [12,31]. These states are exponentially localized around the flux [12,31]. The existence of these states is ensured even in the presence of interactions or disorder by time-reversal symmetry, and has been suggested as a bulk probe of the Z 2 index [12,31]. The concept of fluxons can also be generalized to situations where spin is not conserved, such as in the presence of Rashba coupling. In three dimensions, a magnetic flux gives rise to the wormhole effect [32]. Electron-electron repulsion lifts the degeneracy of charge and spin fluxons, but the two degenerate spin fluxon states constitute a localized spin with S z = ±1/2 [12]. Dynamical π fluxes emerge in the context of fractionalized topological insulators [12,13]. Previous work on π fluxes in noninteracting quantum spin Hall insulators [12,30,31] was based on squarelattice models such as that for HgTe quantum wells [8].
Here we consider the half-filled Kane-Mele model on the honeycomb lattice [3], historically the first model with a Z 2 topological phase, which has close connections to graphene [3], the integer quantum Hall effect [33], and when including interactions, to correlated Dirac fermions [15]. Topological phases of interacting fermions on honeycomb lattices may be realized in transition metal oxides [34], semiconductor structures [35], graphene [36], or cold atoms [37], see also Ref. [10].
Here we use π fluxes in combination with exact quantum Monte Carlo simulations, and show that they can be used efficiently to probe the topological invariant of correlated topological insulators. In particular, this method does not rely on an adiabatic connection to a noninteracting state, and may also be used for fractional states. In addition, we demonstrate that π fluxes permit the construction of quantum spin models of almost arbitrary FIG. 1. (Color online) For a lattice with periodic boundaries, π fluxes can be inserted in pairs. Each flux threads a hexagon (highlighted in blue) of the honeycomb lattice, and the pair is connected by a branch cut (blue line). Hopping processes crossing the branch cut acquire a phase e iπ = −1.
geometry and with tunable, dynamical interactions. The spins correspond to the spin fluxons created by inserting π fluxes, and the interaction is mediated by magnetic excitons corresponding to collective magnetic fluctuations of the topological insulator. These spin models can be studied theoretically with the quantum Monte Carlo method, or experimentally. As examples, we show that a ring of three spins has a ground state with magnetization 1/2, and that a one-dimensional chain of fluxons undergoes a Mott transition and is described at low energies by an XXZ model.
The article is organized as follows. In Sec. II, we introduce the Kane-Mele and Kane-Mele-Hubbard models. Section III provides details about the methods. The use of π fluxes as a probe for topological states is discussed in Sec. IV, whereas the construction of quantum spin models is the topic of Sec. V. Conclusions are given in Sec. VI, and we provide three appendices.

II. MODEL
The half-filled Kane-Mele model with additional electron-electron interactions can be studied with powerful quantum Monte Carlo methods [16,38]. Using the spinor notationĉ † i = ĉ † i↑ ,ĉ † i↓ , whereĉ † iσ is a creation operator for an electron in a Wannier state at site i with spin σ, the Hamiltonian reads The notation i, j and i, j indicates that the sites i and j are nearest neighbors and next-nearest neighbors, respectively, and implicitly includes the Hermitian conjugate terms.
The first term describes the hopping of electrons between neighboring lattice sites. The second term represents the spin-orbit coupling which reduces the SU (2) spin rotation symmetry to a U (1) symmetry. The third term is an additional Rashba spin-orbit coupling [39]. The additional factors τ ij = ±1 take into account any π fluxes present, whereas the original Kane-Mele model (without π fluxes) is recovered from Eq. (1) by setting τ ij = 1.
The spin-orbit term corresponds to a next-nearest neighbor hopping with a complex amplitude iλ, and has been derived from the spin-orbit coupling in graphene [3]. This hopping acquires a sign ±1 depending on its direction, the sublattice, and the electron spin. This sign is encoded in (ν ij · σ), where d ik is the vector connecting sites i and k, and k is the intermediate lattice site involved in the hopping process from i to j. For a coordinate independent representation, the vectors d αβ are defined in three dimensional space, although the z component vanishes. The vector σ is defined by σ = (σ x , σ y , σ z ), with the Pauli matrices σ α . The last term in Eq. (1) is the Rashba spin-orbit interaction [3,5]. It is defined in terms of the spin vector s = σ/2, and the unit vectord ij which can be expressed in terms of the nearest-neighbor vectors δ 1 , δ 2 , δ 3 [40]. The Rashba coupling breaks the z → −z inversion symmetry, and has to be taken into account, for example, in the presence of a substrate. Because this term includes spin-flip terms, spin is no longer conserved. The Rashba term has been included in the results for the noninteracting model (1), but cannot be included in quantum Monte Carlo simulations of the interacting model (3) due to a minus-sign problem.
The model (1) can be solved exactly [3,5,41]. In the absence of Rashba coupling, α = 0, the Kane-Mele model describes a Z 2 quantum spin Hall insulator for any λ > 0. This state is characterized by a bulk band gap ∆ sp = 3 √ 3λ, a spin gap ∆ s = 2∆ sp , and a quantized spin Hall conductivity σ s xy = e 2 2π . The topological state survives for Rashba interactions α < 2 √ 3λ (for chemical potential µ = 0), and has protected, helical edge states for geometries with open boundaries [3,5,41]. We use t as the unit of energy ( = 1), take λ/t = 0.2, and consider periodic lattices with L × L unit cells.
To investigate the effect of electron-electron repulsion, we consider the paradigmatic Hubbard interaction [42] and arrive at the Kane-Mele-Hubbard model [43], Hamiltonian (3) without Rashba coupling has been studied intensely [16,38,[43][44][45][46][47][48]. In particular, its symmetries permit the application of exact quantum Monte Carlo methods without a sign problem [16,38,48]. On a lattice with periodic boundaries, π fluxes can only be inserted in pairs, as illustrated for the minimal number of two fluxes in Fig. 1. The flux pair is connected by a branch cut (or string), and every hopping process

III. METHOD
We have used the auxiliary-field quantum Monte Carlo method [49], which was previously applied to the Hubbard model on the honeycomb lattice [15], and the Kane-Mele-Hubbard model [16,38,48]. The central idea of this stochastic method is to use a path integral representation of the interacting model (3). By means of a Hubbard-Stratonovich transformation, the Hubbard term is decoupled, leading to a problem of noninteracting fermions in an external, space and imaginary-time dependent field. The sampling is over different configurations of these auxiliary fields in terms of local updates. For a given configuration of fields, Wick's theorem can be used to calculate arbitrary correlation functions from the singleparticle Green function. We refer to a review [50], and previous work [15,16,48] for technical details such as the calculation of energy gaps.
Here we have used a projective formulation (with projection parameter θt = 40) to obtain ground state results starting from a trial wavefunction (the ground state of the U = 0 case) [48], and a finite-temperature formulation to calculate thermodynamic properties. Both variants rely on a Trotter discretisation of imaginary time (we used ∆τ = β/L = θ/L = 0.1), but the associated systematic error is smaller than the statistical errors. At half filling, time-reversal invariance ensures that simulations can be carried out without a minus-sign problem even in the presence of π fluxes.  2) with two π fluxes at the maximal distance on an 18 × 18 lattice, for different Rashba couplings α. At temperatures kBT 0.1t, each π flux contributes 1 2k B T to the susceptibility, leading to χs ≈ 1 k B T . Also shown is the spin gap energy scale 0.2∆s for T = 0, α = 0. For α > 0, the chemical potential is adjusted to retain a half-filled band.

IV. USING π FLUXES TO PROBE CORRELATED TOPOLOGICAL STATES
A. Thermodynamic signature of π fluxes In the topological phase of the model (1), each π flux gives rise to four fluxon states which are exponentially localized (due to the bulk energy gap ∆ sp ) near the corresponding flux-threaded hexagons [12,31], see Fig. 1. The states correspond to the spin fluxons |↑ , |↓ with energy E ↑↓ , forming a Kramers pair related by time reversal, and the charge fluxons |+ , |− (with energies E + , E − ) related by particle-hole symmetry. As we show in Fig. 2, the fluxon states lie inside the bulk band gap, and for noninteracting electrons The fluxons leave a characteristic signature in the static spin and charge susceptibilities which are defined in terms of the operators of total spin in the z direction,M z = iĉ † i σ zĉ i , and of the total charge,N = iĉ † iĉi ; the inverse temperature is given by β = 1 kBT . At low temperatures k B T ∆ sp , we can restrict the Hilbert space to {|↑ , |↓ , |+ , |− }. If the spin fluxons are independent, and for α = 0, we expect a Curie law χ s = χ c = 1 2kBT per π flux, and hence χ s = χ c = 1 kBT for two independent π fluxes (see Appendix A). The prefactor of the Curie law follows from the quantized spin Hall conductance in the absence of Rashba coupling [3]. Similarly, a Curie law was also predicted for topological excitations in polyacetylene [51]. Figure 3 shows results for χ s as a function of temperature for the Kane-Mele model with two π fluxes ∆ sp /t ,  located at the largest possible distance. At temperatures k B T ≈ ∆ s , χ s is dominated by bulk effects. For k B T 0.1t, we observe the expected Curie law. The latter is robust with respect to Rashba coupling, which is crucial for possible experimental realizations. Figure 3 establishes the existence and thermodynamic signature of degenerate spin and charge fluxons in a quantum spin Hall insulator threaded by a pair of π fluxes. We now consider the effect of electron-electron interaction in the framework of the Kane-Mele-Hubbard model (3). For λ > 0, the phase diagram of the latter includes a correlated quantum spin Hall insulating phase that is adiabatically connected to that of the Kane-Mele model (i.e., U = 0), and a Mott insulating phase with longrange antiferromagnetic order [43,48]. Figure 4(a) shows the quantum phase transition between these two phases as a function of U/t at λ/t = 0.2. At the transition, the spin gap ∆ s -as obtained from finite-size scaling (see ref. 48 for details)-closes, corresponding to the condensation of magnetic excitons [47,48]. The magnetic order is of the easy-plane type, and the transition has 3D XY universality corresponding to the ordering of local mo- ments [47,48]. For U ≥ U c , time-reversal symmetry is spontaneously broken, and the single-particle gap ∆ sp remains open across the transition [48], see Fig. 4(a). The location of the critical point can be estimated from the scaling behavior of the real-space spin-spin correlation function

B. Probing correlated topological insulators
at the largest distance r = L/2. Here we consider correlations between spins on the A sublattice, but results are the same for the B sublattice. The limit lim L→∞ S xx (L/2) is identical to m 2 , with m being the magnetization per site. This critical value can be obtained by considering the 3D XY scaling behavior at the transition. Following ref. 48, we plot L 2β/ν S xx (L/2) as a function of U for different system sizes using the critical exponents z = 1, ν = 0.6717(1) and β = 0.3486(1) [52]. The well-understood magnetic transition of the model (3) provides a test case for the use of π fluxes to probe a correlated quantum spin Hall state, as well as to track the interaction-driven transition to a topologically trivial state. We solve the interacting model with two π fluxes using exact quantum Monte Carlo simulations. Spin fluxons can be detected by calculating the lattice-  site resolved dynamical spin structure factor at T = 0, defined as Here H KMH |n = E n |n , and |0 denotes the ground state. S(i, ω) corresponds to the spectrum of spin excitations at lattice site i. The value of S Ω (i) is about three orders of magnitude smaller at lattice sites which are further away from a flux, so that the spin fluxons can easily be detected numerically. In Fig. 5(b), we show results for the magnetic insulating phase at U/t = 6. As expected for this topologically trivial state, no well-defined spin fluxons exists.
The dependence of S Ω (i) on U/t across the magnetic quantum phase transition is shown in Fig. 5(c). A clear signal is found deep in the topological insulator phase, whereas a strong drop is observed on approaching the critical point at U c /t = 5.70 (3). Hence, the spin fluxon signal can be used in quantum Monte Carlo simulations to distinguish topological and nontopological phases.
As for the noninteracting case (Fig. 3), the spin fluxons created by the π fluxes give rise to a characteristic Curie law in the spin susceptibility. Figure 6(a) shows quantum Monte Carlo results for the spin and charge susceptibilities defined in Eq. (4) in the topological phase (U/t = 4). We again consider two π fluxes at the maximal separation. At low temperatures, k B T ∆ s , χ s is well described by χ s = 2 kBT , or 1 kBT per π flux. The additional factor of 2 compared to the case U = 0 comes from the splitting of spin and charge states which only leaves the Kramers doublet {|↑ , |↓ } at low energies (see Appendix A). The Curie law holds down to the lowest temperatures considered in Fig. 6(a). Finally, the charge susceptibility is strongly suppressed at low temperatures, and reveals the absence of low-energy charge fluxons as a result of the Hubbard repulsion.

C. Interaction between spin fluxons
So far, we have exploited the thermodynamic and spectral signatures of independent spin fluxon excitations (i.e., free spins). On periodic lattices, spin fluxons can only be created in pairs, and it is therefore interesting to consider their mutual interaction. Such interactions will play a key role in Sec. V, where we consider quantum spin Hall insulators with multiple π fluxes to create and simulate systems of interacting spins.
Interaction effects due to a coupling between two spin fluxons in a lattice with one pair of π fluxes become visible for larger U/t, i.e., closer to the magnetic transition. Figures 6(b) and (c) show a deviation from the Curie law below a temperature scale determined by the interaction between spin fluxons. In the model (3), this interaction is mediated by the exchange of collective spin excitations (magnetic excitons), which are the lowest lying excitations of the correlated topological insulator phase, and evolve into the gapless Goldstone mode of the magnetic state. Since magnetic order is of the easy-plane type, the dominant contribution of the resulting interaction is expected to have the general form fluxon at position r at time τ , D(r, τ ) is the Fourier transform of the exciton propagator D(q, iΩ m ) (q: momentum, Ω m = 2nπ/β: bosonic Matsubara frequency), and g is a coupling constant. At long wavelengths, the dispersion relation of the collective spin mode can be written as ω(q) = v 2 |q − Q| 2 + ∆ 2 s , where v is the spin velocity, ∆ s is the spin gap, and Q is the magnetic ordering wavevector. The minimal exciton energy is given by ω(q = Q) = ∆ s . Fourier transformation of the propagator (see Appendix C) gives in the limit of low energies and long wavelengths The first term determines the sign of the interaction. The decay at large imaginary time τ is governed by the spin gap ∆ s . The fast decay as a function of |r| underlies the clear Curie law seen, e.g., in Fig. 6. The interaction range and strength can be tuned via the spin gap and hence [cf. Fig. 4(a)] by varying U/t. From Eq. (8) we expect the interaction range to increase with increasing U/t due to the decrease of ∆ s , see Fig. 4(a). Indeed, Figs. 6(b) and (c) reveal an enhanced effect of the spin fluxon separation at low temperatures with increasing U/t. In particular, for U/t = 5.5 (close to the magnetic transition), Fig. 6(c) shows a Curie law corresponding to two free spin fluxons only for the largest system sizes (L = 18). As U → U c , the interaction range diverges, and free spin fluxons can no longer exist. For U > U c , time-reversal invariance is broken and π fluxes do not create spin fluxons. Instead, the spin susceptibility in Fig. 6(d) is that of an antiferromagnet; the finite value of χ s /L 2 at T = 0 reflects the density of spin wave excitations. To illustrate the dependence of the interaction strength on ∆ s , we consider two fluxes at a fixed, small distance as illustrated in the inset of Fig. 7. We show the spin susceptibility for different values of U/t in Fig. 7. For U/t = 3, a Curie law χ s ≈ 2 kBT may be inferred at temperatures k B T ≈ 0.1t. Increasing U/t, the interaction between the spin fluxons becomes too large to observe free spin fluxons below the temperature range set by the bulk spin gap ∆ s . The downturn of χ s occurs at higher and higher temperatures with increasing U/t, and reflects a tunable, correlation-induced energy scale for the interaction between spin fluxons that is absent in Fig. 3.

V. π-FLUX QUANTUM SPIN MODELS
The possibility of inserting π fluxes to create localized spin fluxons with a tunable interaction mediated by magnetic excitons provides a toolbox to engineer interacting spin models in correlated topological insulators. The computational effort for quantum Monte Carlo simulations does not depend on the number of π fluxes, and the latter can be arranged in almost arbitrary geometries on the honeycomb lattice.
A. Three-spin system As a first extension of the two-spin cases considered so far, we consider four spin fluxons emerging from two pairs of π fluxes. The fluxes are arranged so that three spin fluxons form a ring, and the fourth spin fluxon is located at the largest distance from the center of the ring. For large enough lattices, the separated spin fluxon will not couple to the other three, and the physical problem is π π π π π π a π 10 −2 10 −1 similar to experiments on coupled quantum dots [53] or flux qubits [54] in the context of quantum computation. The three spin fluxons experience a transverse interaction of the form (7), and behave as an effective spin with S z = ±1/2 at low temperatures (see Appendix B). The spin susceptibility for U/t = 4 shown in Fig. 8 reveals that, at low temperatures, the two independent spins indeed give rise to the expected Curie law χ s = 2 kBT . At higher temperatures k B T ≈ 0.1t, we find χ s = 4 kBT , corresponding to four independent spin fluxons. In the regime where χ s = 2 kBT , the sign of the interaction between the spin fluxons determines the ground state degeneracy of the three-spin cluster. A net ferromagnetic interaction results in a spin-1/2 doublet, whereas an antiferromagnetic coupling gives rise to a four-fold degenerate, chiral ground state (see Appendix B). In principle, the sign of the exchange coupling for the case of Fig. 8 can be determined from entropy measurements. Since Q = 0 for the model (3), Eq. (8) suggests that the interaction is ferromagnetic.

B. Simulation of one-dimensional fluxon chains
Whereas the study of systems with a small number of spins is relevant for applications such as quantum computing, many questions in condensed matter physics are related to periodic spin lattices. In this section, we therefore consider one-dimensional chains of π fluxes in the honeycomb lattice with periodic boundary conditions.
We begin with the noninteracting Kane-Mele model with a periodic flux chain. The fluxon excitations are visible in Fig. 9 which shows the integrated local density of states A Ω (i) = Ω 0 dωA(i, ω); the single-particle spectral function A(i, ω) is defined as usual in terms of the singleparticle Green function, A(i, ω) = −Im G(i, ω). Whereas the fluxons are well localized in the direction normal to the chain, the overlap of neighboring fluxons in the chain gives rise to a tight-binding band inside the topological  Fig. 9). Here λ/t = 0.2, α = 0, and the honeycomb lattice has dimensions 72 × 12. Points correspond to eigenvalues, and lines to the band structure (k) = ±2 t sin(2ka) with t ≈ 0.126t and a ≡ 1.
band gap which can be seen in the spectrum shown in Fig. 10. The specific form of the band structure can be attributed to the fact that the smallest unit cell for the fluxon chain contains two flux-threaded hexagons (and is four hexagons wide), see Fig. 9. Exploiting that the four possible fluxon states per hexagon, {|↑ , |↓ , |+ , |− }, can formally be written in terms of the fermion Fock states {|↑ , |↓ , |0 , |↑↓ }, and assuming nearest-neighbor hopping, a suitable Hamiltonian is given by where φ, ψ refer to the two flux-threaded hexagons in the unit cell, and i numbers the unit cells. The resulting band dispersion (k) = ±2 t sin(2ka) matches the low-energy bands in the spectrum (Fig. 10). The form of the effective low energy Hamiltonian and especially the gapless nature of the spectrum stems from the fact that the unit cell is a gauge choice; a translation by half a lattice vector, a π /2, corresponds to a gauge transformation. This symmetry allows the intercell and intracell hopping integrals to differ only by a phase e iθ . Imposing time-reversal symmetry pins the phase factor to θ = 0 and θ = π, thus leading to the dispersion relations ±2 t cos [(k + θ/a π )a π /2]. The choice θ = π produces the above mentioned dispersion relation, and the choice θ = 0 merely corresponds to translating the reciprocal lattice by half a reciprocal lattice unit vector. In contrast to the helical edge states of a quantum spin Hall insulator, each of the two fluxon bands is spin degenerate. As a result, and because the system is half filled, we expect a Mott transition of charge fluxons for any nonzero electron-electron repulsion. Figure 11 shows the spin and charge susceptibilities of the Kane-Mele-Hubbard model on L × 12 lattices with L/2 π fluxes and U/t = 4. The Hubbard U causes an exponential sup- pression of the charge susceptibility at low temperatures [see Fig. 11(b) and inset], whereas low-energy spin fluxon excitations remain [ Fig. 11(a)]. Hence, similar to the onedimensional Hubbard model, the fluxon chain undergoes a Mott transition to a state with a nonzero charge gap but gapless spin excitations.
In the Mott phase of the fluxon chain, the low-energy physics is expected to be described by spin fluctuations, and hence by an effective spin model with spins corresponding to Kramers doublets of localized spin fluxons. Because the interaction range depends exponentially on the spin gap, we expect nearest-neighbor interactions J xy , J zz between spin fluxons to dominate, except for the close vicinity of the magnetic transition. As argued before, the magnetic exciton is of predominantly easyplane type, and we therefore expect anisotropic interactions, |J xy | |J zz |. The minimal model for the spin chain is the one-dimensional XXZ Hamiltonian, Using the ALPS 1.3 implementation [55] we can simulate this model in the stochastic series expansion repre- sentation to calculate the spin susceptibility as a function of temperature. There is one free parameter, J xy /J zz , which is varied to obtain the best fit to the lowtemperature susceptibility (at high temperatures, bulk states of the topological insulator begin to play a role) of the Kane-Mele-Hubbard model. For example, considering six spins, a rather good match between the spin fluxon data and the XXZ model is obtained for J zz /|J xy | = −0.1 (the sign of J xy is irrelevant), see Fig. 12. Importantly, taking the same parameters, and simulating ten spins with both spin fluxons and the XXZ model, equally good agreement is found in Fig. 12. These results demonstrate that the spin fluxons form a onedimensional spin system with well-defined interactions, and that a quantum spin Hall insulator with π fluxes can indeed be used as a quantum simulator.

VI. CONCLUSIONS
In this work, we have presented quantum Monte Carlo results for a correlated quantum spin Hall insulator with topological defects in the form of π fluxes. Such fluxes represent a universal probe for the topological index that can be used in the presence of electronic correlations, and does not rely on spin conservation or an adiabatic connection to a noninteracting topological insulator. Our results demonstrate that π fluxes can be combined with exact numerical simulations, and lead to clear signatures of nontrivial topological properties in spectral and thermodynamic properties. As a concrete example, we have studied the magnetic quantum phase transition of the Kane-Mele-Hubbard model at which time-reversal symmetry is spontaneously broken. In principle, π fluxes can also be used in connection with fractional quantum spin Hall states.
More generally, π fluxes in correlated topological insulators allow one to construct and simulate quantum spin models, and hence lead to a novel class of quantum simulators. This finding is not restricted to the Kane-Mele-Hubbard model considered here. In particular, magnetism driven by electronic correlations-the origin of the interaction between spin fluxons-is a common phenomenon. The physics described here relies on the coexistence of magnetic correlations and time-reversal symmetry, and cannot be captured by static mean-field descriptions. The spin models share the topological protection of their host against, for example, disorder. In general, they are characterized by a dynamical, time-dependent interaction reminiscent of spin-boson problems. The detailed form and sign of the interaction, whose strength and range can be tuned via the electronic interactions, depends on the electronic Hamiltonian and the lattice geometry of the underlying topological insulator. Because of the spin-orbit interaction, the spin symmetry is U (1), and-similar to cold atom realizations of the quantum Ising model [56]-the spin-spin interaction is generically anisotropic. We have provided explicit evidence for the feasibility of our idea in terms of simulations of two and four spins, as well as of one-dimensional spin chains. With additional Rashba terms, spin models with a discrete Z 2 Ising symmetry result. Although spin fluxon states are still well defined [12,31], it is a priori not clear which operators have to be measured in the numerical simulations. Finally, the concept of fluxons originating from π fluxes carries over to three-dimensional topological insulators [12,32].
An open question of central importance is whether the use of π fluxes will enable us to study quantum spin systems which are currently not accessible to numerical methods, for example due to a sign problem in the presence of frustrated interactions. Whereas we have provided evidence for the possibility to simulate arrays and chains of quantum spins, and to tune the interaction strength and range, entropy measurements are required to determine the sign of the interactions. However, the latter are extremely demanding to carry out using discrete-time quantum Monte Carlo methods. A systematic effort to study spin fluxon chain and ladder geometries is currently in progress.
Our idea may potentially also be used in experiments. A strongly correlated topological insulator on the honeycomb lattice may be realized with Na 2 IrO 3 [34], or with molecular graphene [57]. It has been suggested that π fluxes can be created in a quantum spin Hall insulator by means of an adjacent superconductor and a magnetic field [31]. This idea can be generalized to arrays of π fluxes using Abrikosov lattices. Alternatively, π fluxes may be realized using SQUIDs. A potential problem is that the diameter of the π fluxes will not be of the order of the lattice constant. Other exciting recent proposals which are relevant for the realization of our idea include artificial semiconductor honeycomb structures [35], cold atoms in optical lattices [19], and cold atoms on chips [58]. In solid state setups, π fluxes can also be created by dislocations [29,30] or wedge disclinations [59].
For U k B T , the spin fluxons |↑ , |↓ are the only low-energy excitations, and χ s can be calculated by restricting ψ to {↑, ↓}. Since E ↑ = E ↓ = E ↑↓ due to timereversal symmetry, we get For the case of two independent π fluxes, the above results imply χ s = 1 kBT at U = 0, and χ s = 2 kBT for U > 0. This agrees with the numerical results shown in Fig. 3 for U = 0, and in Figs. 6,7 for U > 0.
Our derivation is only valid in the absence of Rashba spin-orbit coupling α. However, the numerical results in Fig. 3 show that the low-temperature Curie law in χ s is the same also for α = 0. The results for the Kane-Mele-Hubbard model with four π fluxes shown in Fig. 8 reveal a 2 kBT Curie law at low temperatures, and a 4 kBT Curie law at higher temperatures. This finding can be understood as corresponding to either two or four noninteracting spins. The latter case corresponds to the spatially separate spin fluxon and an effective spin-1/2 Kramers doublet (formed by the three nearby spin fluxons) in the regime where χ s ≈ 2 kBT , and to four noninteracting spin fluxons in the regime where χ s ≈ 4 kBT . The cluster formed by the three nearby spin fluxons has the possible configurations where M z denotes the total spin in the z direction.
Since the exciton-mediated interaction in the Kane-Mele-Hubbard model has the form given in Eq. (7) and hence promotes spin-flips, the ground state can be expected to have |M z | = 1/2. The above mentioned effective spin-1/2 then corresponds to the two possible values M z = ±1/2.
The degeneracy of the ground state depends on the sign of the interaction. Considering M z = 1/2, we have the allowed states |↓↑↑ , |↑↓↑ and |↑↑↓ . The spinflip terms which connect these states are of the form J(S + i+1 S − i +S + i S − i+1 ), with periodic boundary conditions. An equivalent representation is given by the Hamiltonian which describes the hopping of a particle (the spin down) on a three-site ring, with |1 = |↓↑↑ etc. The eigenstates are obtained by Fourier transformation, and have the form The eigenvalues are given by E(k) = 2J cos k .
For J < 0, the ground state has k = 0 and energy 2J. For J > 0, the ground state is chiral, with k = ±2π/3 and energy −J. Taking into account the sector M z = −1/2, we find a total ground state degeneracy of two in the ferromagnetic case J < 0, and four in the antiferromagnetic case J > 0.