Large-$S$ and tensor-network methods for strongly-interacting topological insulators

The study of correlation effects in topological phases of matter can benefit from a multidisciplinary approach that combines techniques drawn from condensed matter, high-energy physics and quantum information science. In this work, we exploit these connections to study the strongly-interacting limit of certain lattice Hubbard models of topological insulators, which map onto four-Fermi quantum field theories with a Wilson-type discretization, and have been recently shown to be at reach of cold-atom quantum simulators based on synthetic spin-orbit coupling. We combine large-S and tensor-network techniques to explore the possible spontaneous symmetry-breaking phases that appear when the interactions of the topological insulators are sufficiently large. In particular, we show that varying the Wilson parameter $r$ of the lattice discretizations leads to a novel Heisenberg-Ising compass model with critical lines that flow with the value of $r$.

I. Introduction Our most accurate description of Nature is based on a fourdimensional quantum field theory (QFT) of fermionic matter coupled to gauge fields: the standard model of particle physics [1]. In this context, current challenges arise in the understanding of effects that cannot be treated perturbatively, such as quark confinement in quantum chromodynamics [2]. To advance our understanding of non-perturbative phenomena, quantum field theorists have introduced other simplified models that share some important aspects with the standard model but, at the same time, avoid the intricacies of non-Abelian gauge theories. Such toy QFTs, which are typically defined in reduced dimensionalities, have played an important role in elucidating phenomena such as asymptotic freedom, dynamical mass generation, chiral symmetry breaking, or the role of topological solutions and instantons. Paradigmatic examples of such toy QFTs are the two-dimensional Thirring [3] and Gross-Neveu [4] models, which describe self-interacting Dirac fields, and the two-dimensional nonlinear sigma model [5], which consists of scalar fields coupled through a non-linear constraint. These models serve to develop and test tools, such as bosonization [6] and large-N expansions [7], the predictions of which can be benchmarked with efficient numerical methods for such low-dimensional QFTs. Nonetheless, our most accurate experiments are consistent with a four-dimensional spacetime such that, in a strict sense, the specific predictions of these toy models are not solving a specific real problems in high-energy physics that can be falsified experimentally. Within high-energy physics, these toy QFTs qualify instead as theoretical laboratories.
Remarkably, during the last decades, we have witnessed a change of status for such low-dimensional QFTs. Rather than looking at high energies and small length scales, one may instead focus on non-relativistic condensed-matter systems which, at long wavelengths and small energies [8], display certain universal behaviour determined by emergent relativistic QFTs. Such effective descriptions [9] provide a more flexible framework in comparison to the standard model, as realisation properties such as the effective dimensionality or the emergent symmetries are not fixed a priori, but depend instead on the family of materials at hand. Some characteristic examples where relativistic Dirac fields emerge include graphene [10], Weyl semimetals [11], or topological insulators and superconductors [12]. For the emergence of scalar fields subjected to non-linear constraints, quantum magnets play a prominent role [13].
A further step in this direction is provided by the so-called quantum simulators (QSs) [14]: well-isolated quantum manybody systems with unparalleled levels of control down to the single-particle level that can directly mimic a specific target model [15]. In the context of low-dimensional QFTs, QSs can be tailored such that one has full control of the microscopic arXiv:2203.03434v2 [cond-mat.str-el] 8 Mar 2022 parameters and, moreover, can access the continuum limit in a controlled fashion. In order to do so, QSs of QFTs [16][17][18][19][20][21] typically follow the approach of lattice field theories [22] in their Hamiltonian formulation [23]. Rather than reducing the lattice spacing to recover the continuum limit, one may tune the microscopic couplings of these QSs to approach a critical point where the correlation length is much larger than that spacing, and the continuum description sets in. In recent years, we have seen very promising experimental steps in this direction for Dirac fermions [24,25] and Dirac QFTs [26][27][28][29][30] and gauge theories [31][32][33][34][35][36] in low dimensions.
In this work, we explore the strong-coupling limit of four-Fermi models, namely QFTs of self-interacting Dirac fermions. As discussed below, these QFTs are inspired by the Thirring and Gross-Neveu models, the origin of which can be traced back to seminal contributions of E. Fermi [37,38] and Y. Nambu and G. Jona-Lasinio [39,40]. In particular, we explore specific lattice discretizations based on the so-called Wilson fermions [41], which make direct connections of this four-Fermi QFTs with the aforementioned topological insulators [12,[42][43][44][45][46][47][48], allowing to study the effect of electronelectron interactions. We note that recent advances in coldatom QSs based on schemes of synthetic spin-orbit coupling in atomic gases with negligible interactions [29,30] connect directly to these Wilson-regularised lattice field theories and, furthermore, motivate a careful study of the regime of strong interactions. This is not only relevant from the perspective of QFTs, where one can find novel strongly-coupled fixed points that can only be characterised non-perturbatively [49,50], but also from the perspective of strongly-correlated effects in topological phases of matter, a topic that has received significant attention in recent years [51][52][53][54].
As we have discussed in a series of recent works, the native Hubbard interactions [55,56] of cold-atom QSs of spinorbit coupling in two [57][58][59] and three [60,61] dimensions can be understood as the single-flavour-limit of Four-Fermi QFTs with Lorentz-invariant self-interactions, and regularised via a Wilson-type discretization. Such a discretization introduces the Wilson parameter r ∈ (0, 1], which is customarily set to unity r = 1 in most lattice studies. As briefly discussed in [60,61], setting r < 1 has no important effect in the absence of interactions, as one can simply rescale the axes of the phase diagram in a simple manner to maintain the same layout: topological insulators are separated from trivial band insulators by critical lines in parameter space. The situation is not so clear as one switches on the interactions. Here, as a consequence of spontaneous symmetry breaking, one can find phases with long-range order corresponding to different fermion condensates in the context of relativistic QFTs, as discussed in detail for r = 1 [58][59][60][61]. In this work, we explore the nature of these fermion condensates as the Wilson parameter is allowed to take values in r < 1. To identify the possible condensates and chart the phase diagram, we explore the strong-coupling limit by deriving an effective Heisenberg-type compass model with directional spin-spin interactions. Using the path-integral representation of the partition function, we derive a version of the aforementioned non-linear sigma model with a discrete Z 2 symmetry, a constrained QFT amenable to a large-S expansion in the limit where the effective spin S 1. We then benchmark these predictions for the two-and three-dimensional lattice field theories with numerical simulations based on tensor networks.
Besides the fundamental interest in understanding the role of the Wilson parameter in non-perturbative phenomena of four-Fermi models, we note that the specific cold-atom proposals based on synthetic spin-orbit coupling [60,61] lead to an effective Wilson parameter R that is controlled by the ratio of spin-preserving and spin-flipping tunnelings, each of which can be independently controlled by the lasers that form a socalled optical Raman lattice [29,30]. Accordingly, reaching the regime of r = 1 would require additional fine tuning, as the generic QS would instead lead to r = 1. Regarding the possible experimental realisation of the four-fermi-Wilson model with cold atoms, it is thus an interesting and useful question to understand the effects of 0 < r < 1.

B. Constrained quantum field theories
Let us start by discussing the nature of the constraints in representative QFTs, which will allow us to frame the results of our work appropriately. A well-known QFT where an effective constraint arises is the O(N) model, which describes a real scalar field Φ Φ Φ(x) = (φ 1 (x), · · · , φ N (x)) t with N flavours. In the Hamiltonian formulation, and in the absence of interactions, the free fields evolve under a Klein-Gordon Hamiltonian where we use natural unitsh = c = 1, and Einstein's convention of repeated-index summation. Here, the fields and conjugate momenta Π Π Π(t, , and j ∈ {1, · · · , d} labels the spatial coordinates of a D = (d + 1)-dimensional Minkowski spacetime with metric η = diag(1, −1, . . . , −1). This QFT describes N uncoupled scalar bosons, and is invariant under a continuous internal symmetry is an arbitrary rotation. In order to couple the different flavours, one introduces a quartic self-interaction that respects this internal symmetry where λ 0 is the bare coupling strength. Here, we have introduced Φ 0 as the vacuum expectation value attained by one of the scalar-field flavours, e.g. Φ f (x) = δ f,1 Φ 0 , which corresponds to the spontaneous symmetry breaking (SSB) of the continuous O(N) symmetry in the classical limit, such that O(N) → O(N − 1). This corresponds to the meson sector of the linear sigma model [62], which describes the coupling of (N − 2)(N + 1)/2 pions π π π to an additional heavy scalar corresponding to the symmetry-breaking and Goldstone components respectively. Instead of expanding around the SSB groundstate, one may focus on the strong-coupling limit λ 0 → ∞, where the ground-state minimises the interaction energy (2) by imposing a non-linear constraint on the fields which are thus forced to take values on the unit sphere S N−1 . The Klein-Gordon field theory (1) subjected to this constraint (3) belongs to the family of sigma models [62], which describe particles forced to move on a specific manifold. In this particular case, this constrained model is called the O(N) non-linear sigma model [1]. As noted in the introduction, in two-dimensional spacetimes where the O(N) symmetry cannot be spontaneously broken [63], the O(N) non-linear sigma model shares important features with non-Abelian gauge theories such as asymptotic freedom for N > 2 [5], existence of topologically non-trivial solutions called instantons [64], or large-N methods and dimensional transmutation [7,65].
In this work, we will discuss how similar constraints can appear also in purely fermionic QFTs even in the absence of any continuous internal symmetry. In the most general setting, we consider a spinor field Ψ Ψ Ψ(x) = (ψ 1 (x), · · · , ψ N (x)) t with N flavours evolving under a Dirac Hamiltonian density where we have introduced the gamma matrices {γ µ , γ ν } = 2η µν for spacetime indexes µ, ν ∈ {0, 1, · · · , d}, and the adjoint spinor Ψ Ψ Ψ(x) = Ψ Ψ Ψ † (x)(I N ⊗ γ 0 ). This model describes N uncoupled Dirac fermions, and is invariant under a continuous unitary transformation Ψ Ψ Ψ(x) → u ⊗ I s Ψ Ψ Ψ(x), where u ∈ U(N), and the identity in the spinor components I s depends on the dimensionality of the representation of the gamma matrices.
Paralleling the the discussion around Eq. (2), we can now couple the flavours via a four-Fermi [37][38][39][40] contact interaction where g 2 is the bare coupling strength. Although not directly apparent, as in the bosonic case, we will show below that the strong-coupling limit g 2 → ∞ leads to a constraint similar to the one of Eq. (3), where the σ and π π π fields will be related to particular SSB channels of the above QFT related to fermion condensates. In contrast to the bosonic case, the non-linear constraint appears down to the N = 1 level, as the symmetry being broken is not the U(N) symmetry, but rather a discrete Z 2 symmetry involving the spinor degrees of freedom. In the following section, we will introduce a particular lattice discretization, which plays an important role in determining the specific Z 2 SSB, and its connection with topological insulators.

C. Four-Fermi interactions in topological insulators
In this section, we describe in more detail the Wilson regularisation [41] of the above fermionic QFT (4)-(5), and how it yields a playground to explore interactions in topological insulators. We consider the Hamiltonian lattice formulation [23] obtained by discretising the spatial coordinates x x x ∈ Λ d , focusing on the d = 1, 2 cases n j a j e e e j : n j ∈ Z N j , where {a j } are the lattice spacings along the {e e e j } unit vectors, and N j are the corresponding number of lattice sites along each axis (see Fig. 1). Let us also note that for d = 1, 2 spatial dimensions, one can use the following irreducible representations of the gamma matrices which are proportional to the Pauli matrices. Discretising the spatial derivatives appearing in Eq. (4) using central differences leads to the so-called naive fermions [66], the continuum limit of which contains N D = 2 d Dirac fermions due to fermion doubling [67,68]. We follow Wilson's prescription [41], which introduces additional terms that are responsible for giving different masses to each of these Dirac fermion species: where we have introduced the bare mass m, and the aforementioned dimensionless Wilson parameter r. In Fig. 1, we present a schematic diagram of this Wilsonian discretization by means of tunnelling processes and density-density interactions with strengths obtained after rescaling the fields in terms of dimensionless creation-annihilation operators. In lattice field theories (LFTs), one typically works directly with the Euclidean action associated to the above Hamiltonian by also discretising the Wick-rotated temporal coordinate x = ∑ d µ=0 a µ n µ e e e µ , such that recovering the time-continuum limit requires temporal anisotropies that permit the limit a 0 → 0. In case one is not interested in making contact with the Hamiltonian formulation, it is possible to focus directly on the isotropic regime a µ = a, and consider |r| ≤ 1, as imposed by the reflection positivity of the Euclidean action for g 2 = 0 [69]. A standard choice in the literature is to set r = 1, such that one recovers a single massless Dirac fermion in the limit of m → 0 and at long wavelengths a → 0, while the remaining doublers acquire a very large mass proportional to 1/a, and thus lie at the UV cutoff of the regularised QFT. The choice r = 1 brings the technical advantage that tunnelling terms are proportional to projection operators P j± ≡ 1 2 (1 ± iγ j ) with P 2 j± = P j± , P j± P j∓ = 0. The goal of the present work is to explore regimes with 0 < r < 1, and make connections with effective constrained QFTs in the strong-coupling limit. Likewise, we will also explore anisotropic lattice constants a µ = a ν . We remark that isotropy is not required a priori, since the continuum limit yields a QFT invariant under the full Lorentz group SO(1, d) even when the anisotropic lattice formulation breaks translational, rotational and Lorentz symmetries explicitly. In fact, temporal [70] and spatial anisotropies [71] can actually increase the accuracy of lattice computations. To understand the effect of non-unity Wilson parameters and anisotropic lattice constants, we can start by focusing on the non-interacting limit g 2 = 0. In this case, it is straightforward to compute the half-filled groundstates |ε(k k k) corresponding to the Dirac vacua [58,60,61], where we have introduced the quasi-momentum within the first Brillouin zone k k k ∈ BZ = × j (−π/a j , π/a j ]. Associated to this band structure, one finds a Berry connection A (k k k) = ε(k k k)| i∇ ∇ ∇ k k k |ε(k k k) [72,73] , which characterises the principal fibre bundle associated to the occupied energy band [74]. Such fibre bundles can be characterised by topological invariants which depend on dimensionality.
For d = 1 spatial dimensions, Zak's phase [75] allows to define the Wilson loop for a cycle in momentum space where θ (x) is the Heaviside step function. Therefore, one finds a trivial band insulator with W Z = +1 for ma 1 > 0 or ma 1 < −2r. Alternatively, a non-trivial topological insulator W Z = −1 arises when −2r < ma 1 < 0 and the number of flavours N is odd, which actually lies in the symmetry class BDI [76,77]. In comparison to the previous results of [58], which focused on the standard choice r = 1, we observe that the structure of the non-interacting phase diagram is completely equivalent if one simply rescales the dimensionless mass with the Wilson parameter ma 1 → ma 1 /r.
For d = 2, the Berry curvature B(k k k) = ∇ ∇ ∇ k k k ∧ A (k k k) [72] allows to define the first Chern number where ξ 2 = a 1 /a 2 is an anisotropy ratio, and we have assumed ξ 2 ≤ 1. Here, some comments are in order. In the isotropic limit ξ 2 = 1, and setting r = 1, the groundstate corresponds to quantum anomalous Hall (QAH) phase with N Ch = −N for −2 < ma 1 < 0, and N Ch = +N for −4 < ma 1 < −2, whereas it is a trivial band insulator with N Ch = 0 for ma 1 > 0 or ma 1 < −4. This limit can be readily mapped onto the Qi-Wu-Zhang model [42,43] of the QAH [78], with a central region along the ma 1 axis comprising the QAH phases, surrounded by trivial band insulators at both sides. It is worth noting that both the QAH and the BDI topological insulators can be understood as the bulk of a lower-dimensional version of the domain-wall-fermion construction of lattice field theories [79,80], in which the non-zero topological invariants would give rise to effective Chern-Simons-type terms in the response of the fermions to external gauge fields [43,81]. As discussed in [60,61,82], allowing for spatial anisotropies ξ 2 < 1, and fixing the Wilson parameter to the standard value r = 1, one finds an additional trivial band insulator for −2 < ma 1 < −2ξ 2 separating the two QAH phases with N Ch = ±N. Something completely analogous occurs for spatial anisotropies ξ 2 > 1. From the above expression (10), we see again that the effect of a non-unity Wilson parameter 0 < r < 1 is rather trivial in the non-interacting case, one can simply rescale ma 1 → ma 1 /r, and obtain the same structure and phases as in the limit of r = 1. However, as one switches on interactions g 2 > 0, the situation need not be so simple: there can be SSB processes that lead to long-range-ordered phases different from the above trivial and topological insulators. In the following sections, we will extend our previous studies in [58,60,61], studying the nature of these SSB processes for arbitrary Wilson parameters 0 < r < 1 by exploring the strong-coupling limit, in which the four-Fermi interaction strength is the leading parameter g 2 → ∞. As discussed below, a different constrained QFT controls that limit, which will be exploited to predict the shape of the phase diagram and possible phase transitions. Strong-coupling Heisenberg-Ising spin models: (a) For d = 1, the half-filled chain Λ 1 in the limit of strong interactions can be described by localised spins S = 1/2, here depicted with black arrows. The spins interact with the neighbours via spin -spin couplings S a (x x x)S a (x x x + a 1 e e e 1 ) of strength J 1,a for a ∈ {x, y, z}, here represented by solid lines in a red scale that determines their relative magnitude for a generic Wilson parameter 0 < r < 1, such that |J 1,y | dominates. (b) For the d = 2 half-filled lattice Λ 2 , the effective spins S = 1/2 of the strong-coupling limit are arranged in a rectangular lattice, and the spin-spin couplings S a (x x x)S a (x x x + a j e e e j ) have a directional character J j,a , where j ∈ {1, 2} labels the horizontal and vertical neighbours, leading to a compass-type model. In addition to the horizontal interactions in (a), the spins now interact vertically with strengths J 2,a , here represented by solid lines in a blue scale that determines their relative magnitude for 0 < r < 1, such that |J 2,x | dominates. The anisotropy parameter ξ 2 = a 1 /a 2 controls the directionality of the compass Heisenberg-Ising model, i.e. if the vertical (|J 1,a | < |J 2,a |) or horizontal (|J 1,a | > |J 2,a |) couplings dominate, which occurs for ξ 2 < 1 and ξ 2 > 1, respectively.

A. Ising order and fermion condensates
In order to understand the strong-coupling limit, let us first note that for d = 1, 2 spatial dimensions, the irreducible representations of the gamma matrices (7) imply that the Dirac spinors have two components ψ f (x) = (ψ f,1 (x), ψ f,2 (x)). In the single-flavour limit f = 1 = N, and for a fixed total number of fermions, the four-Fermi term in Eq. (8) can be rewritten as up to an irrelevant shift of the ground-state energy. Accordingly, the strong-coupling limit g 2 → ∞ will give rise to a large energy penalty for configurations in which a pair of fermions occupy the same lattice site. The Dirac vacuum corresponding to the half-filled groundstate will have a single fermion per site n n n, which has the freedom to select one of the two possible spinor configurations |↑ n n n , |↓ n n n . Within this subspace, the operators that fulfil the x x x )] = iδ n n n,n n n ∑ c ε abcŜc (t, x x x) at spatial points x x x = ∑ j n j a j e e e j , x x x = ∑ j n j a j e e e j , becomê Here, S = 1/2, and σ σ σ n n n is an operator acting on the projected Hilbert space H eff = ⊗ n n n span{|↑ n n n , |↓ n n n }, and defined by the tensor product of the identity I 2 on all sites except for n n n, where one applies the vector of Pauli matrices σ σ σ = (σ x , σ y , σ z ). As discussed in [58,60,61] for d = 1, 2, and for unit Wilson parameter r = 1, the lattice model that controls this strong-coupling limit corresponds to an effective spin model, where the spins reside on the sites of the spatial lattice regularisation (6). These spins are subjected to local on-site terms, and interact with each other via nearestneighbour couplings depicted in Fig. 2. The physical mechanism underlying these nearest-neighbour spin-spin interactions is the so-called super-exchange [83,84], and the mostgeneral effective Hamiltonian can be written as follows Here, we have used the label a ∈ {x, y, z} to distinguish the internal spin components from the spatial coordinates j ∈ {1, · · · , d}, and introduced a set of spin-spin couplings J j,a describing the strength of the interactions between neighbouring spins connected by a e e e j link, and coupling their internal spin components via a S a S a interaction (see Fig. 2). In previous works [58,60,61] where we used the standard choice r = 1, the nature of the spin-spin couplings was restricted to be of Ising type. For d = 1, where the coupling strength g 2 is dimensionless, the spin couplings found were which have units of inverse length, such that the effective Hamiltonian (13) with dimensionless spin operators (12) has the correct units of energy. The strong-coupling Hamiltonian for r = 1 thus coincides with a quantum Ising model [22,85] with S y S y ferromagnetic interactions subjected to a transverse field along the internal z axis [58]. This is an exactly-solvable model with a quantum phase transition at |h z | = |J 1,y |/2, marking the onset of SSB of an underlying Z 2 symmetry , ∀x x x ∈ Λ 1 . This symmetry can also be combined with a reflection about the lattice centre, such that It is interesting to note that, for the representation of the gamma matrices in Eq. (7), this Z 2 symmetry corresponds to which is precisely the parity symmetry on Dirac spinors [1]. Accordingly, the ferromagnet with all spins aligned with the internal y axis (FM y ) can be readily identified with a paritybreaking pseudo-scalar π condensate where γ 5 = γ 0 γ 1 is the chiral gamma matrix. Note that, in the quantum Ising model [85], the groundstate always displays a non-zero magnetisation along the direction of the transverse field for any h z = 0. In the language of Dirac spinors, this corresponds to a non-zero value of the so-called scalar σ condensate and the fact that is generically non-zero can be traced back to the explicit breaking of the discrete chiral symmetry ψ f (x) → γ 5 ψ f (x), by the Wilson discretisation (8). The appearance of the scalar condensate is typical of four-Fermi models with dynamical mass generation, such as the Gross-Neveu model [4], whereas the pseudo-scalar one depends on the specific lattice regularisation. A non-zero pseudo-scalar condensate has also been discussed in the context of lattice gauge theories in 3+1 dimensions [86,87], and is known as the Aoki phase. We thus see that the strong-coupling limit captures nicely this condensate even in the N = 1 limit and that, moreover, it provides an analytical expression for the critical line that was proved to be very accurate by comparing with matrix-product-state numerics [58]. In this manuscript, we will explore how this situation changes as the Wilson parameter is modified 0 < r < 1. For d = 2 spatial dimensions, where the coupling strength g 2 has units of length, setting r = 1 [60,61] leads to the following spin-spin couplings which again have the correct units of energy as in Eq. (14). The corresponding spin model (13) is an instance of the socalled 90 o compass model [88,89] with directional spin couplings S y S y (S x S x ) along the e e e 1 (e e e 2 ) spatial axis, and a transverse magnetic field again directed along the internal z axis. In contrast to the quantum Ising chain (14), the compass model is no longer exactly solvable, and presents two different types of phase transition. For h z = 0, which is achieved for a negative bare mass m = −1/a 1 − 1/a 2 , there is a well-studied first-order phase transition at J 1,y = J 2,x [90][91][92]. This critical point separates two different ferromagnets: a FM x characterised by the order parameter | S x (x x x) | > 0 for |J 2,x | > |J 1,y |, achieved for a 1 > a 2 , and a FM y characterised by | S y (x x x) | > 0 for |J 1,y | > |J 2,x | for a 2 > a 1 . Both ferromagnets break the aforementioned Z 2 symmetry (15). We remark that in this d = 2 case, the expression of this symmetry in fermion operators (16) does not correspond to parity, as x → −x is generated by a rotation in the connected component of the Lorentz group SO(1, 2). Rather than breaking parity, a non-zero value of the corresponding fermion π condensates breaks inversion symmetry. We note that taking a continuum long-wavelength limit around the critical lines that separate these ferromagnets from the symmetric paramagnet would lead to a QFT where Lorentz symmetry cannot be recovered when approaching from the condensed phase. Accordingly, these d = 2 FM x , FM y phases were referred in [60,61] as Lorentz-breaking condensates, which contrast the paritybreaking pseudo-scalar condensate (17) of d = 1.
In contrast, the regime with a non-vanishing transverse field h z = 0 has not been studied in so much detail. In the limit of very large spatial anisotropies ξ 2 = a 1 /a 2 → 0 (ξ 2 = a 1 /a 2 → ∞), the compass model (19) reduces to a collection of uncoupled rows (columns), each described by an Ising model in a transverse field with a second-order phase transition at |h z | = |J 1,y |/2 (|h z | = |J 2,x |/2). This critical point separates a paramagnet, which has all spins aligned along the internal z axis, from the aforementioned ferromagnet with | S y (x x x) | > 0 ( | S x (x x x) | > 0) for each row (column). We note that there is an accidental exponentially-large degeneracy in the number of rows (columns) in these large-anisotropy limits. For non-zero anisotropies ξ 2 > 0, yet finite ξ −1 2 = 0, these rows (columns) become coupled, lifting the degeneracy and selecting a unique 2-fold degenerate ferromagnetic ground-state FM x (FM y ) where all columns (rows) get locked to the same spin direction when |h z | < |J 2,x |/ζ and |J 2,x | > |J 1,y | (|h z | < |J 1,y |/ζ and |J 1,y | > |J 2,x |). Here, we have introduced a parameter zeta, which serves to locate the critical point, and will be equal to ζ = 2 in the regime of large spatial anisotropies. For other finite and non-zero anisotropies ξ 2 = a 1 /a 2 , the critical point will change, and one should find that ζ = 2, which can no longer be found exactly, but must be estimated using numerical or analytical approximations [60,61]. So far, we have reviewed known results that apply for Wilson parameter r = 1. Moving away from this limit modifies the super-exchange mechanism, leading to additional spinspin couplings depicted in Fig. 2 (a). These still admit the general form of Eq. (13), but lead to different strengths with respect to those expressed in Eqs. (14)- (19). In particular, for d = 1 we find that the expression of the external field is whereas the spin-spin couplings following the super-exchange mechanism of virtual double occupancies now read One can readily see how the previous ferromagnetic Ising model in a transverse field (14) is recovered for r → 1. In this limit the distinction between ferromagnetic and antiferromagnetic couplings is trivial, as one can invert the sign of the spin-spin couplings J 1,y → −J 1,y by a unitary transforma- . Under this transformation, the SSB ferromagnetic groundstate is transformed into a classical Néel pattern of alternating spins. The discussion of the possible SSB orderings for general Wilson parameter 0 < r < 1 is slightly more involved.
In the limit r → 0, where one recovers the naive-fermion regularisation [66], the spin-spin couplings tend to J 1,x = −J 1,y = −J 1,z , which are unitarily-equivalent to a quantum Heisenberg model with antiferromagnetic couplings [93][94][95][96][97]. The important point is that, if there is an even number of spin-spin couplings with negative sign, these can always be inverted by a spin rotation along the remaining axis with an alternating angle. The specific transformation in this case takes y, z}, and the spin-spin interactions clearly displays the SU(2) symmetry of the Heisenberg model. Additionally, if the external field is non-zero, this maps onto a staggered transverse field under the above transformation, such that where we have introduced the wavevector k k k s = π a 1 e e e 1 . This transformation unveils a continuous U(1) symmetry with respect to rotations along the internal z axis, which is not directly apparent in the original formulation (13). It is interesting to note that rewriting this transformed model in terms of Jordan-Wigner fermions maps the spin chain into a staggeredfermion regularisation [23,98] of the D = (1+1)-dimensional Thirring model [99] for a single fermion flavour, provided that the four-Fermi term has a specific coupling strength. We also note that for vanishing transverse field h z = 0, the Heisenberg chain was exactly solved via the Bethe ansatz [100,101], and does not support long-range order, as shown via the inverse scattering method [102].
Given the clear difference between r = 1 and r = 0 limits, one may expect different groundstates with distinct magnetic orders, i.e. fermion condensates, for Wilson parameters 0 < r < 1. In this regime, we recall that the non-interacting phase diagram comprises regions of non-trivial topological insulators and trivial band insulators separated by topological gap-closing phase transitions (9). For strong interactions, the absolute values of the spin-spin couplings (22) are no longer equal, and the mapping to the antiferromagnetic Heisenberg model no longer holds. Interestingly, one can still find a U(1) symmetry by combining a pair of transformations. First, rotate the spins about the internal z axis in an alternate fash- , which effectively changes the spin-spin couplings to J 1,x = J 1,z = J ⊥ , and J 1,y = J ⊥ ∆, where Next, apply a rotation about the internal ; the spin chain maps onto the XXZ model [97,103,104], also known as a Heisenberg-Ising model, under an additional longitudinal field As advanced previously, for g y = h z /J ⊥ = 0, there is a U(1) symmetry with respect to continuous rotations about the new z axis. In this limit, the XXZ model for S = 1/2 is known to display a Berezinskii-Kosterlitz-Thouless (BKT) phase transition [105,106] at the SU(2)-symmetric point ∆ = 1 [107], separating a critical phase at ∆ < 1 from an Ising SSB phase phase at ∆ > 1. This model has also been exactly solved via the Bethe ansatz [108,109], and the inverse scattering method [102] demonstrates the different decay of spin-spin correlations in the critical and Ising phases, yielding long-range order in the latter. Following these results, we expect that such a BKT transition will not appear in the strong-coupling limit of our model (8), as the effective anisotropy (24) always exceeds unity ∆ > 1 for 0 < r < 1. This favours the Ising long-range ordered phase which, after reversing the previous spin transformations, corresponds to FM y order, i.e. pseudo-scalar condensate in the fermion language (17). There may be, however, other types of transition when the longitudinal field is switched on g y = 0. Despite lacking U(1) symmetry in these cases, the global Z 2 symmetry (15) remains intact in the Hamiltonian (25) for any r > 0 provided that we consider its correspondence in terms of the new rotated spin axes. Accordingly, similar second-order quantum phase transitions as those discussed for the quantum Ising model (14) may still appear, albeit at critical points that flow with the Wilson parameter.
C. Heisenberg-Ising compass models for d = 2 Before checking the validity of the above conjecture, let us discuss the effect of non-unit Wilson parameters on the effective spin model for the d = 2 case depicted in Fig. 2 (b). Instead of the microscopic couplings in Eq. (19), the superexchange for 0 < r < 1 leads to an external field whereas the spin-spin couplings transform into In the limit r → 1, we recover the quantum compass model with the directional spin-spin couplings of Eq. (19), which supports Lorentz-breaking fermion condensates (20) as discussed in [60,61]. In the isotropic a 1 = a 2 and naive-fermion r → 0 limits, we recover a model that is unitarily equivalent to the Heisenberg model on a square lattice with antiferromagnetic couplings J j,a → J = 1/g 2 and staggered field h z → h z (−1) n 1 +n 2 . This requires a slightly more involved transformation in comparison to d = 1 (23). For odd rows, the spins must be trans- , ∀x x x = (n 1 a 1 , (2n 2 − 1)a 2 ), whereas for even rows they trans- where the corresponding staggering wave-vector now reads k k k s = π a 1 e e e 1 + π a 2 e e e 2 . For vanishing transverse field h z = 0, one finds that the strong-coupling limit of the naive-fermion r = 0 isotropic limit a 1 = a 2 corresponds exactly to the two-dimensional antiferromagnetic Heisenberg model on a square lattice. This model is no longer solvable via Bethe ansatz [100], and has been the subject of intense research in the past [13]. In contrast to d = 1, all analytic and numerical evidence supports a groundstate displaying long-range antiferromagnetic order in this case.
Let us now discuss the case of non-zero Wilson parameter 0 < r < 1, where Eq. (27) leads to directional Heisenberg-Ising anisotropies |J 1,y | > |J 1,x | = |J 1,z | and |J 2,x | > |J 2,y | = |J 2,z |. In the limit a 1 /a 2 → 0 (a 1 /a 2 → ∞), we again have a collection of uncoupled spin chains along the rows (columns), as discussed for the quantum compass model (19), with the difference that each of these rows (columns) is no longer described by a quantum Ising model (19) but, instead, by an XYZ chain [110,111]. As remarked around Eq. (25), we find an even number of negative spin-spin couplings in these rows (columns), such that a specific spin rotation maps this model to an antiferromagnetic XXZ model where the Ising anisotropy is always larger than unity. In analogy to the d = 1 case discussed above, we also expect to find either a ferromagnetic FM x ( FM y ) ordering for each of the columns (rows), or, otherwise, a symmetric paramagnet (PM) when the transverse field is larger than the leading spin-spin coupling. Also following this analogy, as we depart from the limits of large anisotropies a 1 /a 2 → 0 (a 1 /a 2 → ∞), the additional spin-spin couplings will lock these columns (rows) to the same ferromagnetic order, and we expect that the critical points of the purely 90 o compass model |h z | < |J 2,x |/ζ and |J 2,x | > |J 1,y | for the FM x -PM transition, and |h z | < |J 1,y |/ζ and |J 1,y | > |J 2,x | for FM y -PM transition will change, such that ζ flows with the Wilson parameter. Since these Heisenberg-Ising compass models are no longer solvable, this can only be determined numerically, or using certain approximations that are discussed below. We will start by deriving a path-integral representation of these effective spin models that will connect to variants of the constrained QFTs discussed previously.
In this section, we derive a path-integral representation for the partition function Z = Tr{e −β H eff } of the strong-coupling Heisenberg-Ising model (13), which will allow us to understand how constraints in QFTs such as Eq. (33) below can appear in our fermionic model (8). This requires the use of spin coherent states [112], which can be defined by performing a SU(2) rotation on a fiducial state, which ful- The coherent-state basis, depicted in Fig. 3, can thus be defined by the action of the following operator on the tensor product of fiducial states for each lattice site Here, we have introduced polar θ (x) and azimuthal φ (x) angles, which define a unit vector per spin pointing along the radial outward direction of a unit 2-sphere S 2 . Therefore, |ω ω ω(τ, x x x)| 2 = 1 parametrises all possible spin directions ω ω ω(x) = sin θ(x) cos φ(x)e e e x + sin θ(x) sin φ(x)e e e y + cos θ(x)e e e z .
(30) Note that x = (τ, x x x) now represents the Wick-rotated spacetime points, where imaginary time τ extent is related to inverse temperature via τ ∈ [0, β ] [113]. One readily checks that in this basis S S S(x x x) = Sω ω ω(x x x), such that the components of this unit vector field contain information about the fermionic σ and π condensates mentioned above σ (x) = ω z (x), π π π(x) = ω x (x)e e e x + ω y (x)e e e y .
One can now rewrite the partition function as a path integral where the vector field Ω(τ, x x x) is constrained to lie on the S 2 sphere through the specific form of the integral measure. As discussed in [13,96,113], the Euclidean action contains a geometric contribution proportional to the sum of the Berry phases of each spin history as it moves along the corresponding trajectory Γ x x x : τ → Ω Ω Ω(τ, x x x) on its respective sphere, these trajectories being closed due to the periodic boundary conditions along the τ direction Ω Ω Ω(0, x x x) = Ω Ω Ω(β , x x x). Altogether, the action is expressed as where the first contribution corresponds to the aforementioned Berry phase [114], and is known as a Wess-Zumino term. For each spin, this term can be understood as an effective Aharonov-Bohm phase gained by a unit test charge q e = 1 moving on the sphere, and subjected to the magnetic field of a monopole of charge q m = 4πS located at its centre [74]. Using Stokes' theorem, this phase can be rewritten as the magnetic flux across the spherical cap enclosed by each spin trajectory containing the north pole of S 2 , where the fiducial state |S, +S x x x points to. Hence, the effective vector potential and magnetic field are those generated by the magnetic monopole e e e z × ω ω ω(τ, x x x) 1 + e e e z · ω ω ω(τ, x x x) , The second term in (34) represents the additional coupling of neighbouring spins due to the spin-spin couplings, as well as their precession under the transverse field. Let us briefly discuss the r → 0 naive-fermion limit, and its relation to an O(3) non-linear sigma model [96]. In this limit, the spin-spin couplings along the different internal directions are all equal |J j,a | = J j , and there is a continuous SU(2) symmetry. In light of the field constraint in the integral measure (32), and up to an irrelevant constant term, the nearest-neighbour couplings can be rewritten as follows −J j a j 1 2a j (Ω Ω Ω(τ, x x x + a j e e e j ) − Ω Ω Ω(τ, x x x)) 2 → J j a j 2 ∂ j Ω Ω Ω · ∂ j Ω Ω Ω, which clearly resembles the spatial derivative-terms in Eq. (1), and can be understood as the energy contribution due to the strain caused by a field deformation. The kinetic part, which depends on the canonical momenta of the scalar vector field in Eq. (1), would appear in the form of Euclidean time derivatives in the corresponding action, and is not readily apparent in equation (34). However, a specific parametrisation of the spin trajectory shows that the Berry phase indeed contains these time derivatives β 0 dτA A A(Ω Ω Ω) · ∂ τ Ω Ω Ω, albeit still being different from the kinetic terms of the O(3) non-linear sigma model. In a seminal work [96], F.D.M. Haldane showed that for antiferromagnetic Heisenberg couplings J j = J > 0, an expansion of Ω Ω Ω(τ, x x x) about the saddle point of the Euclidean action, corresponding to the alternating antiferromagnetic configuration of the classical limit, exactly yields the kinetic term of the O(3) non-linear sigma model. Moreover, in d = 1, the Berry phase also contributes with a topological theta term θ = 2πS which, depending on the half-integer or integer value of the spin S, makes the O(3) non-linear sigma model massless or massive [115].
Let us now explore how this situation changes for our current spin models. Since the effective spin models do not have the continuous SU(2) symmetry of the Heisenberg limit for generic 0 < r < 1, we first need to understand the nature of the saddle point of Eq. (34) controlling the large-S limit, which will eventually lead to a different type of constrained non-linear sigma model. By inspection of the action (34) and the magnetic monopole fields (35), one can readily see that the dΩ Ω Ω · A A A(Ω Ω Ω(τ, x x x)) term can only depend on the equatorial components of the scalar field via ∑ x x x∈Λ d (Ω x ∂ τ Ω y − Ω y ∂ τ Ω x ). This term remains invariant under the Z 2 symmetry (15) which in this context reads Accordingly, the action (34) with the constraint (32) (equivalent to the non-linear sigma constraint (33) via Eq. (31)) describes a Z 2 non-linear sigma model that arises naturally in the strong-coupling limit of the original model (8), even for a single fermion flavour N = 1.

