Quantization of the Interacting Hall Conductivity in the Critical Regime

The Haldane model is a paradigmatic 2d lattice model exhibiting the integer quantum Hall effect. We consider an interacting version of the model, and prove that for short-range interactions, smaller than the bandwidth, the Hall conductivity is quantized, for all the values of the parameters outside two critical curves, across which the model undergoes a ‘topological’ phase transition: the Hall coefficient remains integer and constant as long as we continuously deform the parameters without crossing the curves; when this happens, the Hall coefficient jumps abruptly to a different integer. Previous works were limited to the perturbative regime, in which the interaction is much smaller than the bare gap, so they were restricted to regions far from the critical lines. The non-renormalization of the Hall conductivity arises as a consequence of lattice conservation laws and of the regularity properties of the current–current correlations. Our method provides a full construction of the critical curves, which are modified (‘dressed’) by the electron–electron interaction. The shift of the transition curves manifests itself via apparent infrared divergences in the naive perturbative series, which we resolve via renormalization group methods.


Introduction
One of the remarkable features of the Integer Quantum Hall Effect (QHE) is the impressive precision of the quantization of the plateaus observed in the experiments. While the experimental samples have a very complex microscopic structure, depending on a huge number of non-universal details related to molecular forces and the atomic structure, the conductance appears to be quantized at a very high precision, and the result only depends on fundamental constants. The understanding of this phenomenon, via a connection between the Hall conductivity and a topological invariant [4,40] was a major success of theoretical condensed matter in the 80s. The argument was later generalized to non-interacting disordered systems [1,5,9,10] and to clean multi-particle systems [3,37]: however, the definition of conductivity in the interacting case required the presence of an unphysical averaging over fluxes, expected to be unimportant in the thermodynamic limit, but a proof remained elusive for many years. Arguments based on Ward Identities for Quantum ElectroDynamics in (2 + 1)-dimensions [12,24,30], or on the properties of anomalies [15], offered an alternative view on the QHE: they indicated that quantization should persists in the presence of many body interaction, but such conclusions were based on manipulations of divergent series, or of effective actions arising in a formal scaling limit.
The problem of a mathematical proof of the quantization of the Hall conductivity in the presence of many-body interactions remained open for several years. After the works [1,3,5,9,10,37], it was dormant for more than a decade, and then, in recent years, it was actively reconsidered again, starting from [27], which proved the quantization of the Hall conductance of an interacting electron system using quasi-adiabatic evolution of the groundstate around a flux-torus, under the assumption of a volume-independent spectral gap. In [22] we followed a different approach, and proved the quantization of the interacting Hall conductivity by writing it as a convergent series, and by showing that all the interaction corrections cancel exactly, thanks to Ward Identities. Our result holds for interacting fermionic Hamiltonians of the form H 0 + U V, where H 0 is quadratic and gapped, with the chemical potential in the middle of a spectral gap of width 0 , V is a many body interaction, and |U | 0 ; this smallness condition is assumed to ensure the convergence of the power series expansion in U of the Euclidean correlations. The same result also follows from [27], in combination with a proof of the stability of the spectral gap for such fermionic Hamiltonians [13,26]. See also [6,7] for alternative proofs of the main theorem in [27]. Recently, the bulk-edge correspondence for a class of weakly interacting fermionic systems displaying single-mode chiral edge currents was also proved [2].
Given these results on the quantization of the Hall conductivity in weakly interacting systems (i.e., with interaction strength smaller than the gap), one naturally wonders what happens for stronger interactions. We focus on the interacting extension of the spinful Haldane model [23], which has been recently realized in cold atoms experiments [31] and can be used as the building block of more general topological insulators [25]. Extensions to related systems is straightforward, in particular to the interacting, spin-conserving, Kane-Mele model, for which the quantization of the edge conductivity has been recently established [34]. We recall that, in the absence of interactions, the phase diagram of the spinful Haldane model consists of two 'trivial' insulating phases, with vanishing transverse conductivity, and two quantum Hall phases, with transverse conductivity σ 12 = ± 2e 2 /h, separated by two critical curves. By [22], we know that, away from the critical lines, for interactions U smaller than the spectral gap, the Hall conductivity is quantized and independent of U . However, what happens close to the critical lines, in cases where the interaction is much larger than the gap? This question, and in particular the possible emergence of new quantum phases, for the model at hand, via rigorous renormalization group methods. In Sect. 5 we put things together and complete the proof of our main result.