IV. LARGE-S LIMIT AND SADDLE-POINT EQUATIONS
Let us recall that the spin operators in Eq. (12) correspond to S = 1/2. In this section, nonetheless, we will assume that S is a free parameter, and explore the S → ∞ limit, which can be understood as a mean-field approximation of the effective spin chain. According to our previous discussion of the Berry phase, one can see that time-dependent spin histories ∂ τ Ω Ω Ω = 0 0 0 will get suppressed in this limit, due to the averaging of the rapid oscillations associated to the pure-imaginary Wess-Zumino term. Accordingly, the large-S limit will be controlled by static fields Ω Ω Ω(τ, x x x) = Ω Ω Ω(x x x). This is precisely analogous to the large-N limit in interacting fermion theories; e.g. in the d = 2 Thirring model the induced Chern-Simons term resulting from the leading quantum correction [116] plays no role in determining the ground state in the large-N limit. Moreover, the Euclidean action can be rewritten as S E = 1 h eff s E wherē h eff ∝ 1/S plays the role of an effective Planck constant. In the absence of dynamics and kinetic terms, the Euclidean action per spin s E = βV eff can be expressed in terms of the following effective potential where we note that, in analogy to the large-N limit [7] of the four-Fermi term (5), the spin-spin couplings (22)-(27) must be rescaled so as to give finite contributions for S → ∞ whereJ j,a are finite and non-zero coupling strengths. Accordingly, the large-S limit is controlled by the saddle point of the potential (40) given below, in which quantum fluctuations are suppressedh eff → 0, bearing in mind that the fields are subjected to an extensive number of constraints, i.e. one per spatial coordinate, as they must lie on their corresponding S 2 spheres (33). At this level, one can either introduce a Lagrange multiplier to deal with the non-linear constraint, or work directly with the generalised 'coordinates' {Ω Ω Ω(x x x)} → {θ (x x x), φ (x x x)} satisfying the constraints (30). In this second approach, the saddle-point equations are given by a set of 2 ∏ j N j non-linear equations ∀x x x ∈ Λ d , namely A. Large-S Ising magnetism for d = 1 Let us start by discussing the solutions of these saddle-point equations for d = 1. Following the discussion in Sec. II B, where we argued that the spin couplings (22) for 0 < r < 1 can be mapped onto an antiferromagnetic Heisenberg-Ising chain (25), we can simplify the set of non-linear equations (39) by restricting to translationally-invariant configurations within a 2-site unit cell {θ (x x x), φ (x x x)} → {θ s , φ s }, which may capture a possible alternating order, where s ∈ {A, B} stands for the odd/even sites of the chain (see Fig. 3(a)). Under this simplification, the effective potential reads which leads to a system of four non-linear equations via (39).
To gain some knowledge about the groundstate, one can numerically search for a global minimum of this potential using a coarse-grained discretization of the angles, and performing a grid search restricting the search space to account for the Z 2 symmetry. Once we have a rough estimate of the minimum, for further accuracy we directly minimise the effective potential as a non-linear constrained problem using an interior-point algorithm of non-linear programming [117], choosing as initial points the outcomes of the global coarser minimisation. In practice, we have also initialised the minimisation in all the different combinations of cardinal states of the two S 2 spheres associated to the 2-site unit cell, and checked for consistency, ensuring that the solution found with the nonlinear programming algorithm that starts with the grid-search minimum yields the minimum potential among the different  Fig. 4 (a), we use these numerical large-S solutions to represent the corresponding SSB order parameter, which corresponds to the pseudo-scalar condensate S y (x x x) ∝ Π 5 in Eq. (17). This figure clearly depicts a SSB region hosting an Ising ferromagnet FM y , corresponding to the paritybreaking Aoki phase with a non-zero pseudo-scalar condensate. This phase is separated from the one invariant under parity, namely a disordered paramagnet PM in the language of the spin model, via a critical point where the order parameter behaves non-analytically. To find the accurate location of this point, we represent in Fig. 4 (b) the corresponding susceptibilities χ y = ∂ S y (x x x) /∂ h z , which clearly peak at the corresponding points. These figures show, as qualitatively argued in the previous section,that the critical point h z /J 1,y c flows with the value of the Wilson parameter r. For r = 1, the critical point obtained from the numerical minimisation is |h z /J 1,z | c ≈ 1. This coincides with the analytical solution of the saddle-point equations which, in this limit, can be found We thus find that, for |h z /J 1,y | < |h z /J 1,y | c = 1, the spins align according to where the two possible signs ± account for the two-fold degeneracy associated to the Z 2 parity SSB.
As one now varies 0 < r < 1, it is simple to understand the flow of the critical point by performing a self-consistent mean-field decoupling of the effective Hamiltonian (13). In light of the groundstate expectation values (42), a mean-field decoupling of the additional terms of the Heisenberg-Ising chain (25) that arise when r = 1 (43) would only contribute with terms along the internal z direction, effectively shifting the transverse field to h z →h z = h z + 2J 1,z S z (x x x) . The saddle-point solution (42) now yields a self-consistent equation, which can be readily solved for the couplings in Eq. (22): In order to test the validity of this prediction, we present a contour plot of the Ising order parameter in Fig. 4 (c) as a function of the relative coupling strengths |J 1,x /J 1,y | = |J 1,z /J 1,y | = (1 − r 2 )/(1 + r 2 ) and |h z /J 1,y |. We also plot the critical points extracted from the numerical maxima of the susceptibility in Fig. 5, depicted with red stars, and their analytical prediction (44), represented with a white dashed line. The agreement is very good, and identifies the main source for the flow of the critical coupling with the Wilson parameter. In the d = 1 case, the critical point flows towards smaller values |h z /J 1,y | c < 1 due to the effective renormalisation of the transverse field h z →h z , which increasesh z > h z due to the nonvanishing ferromagnetic couplings J 1,z < 0 along the internal z axis, and the specific alignment of the spins in Eq. (42). Coming back to the context of four-Fermi-Wilson lattice field theories (8), we see that the Aoki phase extends for arbitrary values of the Wilson parameter 0 < r < 1, provided that one goes sufficiently close to the so-called central Wilson branch m = −1/a 1 [118,119]. For strictly vanishing r = 0, where we recover the naive-fermion discretization, the Aoki phase is no longer present. Right at the central branch, where h z = 0, the strong-coupling groundstate lies exactly at the critical line, and would correspond to the gapless groundstate of the antiferromagnetic Heisenberg chain. Before moving towards the d = 2 case, let us note that the large-S approximation is not expected to provide an accurate estimate of the exact critical point but, at least, it captures  Fig. 4(b), as well as the white dashed line for the predictions in Eq. (44), which correspond to a straight line with negative unit slope h z /J 1,y c = 1 − |J 1,x /J 1,y |.
qualitatively the main sources for the flow of the critical point. In fact, for the quantum Ising model at r = 1, where large-S predicts a critical point |h z /J 1,z | = 1 in Fig. 4(a)-(b), the exact solution gives instead instead |h z /J 1,z | = 1/2 [85]. In section V, we will present more accurate predictions of the critical points using the quasi-exact DMRG algorithm based on matrix product states.