The Model
The Haldane model describes spinless fermions on the honeycomb lattice hopping on nearest and next-to-nearest neighbours, in the presence of a transverse magnetic field, with zero net flux through the hexagonal cell, and of a staggered potential. In this section we introduce an interacting, spinful, version of the Haldane model. Note that, in the presence of interactions, the spin could induce a qualitatively different behaviour, as compared with the spinless case (this is a well known fact in the context of one-dimensional fermions [14], including the edge theory of 2D topological insulators [2,34]). Note also that the experimental realization of the interacting Haldane model involves, indeed, spin-1/2 particles, see [31].
Let = x | x = n 1 1 + n 2 2 , n i ∈ Z} ⊂ R 2 be the infinite triangular lattice generated by the two basis vectors 1 . Given L ∈ N, we also let L = /L be the corresponding finite periodic triangular lattice of side L, which will be identified with the set L = x | x = n 1 1 + n 2 2 , n i ∈ Z ∩ [0, L) with periodic boundary conditions. The lattice is endowed with the Euclidean distance on the torus, denoted by | x − y| L = min m∈Z 2 | x − y + m 1 1 L + m 2 2 L|. The number of sites of L is | L | = L 2 . The periodic honeycomb lattice can be realized as the superposition of two periodic triangular sublattices A L ≡ L , B L ≡ L + e 1 , with e 1 = (1, 0) the first Euclidean basis vector. Equivalently, we can think the honeycomb lattice as a triangular lattice, with two internal degrees of freedom corresponding to the A, B sublattices.
It is convenient to define the model in second quantization. The one-particle Hilbert space is the set of functions h L = { f : L × {↑, ↓} × {A, B} → C} C L 2 ⊗ C 4 . We let the fermionic Fock space F L be the exterior algebra of h L . Notice that for fixed L, F L is a finite-dimensional space. For a given site x ∈ L , we introduce fermionic annihilation operators ψ x,ρ,s , with ρ ∈ {A, B} the sublattice label and s ∈ {↑, ↓} the spin label, and we denote by ψ † x,ρ,s their adjoint, the creation operators. They satisfy the standard canonical anticommutation relations {ψ † x,ρ,s , ψ y,ρ ,s } = δ ρ,ρ δ s,s δ x, y and {ψ † x,ρ,s , ψ † y,ρ ,s } = {ψ x,ρ,s , ψ y,ρ ,s } = 0. The operators ψ x,ρ,s are consistent with the periodic boundary conditions on L , ψ x+n 1 L+n 2 L,ρ,s = ψ x,ρ,s .
The reciprocal lattice * L of L is the triangular lattice generated by the basis vectors G 1 , . We define the finite-volume Brillouin zone as . We define the Fourier transforms of the fermionic operators as: The Hamiltonian of the model is: H = H 0 +U V, with H 0 the noninteracting Hamiltonian and V the many-body interaction of strength U . We have: with γ 1 = 1 − 2 , γ 2 = 2 , γ 3 = − 1 and n x,ρ = s=↑,↓ ψ † x,ρ,s ψ x,ρ,s , with ρ ∈ {A, B}. For definiteness, we assume that t 1 > 0 and t 2 > 0. The term proportional to t 1 describes nearest neighbor hopping on the hexagonal lattice. The term proportional to t 2 describes next-to-nearest neighbor hopping, with the complex phases e ± iφ modeling the effect of an external, transverse, magnetic field. The term proportional to W describes a staggered potential, favoring the occupancy of the A or B sublattice, depending on whether W is negative or positive. Finally, the term proportional to μ is the chemical potential, which controls the average particle density in the Gibbs state. See Fig. 1. Concerning the manybody interaction, we assume it to be a density-density interaction of the form: with v a finite range, rotationally invariant, potential (we recall that e 1 is the first Euclidean basis vector). The noninteracting Hamiltonian can be rewritten as: where: (2.6) The corresponding energy bands are The size of the bands can be bounded by max k ε + ( k) − min k ε − ( k), which we call the bandwidth. To make sure that the energy bands do not overlap, we assume that t 2 /t 1 < 1/3. For L → ∞, the two bands can touch only at the Fermi points k ± If, instead, m + and m − are both different from zero, then the spectrum ofĤ ( k) is gapped for all k, corresponding to an insulating phase.

Lattice Currents and Linear Reponse Theory
Let n x = ρ=A,B n x,ρ be the total density operator at x. Its time-evolution is given by n x (t) = e iHt n x e −iHt , which satisfies the following lattice continuity equation: with j x, y the bond current: Notice that j x, y = − j y, x . Thus, using that H ( x) = 0 if and only if x = 0, ± 1 , ± 2 , ±( 1 − 2 ), Eq. (2.9) implies: is the lattice derivative along the i direction, and: The operatorsJ i, x are the components along the i directions of the total vectorial current, defined as Note that, given the definitions of 1,2 , the components of the lattice current along the two reference, orthogonal, coordinate directions are: (2.14) We are interested in the transport properties of the Haldane-Hubbard model, in the linear response regime. The Gibbs state of the interacting model is defined as: · β,L = Tr F L · e −βH /Z β,L with Z β,L = Tr F L e −βH the partition function. We define the conductivity matrix via the Kubo formula, for i, j = 1, 2: with J = x∈ L j x the total current operator, X the second quantization of the position operator 1 , and where lim β,L→∞ must be understood as lim β→∞ lim L→∞ , i.e., thermodynamic limit first, and then temperature to zero. Note that formally, in the thermodynamic limit, J = i[H, X ], as it should. Equation (2.15) describes the linear response of the average current at the time t = 0 to an adiabatic external potential of the form e ηt E · X , see e.g. [17] for a formal derivation, and [8,36,39] for a rigorous derivation in a slightly different setting. Remark. The indices i, j labelling the elements of the conductivity matrix (2.15) refer to the two reference, orthogonal, coordinate directions. Sometimes, a similar definition of the Kubo matrix is given, where, instead, the indices i, j label the two lattice coordinate directions 1 , 2 ('adapted basis'). The two definitions are, of course, related in a simple way, via the transformation induced by the change of basis. In particular, the transverse conductivities defined in the orthogonal and in the adapted basis are the same, up to an overall multiplicative factor, equal to | 1 ∧ 2 |. The longitudinal conductivities are, instead, related via a matrix relation that mixes up the diagonal and non-diagonal components of the conductivity matrix. For ease of comparison with experimental papers on graphene, or graphene-like materials, we prefer to use the definition involving the orthogonal reference directions, which we find more natural.
In the absence of interactions, the Kubo conductivity matrix of the Haldane model can be computed explicitly. Suppose that m ω = 0, both for ω = + and for ω = −, and let us choose the chemical potential in the spectral gap. For instance, let μ = −2t 2 cos φα 1 (k ω F ), which 1 There is an issue in defining the position operator on the torus. In order to avoid the problem, we interpret corresponds to choosing the chemical potential in the 'middle of the gap'. Then, it turns out that [23]: The integer ν is the Chern number of the Bloch bundle associated toĤ ( k). The zeros of m ω = W + ω3 √ 3t 2 sin φ, with ω ∈ {+, −}, define the critical curves of the Haldane model, which separate the different topological phases, corresponding to different values of ν. On the curves, the spectrum is gapless: the energy bands intersect with conical intersection, and the system displays a quantization phenomenon of the longitudinal conductivity: while σ 11 = σ 22 = 1 4 at the 'graphene points' m + = m − = 0.

Main Result: Interacting Topological Phases and Phase Transitions
Let us now turn on the many-body interaction, U = 0. In previous works, it was proved that the quantization of the conductivity persists, but only for interactions of strength much smaller than the gap of H 0 . Our main result, summarized in the next theorem, overcomes this limitation.
An illustration of how the interaction deforms the critical lines is shown in Fig.2. The main improvement of the result stated in Theorem 2.1 compared to previous works is that it establishes the quantization of the Hall conductivity for values of the coupling constant U that are much larger than the gap of the bare Hamiltonian: it states that the interaction does not change the value of the interacting Hall conductivity, provided we do not cross the interacting critical curves, which we construct explicitly; this universality of the Hall coefficient holds, in particular, arbitrarily close to the critical curves. On the critical curves the system is massless, i.e., correlations decay algebraically at large distances, and we do not have informations on the transverse conductivity coefficient. However, the critical longitudinal conductivity displays the same quantization phenomenon as the non-interacting one: namely, if W = W R ω (φ), for either ω = + or ω = −, and φ = 0, π, then while σ 11 = σ 22 = 1 4 for (W , φ) = (0, 0), (0, π); see [19] for the proof. We remark that the proof of Theorem 2.1 is constructive: therefore, a patient reader can extract from it an explicit bound on U 0 . Such a bound would certainly be far from optimal; optimizing it would be a non-trivial, interesting, exercise, requiring a computerassisted proof (at least if one is interested in getting a physically significant bound). In any case, conceptually, the only important requirement should be that U is sufficiently small, compared to the bandwidth of H 0 , see the definition after (2.7).
Finally, concerning the model: we expect that the specific choice of the interacting Haldane model is not crucial for the validity of the result. The proof extends straightforwardly to strictly related models, such as the spin-conserving Kane-Mele model. An appropriate adaptation should apply, more generally, to any interacting Hamiltonian of the form H = H 0 +U V, with: (i) V a short-range, spin-independent, interaction, (ii) |U | small compared to the bandwidth, and (iii) H 0 a quadratic Hamiltonian that can become gapless as a parameter is varied: in the gapless case, H 0 has a degenerate, point-like, Fermi surface, around which the dispersion relation has a linear, 'graphene-like', behavior. Note that, as discussed in the introduction, the latter condition is needed to guarantee the irrelevance of the interaction. Even if conceptually non problematic, the extension to such a general class of many-body Hamiltonians would require a discussion of their symmetry properties, in connection with the classification of the possible relevant and marginal effective coupling that can be generated under the multiscale Renormalization Group construction of the Euclidean correlations, cf. with Sect. 4 below.
This goes beyond the scopes of this article: for this reason, we restrict to the specific example of the interacting Haldane model, which is physically the most relevant for applications to 2D topological insulators.

Strategy of the Proof
Let us give an informal summary of the main steps of the proof. For simplicity, we limit ourselves to the generic case W = 0, φ = 0, the special, symmetric, complementary case (W = 0 and/or φ = 0) being treatable analogously. Thanks to the symmetries of the model, see Eqs.(4.7)-(4.13) below, we further restrict ourselves, without loss of generality, to the range of parameters which corresponds to the case m + > |m − |, where m ± are defined in (2.8). Note that, under these conditions, the amplitude of the bare gap is given by |m − |.
We expect the interaction to modify ('renormalize') in a non trivial way both the chemical potential and the width of the gap 2 . In order to compute the interacting gap, we proceed as follows. For the purpose of this discussion, let us denote by H 0 (W , φ, μ) the non-interacting Hamiltonian (2.2), thought of as a function of the parameters (W , φ, μ), at fixed t 1 , t 2 . We rewrite μ in the form μ = −2t 2 cos φ α 1 where the parameter d will be chosen in such a way that m R,− = m − − d has the interpretation of renormalized gap. By using these rewritings, we find: Let us now introduce the reference Hamiltonian H R , thought of as a function of the parameters U , m R,− , φ, defined by Note that H in (2.21) has the same form as H R in (2.22), with the important difference that in passing from H to H R , the parameters d and z have been replaced by the two functions δ(U , m R,− , φ) and ξ(U , m R,− , φ); for the moment, these two functions should be thought of as being arbitrary: they will be conveniently fixed below. Therefore, H R is in general different from the original Hamiltonian H. However, by construction, , and m R,− is a solution of the fixed point equation (2.23) Our construction, described below, will allow us to fix the counterterms ξ(U , m R,− , φ) and δ(U , m R,− , φ) in such a way that they are small, of order O(U ), and that, as anticipated above, m R,− has the interpretation of renormalized gap: in particular, the condition m R,− = 0 implies that the system is massive, that is, correlations decay exponentially at large distances, with decay rate m R,− . Given these definitions, the main steps of the proof are the following.
and we show that |d(U , W , φ)| C|U |(W + sin φ). The equation for the interacting critical curve has the form: (iv) Finally, once we derived explicit estimates on the decay properties of the Euclidean correlations, we infer the identity between the Euclidean and the real-time Kubo conductivity, via [2, Lemma B.1].
The key technical difference with respect to the strategy in [22] is the rewriting of the model in terms of the renormalized reference Hamiltonian H R 0 : this allows us to take into account the renormalization of the gap and of the chemical potential, which characterizes the interacting critical point of the theory.

Lattice Conservation Laws and Universality
In this section, we show how lattice conservation laws can be used to prove the universality of the Euclidean Kubo conductivity, see step (i) above. The main result of this section is summarized in Lemma 3.4. Before getting to this lemma, in Sect. 3.1 we introduce the Euclidean formalism and derive the Ward identities, associated with the lattice continuity equation (2.9), for the Euclidean correlations. In Sects. 3.1.1 and 3.1.2 we differentiate and manipulate the Ward identities, under the assumption that the current-current correlations are sufficiently smooth in momentum space, thus getting some important identities, summarized in Lemma 3.1 and 3.2. Finally, in Sect. 3.2, we prove Lemma 3.4, by combining these identities with the Schwinger-Dyson equation.