B. Large-S compass magnetism for d = 2
Let us now discuss the d = 2 case, where the A and B sub-lattices correspond to two interpenetrating square lattices with lattice spacing (a 2 1 + a 2 2 ) 1/2 , rotated with respect to the original rectangular lattice by angles ϕ, π − ϕ, where ϕ = arctan(a 2 /a 1 ) (see Fig. 3(b)). In analogy to d = 1, we could restrict the configurations of spin coherent states to be translationally-invariant within these sub-lattices where s ∈ {A, B}, and still account for anti-ferromagnetic configurations when the respective A, B angles differ, or for ferromagnetic ones when they are equal. However, this choice may not suffice to capture the groundstate ordering of the effective strong-coupling model (13). For illustrative purposes, consider the limit a 1 a 2 and r ≈ 1 so that, in light of the spin-spin couplings in Eq. (27), |J 2,x | {|J j,a |, ∀ j = 2, a = x}. Since this leading spin coupling is negative J 2,x < 0, the spins will want to align ferromagnetically along each of the columns, adopting polar and azimuthal an- Additionally, since the perturbative coupling between neighbouring columns fulfils J 1,x > 0, spins in adjacent columns might minimise the groundstate energy by choosing opposing azimuthal angles φ * A = φ * B ∈ {π, 0}. Since this is inconsistent with the A, B sub-lattice layout, allowing for this possible ordering in the parametrisation of the constrained effective potential requires augmenting the number of configurations by considering a 4-site unit cell {θ (x x x), φ (x x x)} → {θ s,c , φ s,c }, where c ∈ {l, r} labels the left and right corners, as depicted in Fig. 3(b). By directly incorporating the non-linear constraint as we did forEq. (40), the effective potential reads in this case where we have introduced the functionc(2) = c when the spin-spin couplings occur along the e e e 2 spatial direction, and c(1) = c along e e e 1 . For the latter, we define c = r(l) for c = l(r), which swaps the left and right corners. The saddle-point conditions (39) corresponding to this potential lead to a non-linear system of 8 equations which, once again, must be solved numerically for generic cases. The exception is the standard limit r = 1, where the effective spin model reduces to the 90 o compass model [88] in a transverse field, and one finds The SSB order parameters for these solutions are which predict critical points at |h z /J 1,y | c = 1 if a 1 < a 2 , and |h z /J 2,x | c = 1 if a 1 > a 2 . Accordingly, for |h z | < |J 2,x | and a larger horizontal lattice spacing a 1 > a 2 , the SSB order parameter S x (x x x) ∝ Π 1 corresponds to ferromagnetic ordering along the internal x axis FM x , which corresponds to the inversion-breaking π condensate of Eq. (20) for the underlying four-Fermi model. Alternatively, for |h z | < |J 1,y | and a larger vertical lattice spacing a 1 < a 2 , the SSB order parameter describes ferromagnetic ordering along the internal y axis FM y , which corresponds to the other inversion-breaking π condensate S y (x x x) ∝ Π 2 . We note that these large-S solutions for r = 1 coincide with the variational mean-field estimates discussed in [60,61], and recall again that these condensates are different from the Aoki parity-breaking phase. To treat 0 < r < 1, we must solve the problem numerically. We use the same strategy as described for d = 1, which combines a coarse global minimisation with more efficient non-linear programming methods that are consistently initialised to yield accurate estimates of the potential minima. We obtain the SSB order parameters from the numerical saddle points {θ s,c , φ s,c } for various Wilson parameters r ∈ {0.5, 0.6, 0.7, 0.8, 0.9, 1}. The corresponding magnetisations display similar non-analytic behaviours, which can be used to infer the location of the critical points, and how these flow as one varies r. In Fig. 6, we present a stack of twodimensional contour plots that summarises the large-S phase diagram, and shows a clear dependence on both the Wilson parameter and the anisotropy parameter In this contour plot, we represent the difference of the two possible SSB order parameters S y (x x x) − S x (x x x) , such that negative (positive) values signal a FM x (FM y ) phase with a Π 1 (Π 2 ) Lorentz-breaking condensate, and are depicted in red (yellow) scale. In the stacking z direction, we plot the anisotropy parameter for ξ 2 ∈ [0, 1] for a 1 < a 2 (black axis), while we represent 2 − 1/ξ 2 ∈ [1, 0] for a 2 > a 1 (grey axis). The lower stacked contour plots thus represent FM y ordering, whereas the upper ones represent FM x . In the x and y axes of these contour plots, we select the relevant normalised couplings such that the stacked contour plots are completely symmetric as one crosses the isotropic configuration a 1 = a 2 . In general, one observes that the critical point |h z /J 2,x | c = 1 (|h z /J 1,y | c = 1) at r = 1 and ξ 2 > 1 (ξ 2 < 1), changes as one decreases r, increasing in this way the remaining coupling strengths of the compass Heisenberg model (27), namely |J 2,z /J 2,x | = |J 2,y /J 2,x | ( |J 1,z /J 1,y | = |J 1,x /J 1,y |). In addition, one can also observe how the critical points change with anisotropy when r < 1, such that the extent of the Lorentz-breaking fermion condensates in general depends on the anisotropy of the lattice regularisation. To gain some further understanding, we can extend the previous discussion of the d = 1 case around Eq. (43) to cover d = 2 in the regime of very large anisotropy, since in this case the compass Heisenberg-Ising models reduce to very weakly-coupled columns (rows). To derive a self-consistent mean-field decoupling for d = 2, note that we now have to deal with these additional spin-spin couplings between adjacent rows (columns) for a 1 > a 2 (a 1 < a 2 ). Setting r ≈ 1, one can solve the corresponding self-consistent equations approximately, leading to The white dashed lines of Fig. 6 represent these analytical predictions of the critical points. As can be observed, they match nicely the numerical critical lines for large anisotropies ξ 2 1 (ξ 2 1), but the discrepancy increases as one approaches the isotropic case ξ 2 = 1.