Euclidean Formalism and Ward Identities
Given an operator O on F L and t ∈ [0, β), we define the imaginary-time evolution generated by the Hamiltonian H R , Eq. (2.22), as t n on F L , each of which (i) can be written as a polynomial in the time-evolved creation and annihilation operators ψ ± (t, x),ρ = e tH R ψ ± x,ρ e −tH R , (ii) is normal-ordered, and (iii) is either even or odd in ψ ± (t, x),ρ , we define their time-ordered average, or Euclidean correlation function, as: where the (linear) operator T is the fermionic time-ordering, acting on a product of fermionic operators as: x),ρ,s ), and π is a permutation of {1, . . . , n} with signature sgn(π) such that t π (1) . . . t π(n) . If some operators are evaluated at the same time, the ambiguity is solved by normal ordering. We also denote the connected Euclidean correlation function, or cumulant, by Let j μ, x , with μ ∈ {0, 1, 2}, be the three-component operator such that j 0, x := n x , while j i, x , with i ∈ {1, 2}, are the components of the total current along the reference, orthogonal, coordinate directions, see (2.14). Note that j μ, x is the natural current operator, associated both with H and with H R , because i[H, n x ] = i[H R , n x ]. Therefore, its imaginary-time evolution with respect to H R satisfies the analogue of the continuity equation (2.11), cf. with (3.5) below. We define the normalized current-current correlation functions as: for μ i ∈ {0, 1, 2}. We also denote the infinite volume, zero temperature limit of the Euclidean correlations by: where, in the second term, · R ∞ := lim β→∞ lim L→∞ 1 L 2 · R β,L , and the expression [J j , X i ] must be understood as explained in the footnote 1 above. This definition can be obtained via a formal 'Wick rotation' of the time variable, t → −it, starting from the original definition of the Kubo conductivity, (2.15), see, e.g., [17]. A posteriori, we will see that in our context the two definitions coincide, see Sect. 5 below.
The structure correlation functions, and hence the conductivity, is severely constrained by lattice Ward identities. These are nonperturbative implications of lattice continuity equation, which we rewrite here in imaginary time (cf. with (2.11)): where we used the notation div For instance, consider the current-current correlation function 3 , where θ(t) is the Heaviside step function and the correlations in the right side are the timeunordered ones (i.e., they are defined without the action of the time-ordering operator). Using the continuity equation Eq. (3.5): Let us now take the Fourier transform of both sides: integrating by parts w.r.t. x 0 and using (3.7), we find where we used thatJ i,x = j x · G i 2π , with G i , i = 1, 2, the vectors of the dual basis, see definition in Sect. 2.1. More generally, denoting (0, ν 2 , . . . , ν n ) by (0, ν), one has: ] (here, with some abuse of notation, we letĵ μ,(t, p) be the imaginary-time evolution at time t ofĵ μ, p ), and with the understanding that p n = −p 1 − . . . − p n−1 . Even more generally, the identity remains valid if some of the current operators j ν i ,p i are replaced by other local operatorsÔ i,p i : in this case, of course, the operators C ν i must be modified accordingly. In the following, we will be interested in replacing one of the current operators either by the staggered densitŷ where n p,ρ is the Fourier transform of n (t, x),ρ := σ ψ + (t, x),ρ,σ ψ − (t, x),ρ,σ , or by the quartic interaction potential As we shall see below, the combination of the identity (3.9) together with the regularity of the correlation functions has remarkable implications on the structure of the correlations.

Consequences of the Ward Identities for C 1 Correlations
Here we start by discussing the consequences of the Ward identities for continuously differentiable correlations.
We differentiate both sides w.r.t. p i , and take the limit p → 0, thus getting (recall that i · G j = 2πδ i, j ): where we also used that j 0, x = n x . Taking the limit β, L → ∞ and the derivative with respect to (3.4), and the expression [J j , X i ] must be understood as explained in the footnote 1 above. In conclusion, and, if we plug this identity in (3.4), noting that [J j , X i ] R ∞ is even under the exchange i ← → j, we obtain the desired identity.

Consequences of the Ward Identities for C 3 Correlations
Next, we discuss some other implications of the Ward identities for C 3 three-point correlations of the current operator (twice) with either the staggered densityĵ 3,p (see (3.10)), or the interaction potential (see (3.11)), defined as be the new Schwinger terms (recall that C j was defined right after (3.9)). As usual, we denote by K R μ,ν, , S R j, the β, L → ∞ limits of K β,L;R μ,ν, (· · · ), S β,L;R j, B ε (0), and that they are of class C 3 in this domain. Then: In particular, the left side of Eq. (3.18) vanishes as p 0 → 0.
Proof Taking the β, L → ∞ limit of the Ward Identity (3.9) with ν = (0, ), we find Similarly, choosing ν = ( j, ) 20) and, exchanging the roles of p and q, we also get Combining (3.19) with (3.21), we find We now derive w.r.t. p i , q j , and then set p = −q = ( p 0 , 0), thus finding 4 : Finally, notice that ∂ p 1, j S R i, (− p 0 , 0), ( p 0 , 0) is constant in p 0 (recall the definition of Schwinger term, Eq. (3.16), and of C j , Eq. (3.9)). Therefore, after differentiation in p 0 , the final claim follows. 4 In (3.23), we denote by p 1 and p 2 the first and second arguments of K R 0,0, , as well as of S R i, ; correspondingly, we denote by ∂ ∂ p 1,i and ∂ ∂ p 2,i the derivatives with respect to the i-th components of the first and second arguments thereof.

Universality of the Euclidean Conductivity Matrix
Here we prove the universality of the Euclidean conductivity matrix, defined in Eq. (3.4). We restrict to the range of parameters (2.20), as discussed at the beginning of Sect. 2.3.1. In terms of the renormalized parameters, we restate (2.20) as where m R,+ := m R,− + 6 √ 3 t 2 sin φ . (3.25) A key ingredient in the proof is the following regularity result for the correlation functions. The proof of this proposition is postponed to the next section. Its content, combined with the (consequences of the) Ward identities discussed above, immediately implies the universality of the Euclidean conductivity matrix. Proof (Assuming the validity of Proposition 3.3). Thanks to Proposition 3.3, we know that the correlation functions K R μ,ν (p), K R μ,ν, (p, q), and the Schwinger terms S R j (p), S R j, (p, q), with ∈ {0, 3, V }, are C 2 in p, q ∈ B ε (0), for |U | < U 0 . Therefore, we can apply Lemma 3.1 and Lemma 3.2. Using Lemma 3.1, we rewrite the Euclidean conductivity matrix as: Then, we rewrite K R i, j in terms of the non-interacting current-current correlation associated with H R 0 , via the following interpolation formula: is the correlation associated with the (β, L → ∞ limit of the) Gibbs measure with Hamiltonian cf. with Eq. (2.22). Computing the derivative in U : We now take the derivative w.r.t. p 0 and take p 0 → 0. Using Lemma 3.2, we immediately get: is the non-interacting Euclidean conductivity associated with the quadratic Hamiltonian H R 0 at m R,− , which is assumed to be different from zero). The final claim, Eq.

Proof of Proposition 3.3
The proof of Proposition 3.3 is a rather standard application of RG methods for fermions (see, e.g., [11,16,18,32] for reviews). A similar analysis for interacting graphene, which corresponds to the case t 2 = W = 0, has been discussed in [20,21], which we refer to for further details. See also [19], where an application to the Haldane-Hubbard model was discussed. The RG construction of the ground-state correlation functions, uniformly in the gap, is ultimately made possible by the fact that the many-body interaction, in the critical, massless, case, is irrelevant in the RG sense. The only qualitative effect of the interaction, with respect to the non-interacting theory, is a finite renormalization of the gap, of the chemical potential, of the Fermi velocity and of the wave function renormalization.
We recall once more that we restrict the discussion to the range of parameters (3.24). Moreover, we assume that W is not too large, W M 0 , for a pre-fixed constant M 0 , the case of large W being substantially simpler, and left to the reader (for large W , the system is massive and is in a trivial, non-topological, insulating phase, as it follows from the proof of [22]). Finally, for simplicity, we set t 1 = 1, that is, we set the scale of the bandwidth equal to one.
Proof The starting point is the well-known representation of the Euclidean correlation in terms of Grassmann integrals (see, for instance, [20,22]). The generating functional of the correlations is denoted by W( f , A), with f an external Grassmann field coupled to the fermionic fields, and A a (five-component) external complex field conjugated to the lattice currents and the quartic interaction. We have: where: ± x,s , with x = (x 0 , x) ∈ [0, β) × L and s ∈ {↑, ↓}, is a two-component Grassmann spinor, whose components will be denoted by ± x,ρ,s , with ρ = A, B; P(d ) is the fermionic Gaussian integration with propagator where Z L = Z/LZ and, letting and recalling that we set t 1 = 1, with the understanding that, at contact, g(x, x) should be interpreted as lim ε→0 where n x,ρ = s=↑,↓ + x,ρ,s − x,ρ,s is the Grassmann counterpart of the density operator, and n x = ρ=A,B n x,ρ ; finally, whereĴ p,μ = β 0 dx 0 x∈ L e −ip·x J x,μ and: J x,0 = n x is the Grassmann counterpart of the density; J x,1 , J x,2 are the Grassmann counterparts of the two components of the lattice current, and 3 = n x,A − n x,B is the Grassmann counterpart of the staggered density; J x,4 is the Grassmann counterpart of the quartic interaction, The derivatives of the generating functional computed at zero external fields equal the Euclidean correlation functions, cf. with, e.g., [19,Eq. (27), (28)]. Needless to say, the Euclidean correlations satisfy non trivial Ward Identities, following from the lattice continuity equation. For an example, cf. with [19,Eq. (19), (20)].
In order to compute the generating functional W( f , A) in Eq. (4.1), we use an expansion in U , which is convergent uniformly in the volume and temperature, and uniformly close to (and even on) the critical lines m R,± = 0. Note that, in the parameter range (3.24) the propagatorĝ(k) is singular only when m R,− = 0, in which case the singularity is located at k . Due to this singularity, the Grassmann integral has, a priori, an infrared problem, which we resolve by a multi-scale re-summation of the corresponding singularities.
The multi-scale computation of the generating function proceeds as follows. First of all, we distinguish the ultraviolet modes, corresponding to large values of the Matsubara frequency, from the infrared ones, by introducing two compactly supported cut-off functions, χ ± (k), supported in the vicinity of the Fermi points k ± F = (0, k ± F ); more precisely, we let where χ 0 is a smooth characteristic function of the ball of radius a 0 , with a 0 equal to, say, 1/3) and by letting χ uv (k) = 1 − ω=± χ ω (k). We correspondingly split the propagator in its ultraviolet and infrared components: where g (1) (x, y) and g ( 0) ω (x, y) are defined in a way similar to Eq.(4.2), withĝ(k) replaced by χ uv (k)ĝ(k) and by χ 0 (k)ĝ(k + k ω F ), respectively. We then split the Grassmann field as a sum of two independent fields, with propagators g (1) and g ( 0) : x,s,ω and we rewrite the Grassmann Gaussian integration as the product of two independent Gaussians: P(d ) = P(d ( 0) )P(d (1) ). By construction, the integration of the 'ultraviolet' field (1) does not have any infrared singularity and, therefore, can be performed in a straightforward manner, thus allowing us to rewrite the generating function W( f , A) as the logarithm of where V (0) and B (0) are, respectively, the effective potential and the effective source where, denoting the Pauli matrices by σ 1 , σ 2 , σ 3 , we defined that is, T is the spatial rotation by 2π/3 in the counter-clockwise direction.
(2) Complex conjugation: where c ∈ C is a generic constant appearing in P(d ) or in V (ψ) and c * is its complex conjugate.
(3) Horizontal reflections: These symmetries have nonperturbative consequences on the structure of the effective interaction action V (0) . At fixed W , φ, the theory is invariant under the transformations (1), (2)+(4), and (2)+(5). In particular, these transformations leave the quadratic part (4.14) of the effective potential V (0) ( ) invariant (in (4.14), dk 2π |B| is a shorthand for the Riemann sum (β L 2 ) −1 . This means that: The values ofŴ for two real constants ξ ω,0 and δ ω,0 . Let us now discuss the structure of the derivative of the kernel of the quadratic terms. By taking the derivative of Eq. (4.15) w.r.t. k and then setting k = 0, we get: where R v (resp. P) is the diagonal matrix with diagonal elements (1, 1, −1) (resp. (1, −1, −1)). 4.18) implies that: where u ω , z 1,ω , z 2,ω are real constants. The integration of ( 0) ω is performed iteratively. One rewrites ω , supported for quasi-momenta k such that a 0 2 h−1 |k | a 0 2 h+1 , will be defined inductively. We consider two different regimes. The first corresponds to scales h h * 1 , with h * 1 := min{0, log 2 m R,+ }, (4.20) and the rest to scales h * 1 h h * 2 with h * 2 := min{0, log 2 |m R,− | } (recall that we are focusing on the case that m R,+ > |m R,− |.). We describe the iteration in an inductive way. Assume that the fields (0) , (−1) , . . . , (h+1) , h h * 1 , have been integrated out and that after their integration the generating function has the following structure, analogous to the one at scale 0: where V (h) and B (h) are, respectively, the effective potential and source terms, satisfying the conditions that is the Grassmann Gaussian integration with propagator (diagonal in the s and ω indices) where, letting (4.24) and (4.26) and the understanding that (−1) ρ−1 is equal to +1, if ρ = A, and equal to −1, if ρ = B. The quantities Z ρ,ω,h and v ω,h are real, and they have, respectively, the meaning of wave function renormalizations and of effective velocities. Note that r ω ( k ) and s ω ( k ) are both of order O(| k | 2 ), while the mass satisfies (again, recall that m R,+ = m R,− + 6 √ 3t 2 sin φ): By definition, the representation above is valid at the initial step, h = 0. In order to inductively prove its validity at the generic step, let us discuss how to pass from scale h to scale h −1, that is, how to integrate out the field (h) , and how to re-express the resulting effective theory in the form (4.21), with h replaced by h − 1. Before integrating the (h) field out, we split V (h) and B (h) into their local and irrelevant parts (here, for simplicity, we spell out the definitions only in the f = 0 case): we let: By the symmetries of the model,  ( , 0, A). Notice that their structure is constrained by the Ward Identities. E.g., using [19,Eq. (20)], one finds that γ 0,ω,h = − 2 ρ=1 (Z ρ,ω,h + z ρ,ω,h )n ρ (where n ρ = (1 + (−1) ρ−1 σ 3 )/2), γ 1,ω,h = −(v ω,h + u ω,h )σ 2 , and γ 2,ω,h = −ω(v ω,h + u ω,h )σ 1 . However, in the following, we will neither need these identities, nor to identify any special structure of γ μ,ω,h , with μ = 3, 4.
Once the effective potential and source have been split into local and irrelevant parts, we combine the part of LV (h) in the second line of (4.27) with the Gaussian integration P(d ( h) ), thus defining a dressed measureP(d ( h) ) whose propagatorg , y), with the only difference that the functions a ρ,ω,h , b ω,h in (4.25)-(4.26) are replaced bỹ Now, by rewriting the support function χ h (k ) in the definition ofg y) is defined exactly as in (4.25), (4.26), with h replaced by h − 1, and Z ρ,ω,h−1 , v ω,h−1 defined by the flow equations: We are now ready to integrate the fields on scale h. We define: whereP(d (h) ) is the Gaussian integration with propagatorg , we obtain the same expression as (4.21), with h replaced by h − 1. This concludes the proof of the inductive step, corresponding to the integration of the fields on scale h, with h h * 1 . By construction, the running coupling constants τ h = (ξ ω,h , δ ω,h , Z A,ω,h , Z B,ω,h , v ω,h ) ω∈{±} verify the following recursive equations: for suitable functions β ω,h , known as the (components of the) beta function. Note that the initial data ξ ω,0 , δ ω,0 , Z ρ,ω,0 , v ω,0 are analytically close to ξ, δ, 1, 3 2 , respectively; they are not exactly independent of the indices ρ, ω, due to the effect of the ultraviolet integration. However, for small values of m R,+ , the difference between the initial data, for different values of the indices, differ at most by O(U m R,+ ) (note that m R,+ = O(|m R,+ | + sin φ)). As we shall see below, the running coupling constants remain analytically close to their initial data, for all h 0. Similarly, the vertex functions satisfy recursive equations driven by the running coupling constants themselves: (U , τ h , . . . , τ 0 ) , whose solution remains analytically close to the corresponding initial data, for all h 0. From the structure and properties of the effective propagator on scale h, see (4.25) and following lines, one recognizes that the effective theory at scale h is a lattice regularization of a theory of relativistic fermions with masses m R,± . As anticipated above, Z ρ,ω,h and v ω,h remain analytically close to their initial data 1, 3 2 , for all h 0: therefore, it is straightforward to check that the single scale propagator satisfies Moreover, the single-scale propagator admits the decomposition: ω (x, y) by setting m R,ω = 0, and where the remainder term g (h) ω,r satisfies the same bound as g (h) ω times an extra factor m R,ω 2 −h , which is small, for all scales larger than h * 1 . Due to the fact that m R,+ |m R,− |, once we reach the scale h = h * 1 , the infrared propagator of the field corresponding to ω = + satisfies the following bound: that is, it admits the same qualitative bound as the corresponding single scale propagator on scale h = h * 1 . For this reason, it can be integrated in a single step, without any further need for a multiscale analysis. We do so and, after its integration, we are left with an effective theory on scales h h * 1 , depending only on , which we integrate in a multiscale fashion, similar to the one described above, until the scale h = h * 2 is reached. At that point, the infrared propagator g ( h * 2 ) − satisfies a bound similar to (4.32), with h * 1 replaced by h * 2 , and the corresponding field can be integrated in a single step. The outcome of the final integration is the desired generating function. The iterative integration procedure described above provides an explicit algorithm for computing the kernels of the effective potential and sources. In particular, they can be represented as sums of Gallavotti-Nicolò trees, identical to those of [20, Section 3], modulo the following minor differences. The endpoints v on scale h v = +1 are associated either with F  ( ( 0) , f , A), or with one of the terms in RV (0) ( ( 0) ) or in RB (0) ( ( 0) , f , A); the endpoints on scale h v 0 are, instead, associated either with F −1) , f , A). The most important novelty of the present construction, as compared with [20], is the presence of the relevant couplings ξ ω,h , δ ω,h , whose flow must be controlled by properly choosing the counterterms ξ and δ, see discussion below. Recall that the flows of ξ +,h and δ +,h stop at scale h * 1 ; for smaller scales, we let ξ +,h = δ +,h = 0, ∀h < h * 1 . Similarly, we let the other running coupling constants with ω = +, that is, Z ρ,+,h and v +,h , be zero for scales smaller than h * 1 . It turns out that the tree expansion is absolutely convergent, provided that U is small enough and the relevant couplings remain small, uniformly in the scale h 0. More precisely, the kernels of the effective potential satisfy the following bound (a similar statement is valid, of course, for the kernels of the effective source). Notation-wise, we let W (h) n (x 1 , . . . , x n ) be the kernel of the effective potential V (h) ( ) associated with the monomial in of order n; of course, W (h) n is non zero only if n is even. The arguments x 1 , . . . , x n are the space-time coordinates of the Grassmann fields; the kernel implicitly depends also on the ρ, ω indices of the external fields, but we do not spell out their dependence explicitly. We also let W (x 1 , . . . , x n )| (here dx is a shorthand for β 0 dx 0 x∈ L ), which is independent of x 1 , due to translational invariance.