V. TENSOR-NETWORK NUMERICAL SIMULATIONS
In this section, we test the validity of the large-S predictions of Sec. II by means of variational algorithms based on tensor network states (TNSs) [120][121][122].
The quantum state of a lattice model composed of N sdlevel systems, i.e. spins, can be written in the basis of tensor products of local states. The quantum state is then fully characterised by the coefficients of these basis states, which are tensors C i 1 ,i 2 ,··· ,i N of rank N s and dimensiond Here, we have introduced the indexes i n ∈ {1, · · · ,d}, such that the description of the state requires ofd N s complex parameters. This exponential growth makes this generic description unsuitable for numerical analysis. However, in a number of situations, the physically-relevant many-body states admit a more concise description based on TNSs [123,124]. Obtained from a contraction of low-rank tensors on so-called virtual indices, TNSs economically approximate the states with local interactions in thermal equilibrium. The number of required parameters scales only polynomially with system size [125], circumventing the previous exponential growth of the most generic description. In fact, these variational states are based on the powerful insights related to the area law [126,127]. The area law places bounds on quantum entanglement that a many-body system can generate, which translates directly to the number of parameters required to describe a physically-relevant quantum state. Tensor-network calculations benefited from the advent of White's density matrix renormalisation group (DMRG) [128], famous for its extraordinary accuracy in solving onedimensional quantum systems, which is intimately connected with a tensor decomposition known as the matrix product state (MPS) [129,130]. These variational states can be understood in terms of pairs of maximally-entangled states on neighbouring lattice sites, which describe auxiliary degrees of freedom, and get locally projected onto the lower-dimensional subspace of physical spins at each lattice site. In fact, a very useful and intuitive way of thinking about MPS is the following valence-bond construction. Consider N s spins aligned on a ring, the states of which are labelled by the internal index i. One assigns two auxiliary spins of dimension D to each of these physical spins, assuming that each pair of neighbouring auxiliary spins is initially in a maximally entangled state |I = ∑ D α=1 |α, α , often referred to as an entangled bond. Applying the map that plays the role of the aforementioned projector to each of the N s spins, and interpreting A i as a D ⊗ D matrix, we find that the coefficients of the final state can be expressed by a matrix product Tr [A i 1 A i 2 · · · A i N ] (see Fig. 7 (a)). In the d = 1-dimensional models discussed in this work, we consider physical spins S = 1/2, such that i n = s x x x n ∈ {↑, ↓},d = 2, and n ∈ {1, · · · , N s } with N s = N 1 being the number of sites of the spatial chain. In general, the dimension of the entangled state |I can be site-dependent and we write A [x x x n ] s x x xn for the D n × D n+1 matrix corresponding to site x x x n ; the states then have the form (52) and are called MPS. This construction can be mathematically expressed as a network of tensors with multiple indexes corresponding to the physical and auxiliary degrees of freedom, such that those corresponding to the auxiliary ones are contracted as described in Fig. 7 (b). In this case the number of parameters needed to describe a physical state in the MPS language scales as O(N 1 dD 2 ) with d the physical dimension of the spins.
A natural generalization of MPS to two, or even higher, spatial dimensions is represented by projected entangled pair states (PEPS [131]). Again, this kind of state can be understood in terms of pairs of maximally-entangled states of neighbouring auxiliary systems, which are are locally projected into the low-dimensional physical subspace. As represented in Fig. 7 (c), the PEPS describes a state through interconnected tensors. For the two-dimensional spatial lattices considered in this work, which consist of N s = N 1 N 2 sites, we specify the PEPS variational ansatz [131][132][133] as s x x xn , some of which are connected according to the geometry of the lattice and the notion of neighbouring lattice sites. Each tensor of the PEPS has N b so-called bond indices of dimension D n , which describe the aforementioned auxiliary degrees of freedom, and a single physical index of dimension d. The choice of N b in the tensor network can be arbitrary, and typically depends on the geometry of the lattice. For example for a N 1 × N 2 lattice, a PEPS contains N 1 × N 2 bulk tensors with N b = 4 and D n = D. Each tensor depends on dD 4 complex coefficients. Therefore the PEPS is characterized by O(N s dD 4 ) parameters. The function F contracts all the tensors A s x x xm , according to this pattern, and then performs the trace to obtain a scalar quantity such that Eq. (53) can be understood as a parametrisation of a particular set of states in the exponentially-large physical Hilbert space.
In this manuscript, we study the groundstate properties of quantum lattice Hamiltonians using different strategies for d = 1 and d = 2 spatial dimensions. For d = 1 we variationally optimize the MPS tensors, so as to minimise the expectation value of the corresponding Hamiltonian. By contrast for d = 2, in analogy to spectroscopic methods that determine the particle spectrum via the imaginary-time evolution of correlators in Euclidean LFTs [66], we evolve the system in imaginary time until a stationary state corresponding to the groundstate is reached. This assumes that this groundstate is unique, and that the energy gap is non-zero, as done in the time-evolving block-decimation method (TEBD) for onedimensional chains [134,135]. In the following, we will use this method in the thermodynamic limit for the infinite PEPS state (iPEPS) [136,137].
A. Tensor-network Ising magnetism for d = 1 In this section, we analyse the effect of correlations in the phase diagram of the Heisenberg-Ising chain (13) with spinspin couplings defined in Eq. (22), and subjected to an additional transverse field in Eq. (21). All of these parameters depend on the Wilson r, and the goal is to explore the phase diagram as it is varied within 0 < r < 1. In particular, we benchmark the large-S results discussed in Sec. IV A, giving more accurate predictions of the phase diagram and critical points presented in Fig. 4 (c).
As discussed in Sec. IV A, the Heisenberg-Ising chain presents a critical line separating the ferromagnet FM y and the paramagnet PM. In Fig. 8 we present the corresponding MPS phase diagram as a function of relative coupling strengths |J 1,x /J 1,y | = |J 1,z /J 1,y | and |h z /J 1,y |. Our numerical results for the phases of matter are extrapolated using the quasi-exact DMRG algorithm, as discussed in detail below. The lines represent the critical points where the SSB phase transitions occur, either obtained with DMRG based on finite MPS with bond dimension D = 200 (red stars), or by self-consistent mean-field method (yellow dashed lines), which exploits exact solutions of the transverse-field Ising model to derive a selfconsistent equation for the transverse magnetization that can be solved analytically in the limit |J 1,x /J 1,y | = |J 1,z /J 1,y | 1. The ferromagnetic susceptibility χ M y , for fixed coupling strength J y = −1, and for different couplings J x . As magnetic field h z is varied it develops peaks at the critical points. In the inset, we show ferromagnetic magnetization along the y direction. The system develops a non-zero expectation value for transverse fields below a critical value h z < h c z . (b) The paramagnetic susceptibility χ M z , for the same parameters, which develops peaks at those critical points. In the inset, we show ferromagnetic magnetisation along the z direction.