Lemma 4.1
There exist positive constants U 0 , θ, C 0 , such that the following is true. Suppose that max ρ,ω,k h {|Z ρ,ω,k − 1|, |v ω,k − 3 2 |, |ξ ω,k |, |δ ω,k |} C|U |. Then, the kernels of the effective potential on scale h − 1 are analytic in U for |U | U 0 /(C + 1), and satisfy the bound The components of the beta function are analytic in U in the same domain, and satisfy: The proof of the lemma goes along the same lines as the proof of [20,Theorem 2], see also the review [18], and will not be repeated here. Two key ingredients in the proof are: the representation of the iterated truncated expectations in terms of the Brydges-Battle-Federbush determinant formula, and the Gram-Hadamard bound. The factors 2 θ h appearing in the right sides of (4.33), (4.34) and (4.35), represent a 'dimensional gain', as compared to a more basic, naive, dimensional bound, proportional to 2 (3−n)h , which is suggested by the fact that the scaling dimension of the contributions to the effective potential with n external fermionic is equal to 3 − n, in the RG jargon (we use the convention that positive/negative scaling dimensions correspond to relevant/irrelevant operators). Such a dimensional gain is due to the RG irrelevance of the quartic interaction (note that 3−n = −1 for n = 4) and to the so-called short-memory property of the Gallavotti-Nicolò trees ("long trees are exponentially suppressed"): all the contributions to the effective potential associated with trees that have at least one endpoint on scale +1 have this additional exponentially decaying factor. The only contributions not having such a gain are those associated with trees without endpoints on scale +1. The key remark is that, since the running coupling constants are all associated with quadratic contributions in the fermionic fields, such contributions are very simple and explicit: they can all be represented as sums of linear Feynman diagrams with two external legs ('chain diagrams'), obtained by contracting in all possible ways the two-legged vertices corresponding to the running coupling constants ξ ω,k , δ ω,k . Therefore, they only contribute to the quadratic part of the effective potential, and they lead to the first term in the right side of (4.33). Note also that such diagrams do not contribute to the beta function: in fact, the beta function at scale h is obtained by taking the 'local part' of W (h) 2 , which is equal to the value of the Fourier transform W (2) 2 at k = 0. If we compute the chain diagrams at k = 0, we see that the quasi-momenta of all the propagators of the chain diagram are equal to zero; therefore, the value of the diagram is zero, too, due to the compact support properties of the single-scale propagator.
The idea, now, is to use the bound on the beta function to inductively prove the assumption on the running coupling constants, or, more precisely, the following improved version of the inductive assumption: for a suitable C > 0 (recall that, by definition, ξ +,h = δ +,h = Z ρ,+,h = v +,h = 0, ∀h < h * 1 ). Note that the bound on the beta function is already enough to prove the assumption for Z ρ,ω,h and v ω,h . The subtle point is to control the flow of ξ ω,h , δ ω,h , provided the initial data ξ, δ are properly chosen. This is the content of the next lemma.
We want to show that the map τ → T(τ ) admits a unique fixed point in the ball B 0 = {τ ∈ M : τ θ C|U |}, for a suitable C > 0. In order to prove this, we show that, if τ , τ ∈ B 0 , for a suitable C. Once (4.51) is proved, the existence of a unique fixed point in B 0 follows via the Banach fixed point theorem, and we are done: such a fixed point defines the initial data ξ −,0 , δ −,0 generating a solution to the flow equation satisfying (4.36), as desired. Of course, fixing ξ −,0 , δ −,0 is equivalent (thanks to the analytic implicit function theorem) to fixing ξ, δ: therefore, the existence of such a fixed point proves the statement of the lemma. We are left with proving (4.51). If τ ∈ B 0 , by using the bound (4.35) on the beta function, as well as the assumptions (4.41), (4.50) on the initial data (together with their analogues for d + ,v ω ), it is immediate to check that  τ )). Now, the first term in the right side is bounded by 2 −h |x + | 2C 1 |U |, for all h h * 1 , by (4.41) and the very definition of h * 1 , (4.20). In order to bound the sum 0 τ )), we note that β ξ +,k − β ξ −,k can be expressed as a sum over trees with root on scale k, at least an endpoint on scale +1 (recall the discussion after the statement of Lemma (4.1)) and: either an endpoint corresponding to a difference ξ +,k − ξ −,k , or an endpoint corresponding to δ +,k − δ −,k , or a propagator g (k ) − , with k k. The propagator g (k ) − admits a dimensional bound that is the same as g (k ) ω times a gain factor 2 h * 1 −k ; the differences ξ +,k − ξ −,k and δ +,k − δ −,k are proportional to the same gain factor, due to the assumption that τ ∈ B 0 . All in all, recalling the basic bound on the beta function, (4.35), we find a similar bound, improved by the gain factor 2 h * 1 −k : This, together with the bound on 2 −hx + , implies the desired bound, h 0 and C sufficiently large. Exactly the same argument implies the desired bound for δ +,h − δ −,h .
The proof of the second of (4.51) goes along the same lines, and we only sketch it here. A similar argument, discussed in all details, can be found in [11,Section 4]. Let us focus, for simplicity, on the first component of T(τ ) − T(τ ), which reads: can be represented as a sum over trees with root on scale k, at least an endpoint on scale +1, and: either an endpoint corresponding to a difference ξ ω,k −ξ ω,k , or an endpoint corresponding to δ ω,k −δ ω,k , or a propagator corresponding to the difference between g (k ) ω computed at the values (Z ρ,ω,k , v ω,k ) of the effective parameters and the same propagator computed at (Z ρ,ω,k , v ω,k ), for some k k. The difference between the propagators computed at different values of the effective parameters can be bounded dimensionally in the same way as g (k ) ω , times an additional factor max ρ,ω {|Z ρ,ω,k − Z ρ,ω,k |, |v ω,k − v ω,k |}. Therefore, recalling the basic bound on the beta function, (4.35), we find a similar bound, multiplied by the norm of the difference between the running coupling constants: which implies the desired estimate on the first component of T(τ )−T(τ ). A similar argument is valid for the other components, but we will not belabor the details here.
We now have all the ingredients to prove Proposition 3.3. In fact, in view of Lemma 4.1 and Lemma 4.2, we can fix the counterterms ξ, δ in such a way that the kernels of the effective potential on all scales are analytic in U , uniformly in the scale, and satisfy (4.33). A simple by-product of the proof shows that the kernel W (h) n (x 1 , . . . , x n ) decays faster than any power in the tree distance among the space-time points x 1 , . . . , x n , with a decay length proportional to 2 −h . Analogous claims are valid for the kernels of the effective source term and of the generating function. In particular, recalling that the scale h is always larger or equal than h * 2 , we have that the kernels of the effective potential, which are nothing else but the multi-point correlation functions, are analytic in U and decay faster than any power in the tree distance among their arguments, with a typical decay length of the order 2 h * 2 ∼ |m R,− |. Therefore, for any m R,− = 0, the Fourier transform of any multi-point correlation of local operators is C ∞ in the momenta. In the massless case, the correlations are dimensionally bounded like in the graphene case [20,21]: in particular, the two-point density-density, or current-current correlations decay like |x − y| −4 at large Euclidean space-time separation. For further details about the construction and estimate of the correlation functions, the reader is referred to, e.g., [16,21]. This concludes the proof of Proposition 3.3.

Proof of Theorem 2.1
In order to conclude the proof of Theorem 2.1, we need to prove that: there exists a choice of m R,-for which the Euclidean correlations of the reference model with Hamiltonian H R , see (2.22), coincide with those of the original Hamiltonian H; the Euclidean Kubo conductivity coincides with the real-time one. Cf. with the last two items, (iii) and (iv), of the list after (2.23). We also need to prove the regularity and symmetry properties of the critical curves, stated in Theorem 2.1.
Let us start with discussing item (iii), as well as the C 1 regularity of the critical curves. In order to prove the equivalence of H and H R , it is enough to fix the counterterms as discussed in the previous section, and choose m R,− to be the solution of (2.23). Let us then show that (2.23) can be inverted in the form m R,− = m R,− (U , W , φ), with m R,− (U , W , φ) analytic in U and C 1 in W , φ. We want to appeal to the analytic implicit function theorem. For this purpose, we need to estimate the derivative of δ(U , m R,-, φ) w.r.t. m R,− . Recall that δ −,0 = δ −,0 (U , m R,− , φ) satisfies the second of (4.38), and that δ(U , m R,− , φ) and δ −,0 (U , m R,− , φ) are analytically close (they differ only because of the effect of the ultraviolet integration). Therefore, where β δ −,k (U , τ ) accounts for the difference between δ and δ −,0 due to the ultraviolet integration. Differentiating both sides with respect to the mass, we find: which should be looked at as (a component of) a fixed point equation for the derivatives of the running coupling constants, analogous to the ones solved in the proof of Lemma (4.2). When acting on the beta function, the derivative with respect to m R,− can either act on a propagator g (k ) ω , or on a running coupling constant. When acting on a propagator, it replaces g (k ) ω by ∂ g (k ) ω ∂m R,− , which is bounded dimensionally in the same way as g (k ) ω , times an extra factor proportional to 2 −k . On the other hand, the action of the derivative on a running coupling constant should be bounded inductively, in the same spirit as the proof of Lemma 4.2. All in all, recalling also the basic bound on the beta function, (4.35), we get The last estimate is optimal for small φ. For larger values of φ, one can also take advantage of the symmetry under exchange φ → π −φ (the 'magnetic reflections', see (4.13)) to conclude that the derivative of δ with respect to φ vanishes continuously as φ → (π/2) − . Moreover, by the symmetry properties of the model, δ(U , 0, 0) = 0. Therefore, |δ(U , m R,− , φ)| 2C 2 |U |(|m R,− | + sin φ). Using these bounds and the implicit function theorem, we see that (2.23) can be inverted in the form (2.24), with |d(U , W , φ)| C|U |(W + sin φ) for some constant C. The equation for the critical curve in the parameter range we are considering is simply m R,-= 0, that is W = 3 √ 3t 2 sin φ + δ(U , 0, φ), which is C 1 in φ and, thanks to the symmetries of the problem, it satisfies the properties stated in Theorem 2.1.
We are left with discussing item (iv), that is, the equivalence between the Euclidean and real-time Kubo conductivities. Given our bounds on the Euclidean correlations, the equivalence follows from result discussed in previous papers. In fact, our bounds imply that the current-current correlations, at large space-time separations, decay either faster-thanany-power decay, if m R,− = 0, or like |x − y| −4 , otherwise: therefore, we can repeat step by step the proof of [22,Theorem 3.1], as the reader can easily check. For a slightly modified and simplified proof, see also [2,Appendix B] and [35,Section 5].
This concludes the proof of Theorem 2.1.