This yields
which must be compared to the large-S estimate (54), depicted with a white dashed line in the figure.
As can be observed in Fig. 8, for small relative coupling |J 1,x /J 1,y |, attained for r ≈ 1, the self-consistent mean-field and DMRG critical points separating the FM y and PM regions yield two critical lines that are very similar. Increasing |J 1,x /J 1,y |, larger differences appear between the critical lines, since the self-consistent mean-field predicts a smaller FM y region. Note, however, that this prediction is strictly valid only in the vicinity of the exact critical point |h z /J 1,y | c = 1/2. In comparison with the large-S white-dashed line, we see that there is a large quantitative difference of the critical lines, e.g. for r = 1 |h z /J 1,y | c = 1 within large-S, whereas |h z /J 1,y | c = 1/2 with DMRG. This difference is characteristic of large-S methods, and is a consequence of how we neglect quantum fluctuations by taking the saddle-point solution. However, note that the large-S captures the physics of the model qualitatively correctly: it predicts that the FM y region (pseudoscalar condensate) shrinks proportionally to the combination (1 − r 2 )/(1 + r 2 ), whenever 0 < r < 1. Had we solved the lattice model for increasing S using DMRG, we would have found that the two lines approach each other as the spin is increased S ∈ {1/2, 3/2, 5/2, · · · }.
In the background of Fig. 8, we present a contour plot of the SSB order parameter, which corresponds to the pseudo-scalar condensate S y (x x x) ∝ Π 5 defined in Eq. (17). In order to avoid numerical problems due to the incomplete symmetry breaking of the magnetization M y = S y (x x x) = 1 N 1 ∑ x x x∈Λ 1 S y (x x x) , we determine instead the corresponding structure factors The zero-momentum component of these structure factors yield the desired magnetisation in the thermodynamic limit M y = (S yy (0)) 1/2 . The contour plot of the magnetisation clearly identifies the SSB region in yellow-green scale with a non zero pseudo-scalar condensate, namely a Ising ferromagnet FM y , which is separated from a region in blue scale, where the parity is preserved, namely a paramagnet PM with zero magnetization along the internal y axis.
Let us now give some more details on the methodology used to extract numerically the critical points shown in Fig. 8. By calculating the ferromagnetic magnetisations and, particularly, the corresponding susceptibilities, we can identify the critical points occurring at a non-zero external field h z > 0. In Fig. 9 (a), we present the susceptibility χ y = ∂ M y /∂ h z for different values of Wilson parameter r, which clearly peaks at specific values of the ratio |h z /J 1,y | that move to the left as r is decreased. In Fig. 9 (b), we display the ferromagnetic susceptibility χ y for different number of sites N s , fixing r = 0.82, and varying the external magnetic field h z . The finite size scaling (FSS) of the magnetic susceptibility maxima, as a function of N 1 is displayed in the lower inset. As one can see in the inset, the peak of the chiral susceptibility at transverse field h c z diverges with the size of the chain, and fitting the maxima of h c z to h c z (N 1 ) = h c z (1 + aN −1 1 + bN −2 1 ), we can delimit the ferromagnetic region and locate the phase transitions in the thermodynamic limit N s → ∞. Once the critical point is known, in the upper inset we show the data collapse of the magnetization curves when rescaled with the system size using the critical exponent of the 2D Ising universality class. Accordingly, the whole critical line delimiting the Aoki phase belongs to this universality class, in spite of having perturbation to the Ising limit in the form of the additional spin-spin couplings (22) when r < 1. In this section, we show the results obtained by using the above iPEPS algorithm for the Heisenberg-Ising compass model (13) with spin-spin couplings defined in Eq. (27), and subject to an additional transverse field in Eq. (26), working directly in the infinite-lattice limit. In particular, we have computed the ground state wave function |ψ GS of the system by performing the imaginary-time evolution for different values of the spin couplings {J 1,a , J 2,a } a∈{x,y,z} , and the transverse magnetic field h z , and then evaluated observable quantities on it, such as the groundstate energy and the local order parameters related to the ferromagnetic phases.
In Sec. IV B, we used a large-S method to predict a critical line separating the symmetry-broken ferromagnets FM x and FM y from a paramagnet via second-order phase transitions. Under certain approximations, these critical lines can be analytically found (49), corresponding to the white dashed lines of Fig. 6. In order to test the validity of these large-S predictions, we use our iPEPS algorithm for D = 2. By measuring paramagnetic and ferromagnetic magnetisations, we confirm that these quantities can be used to identify the critical points also for a non-zero magnetic field h z > 0.
We start by setting the spatial anisotropy to ξ 2 = 3.16, which would correspond to a specific stacked plane in the large-S phase diagram of Fig. 6, where FM x order competes with the disordered PM. Our numerical iPEPs results for this competition are presented in Fig. 10. The lines correspond to the critical points where the FM x and PM phase transitions occur, either obtained with a imaginary time evolution based on infinite iPEPS with bond dimension D = 2 (red stars), or by the large-S approximate prediction (49) (orange dashed lines). We observe a similar trend as in the case of d = 1; the region of the inversion-breaking condensate shrinks as the value of the Wilson parameter r is reduced within 0 < r < 1. In the vicinity of the standard choice r = 1, and for ξ 2 1, we see once again that the region of non-vanishing condensate decreases proportionally to the ratio (1 − r 2 )/(1 + r 2 )) = |J 2,y /J 2,x |. On the other hand, as r → 0, the iPEPs critical line bends upwards, as also occurred for d = 1 (see the red stars in Fig. 8). Note that, although the large-S and iPEPs critical lines cross for smaller values of r, i.e. larger ratios of |J 2,y /J 2,x |, the analytical predictions (49) are not strictly valid in this regime. On the other hand, for the regime where r ≈ 1, we see that the large-S predictions are closer to the iPEPs results in comparison to the d = 1 case, showcasing that mean-field predictions typically improve as dimensionality increases.
Let us again discuss details on how we extracted such critical points numerically. In the inset of Fig. 11, we present the magnetisation M x = S x (x x x) as a function of transverse magnetic field h z , setting J 1,x = −1 and exploring different values of J 1,y < 1. This figure shows that, for weak transverse fields, the magnetization attains a non-zero value signalling the broken symmetry FM x phase, which corresponds to the inversion symmetry-broken fermion condensate. The main panel shows the corresponding magnetic susceptibility χ x peaking at a specific value of the transverse field, which can be used to locate the corresponding critical points. Note that these peaks are not as pronounced as for d = 1 and that, given that we work with translationally-invariant iPEPS, we cannot perform FSS to see how the peak diverges and extract accurate estimations of the critical point and universality class. In future studies, it would be interesting to push the numerics to explore larger values of bond dimension D > 2, which would permit a more accurate location of the critical points. This question is particularly relevant in the regime r → 0. where d = 1 and d = 2 results seem to differ qualitatively. In the 1D case, the Aoki phase shrinks all the way to zero, whereas in the 2D case it seems to survive. This could be related to the fact that the r = 0 limit maps onto a Heisenberg model on a rectangular lattice, and that this model has long-range order in contrast to the 1D version. We note that these questions could also be addressed with other methods such as Monte Carlo simulation.

VI. CONCLUSIONS AND OUTLOOK
In this work, we have explored the limit of strong Hubbard interactions in models of correlated topological insulators that arise for spin-orbit coupled fermions in lattices of one and two spatial dimensions. These models can be understood as the single-flavour limit N = 1 of four-Fermi quantum field theories with a Wilson-type discretization, making an interesting and fruitful connection between condensed matter and lattice field theories. As discussed in this work, most lattice field theory studies fix the Wilson parameter of this discretization to r = 1, as a non-unity value has trivial consequences in the non-interacting limit. However, understanding the role of 0 < r < 1 in the presence of interactions is not clear a priori, which is the question explored in this work. Moreover, given the fact that these four-Fermi field theories are amenable to study using cold-atom quantum simulators, with spin-orbit coupled fermions in Raman lattices where the effective value of r depends on intensities of the lasers that control the Raman lattice, which are generically different, i.e. r = 1, the question addressed in this work is also relevant for understanding the possible phases that can be explored in possible cold-atom realisations.
To address this question, we have derived an effective spin model in the limit of strong four-Fermi interactions, finding a specific dependence of the couplings on the Wilson parameter r. In d = 1, the resulting model can be related to an XXZ, also known as Heisenberg-Ising, model in a staggered magnetic field, whereas in d = 2 it is related to a compass model with directional spin-spin couplings, each of which is described by different Heisenberg-Ising couplings. We have formulated a path-integral representation of the partition function, which connects the strongly-interacting limit with a constrained QFT: a non-linear sigma model with a discrete Z 2 symmetry. This permits exact solutions in the large-S limit, which enable us to identify the relevant phases of matter, and draw specific predictions about phase transitions and the flow of the critical points with the Wilson parameter r. The validity of these predictions have been tested against tensor-network numerical simulations, which show that the large-S diagrams are qualitatively correct. On the other hand, the numerical results give more accurate estimates of the flow of the critical lines and, in some cases, allow us to infer the correct universality class of the lines that differs from the large-S meanfield-type scaling.
As an outlook, we believe that the present manuscript, together with other recent works [57][58][59][60][61], provide a rich crossdisciplinary toolbox to understand interaction effects in topological matter. It will be interesting to exploit, adapt, and combine all of these different techniques to study harder problems such as, for instance, abandoning half-filling and exploring correlated topological phases at non-zero fermion densities.