Concluding Remarks
In conclusion, the universality of the Hall conductivity (i.e., its independence from the interaction strength) can be seen as a consequence of lattice conservation laws, combined with the regularity properties of the correlation functions. The quantization of the interacting Hall conductivity then follows from its quantization in the non-interacting case: however, an important point in the proof is to compare the interacting system and its conductivity with the right reference non-interacting system, that is, the one with the right value of the mass; this is the reason why we introduce a reference non-interacting system with mass equal to the renormalized mass of the interacting system; in order to fix the correct value of the renormalized mass, we need to solve a fixed point equation for it. The same strategy we proposed in the present context can be easily extended to prove that the Hall conductivity is constant against any deformation of the Hamiltonian, even non-translationally invariant ones, provided that the off-diagonal decay of the Euclidean correlations in space and imaginary time is sufficiently fast, in the sense specified by 5 Proposition 3.3. Note that our universality result is valid as soon as the Fourier transform of the current-current-interaction correlations are C 3 in momentum space, which corresponds to a space-time decay faster than (dist.) −6 (a critical analysis of the proof shows that we need even less: C 2+ε with ε > 0 is a sufficient condition for our construction to work; this translates into a space-time decay faster than (dist.) −5 ). This means that we do not require the existence of a spectral gap, in the strong sense of exponential decay of correlations: sufficiently fast polynomial decay is actually enough. It would be nice to provide a realistic example of a gapless model with fast polynomial decay of correlations, exhibiting a non-trivial, universal behavior of the transverse conductivity; or, in alternative, to exclude the possibility that such a model exists.
A problem connected with the one discussed in this paper, but much more challenging, is to prove universality of the conductivity for clean massless models with slow polynomial decay of correlations: by 'slow', here, we mean that Proposition 3.3 cannot be applied. A first example is the Haldane model, considered in this paper, for values of the parameters on the critical line. In this case, as already recalled after the statement of Theorem (2.1), one can prove the universality of the longitudinal conductivity [19]: the proof, which generalizes the one in [21], uses lattice Ward Identities, combined with the symmetry properties of the current-current correlation functions. It would be very interesting to establish the universality, or the violation thereof, of the transverse conductivity on the critical line.
Another context, where the issue of the universality of the conductivity naturally arises, is the case of bulk massive systems in non-trivial domains with, say, Dirichlet conditions imposed at the boundary. In such a setting, usually, massless edge states appear, and the edge system is characterized by correlations with slow polynomial decay. Nevertheless, universality holds as a consequence of a more subtle mechanism, which relies on the nonrenormalization of the edge chiral anomaly. Using these ideas, two of us proved the validity of the bulk-edge correspondence in lattice Hall systems with single-mode chiral edge currents [2], and in the spin-conserving Kane-Mele model [34]. It would be very interesting to generalize these findings to lattice systems with several edge modes, as well as to continuum systems.
Finally, it would be extremely interesting to include disorder effect, even in the regime where the interaction is smaller than the non-interacting gap. Understanding the combined effects of disorder and interactions in the vicinity of the critical lines is a major open problem, even from a theoretical physics perspectives. We do not expect that the phase diagram will remain qualitatively unchanged in their presence: new quantum phases may in general arise in the vicinity of unperturbed critical lines. In this sense, we expect that the stability of the phase diagram, if valid at all, will depend on the specific features of the model under investigation. However, as far as we know, not even the effects of disorder alone are well understood in the vicinity of the critical lines.